LeukFr: Cox regression with random effects

Freireich et al (1963)'s data presented in the Leuk example actually arise via a paired design. Patients were matched according to their remission status (partial or complete). One patient from each pair received the drug 6-MP whilst the other received the placebo. We may introduce an additional vector (called pair) in the BUGS data file to indicate each of the 21 pairs of patients.

We model the potential 'clustering' of failure times within pairs of patients by introducing a group-specific random effect or frailty term into the proportional hazards model. Using the counting process notation introduced in the Leuk example, this gives


   I_i (t) dt   =   Y_i (t) exp( β' z_i + b_pairi ) dΛ₀(t)   i = 1,...,42;   pair_i = 1,...,21
   b_pairi   ~   Normal(0, τ)

A non-informative Gamma prior is assumed for τ, the precision of the frailty parameters. Note that the above 'additive' formualtion of the frailty model is equivalent to assuming multiplicative frailties with a log-Normal population distibution. Clayton (1991) discusses the Cox proportional hazards model with multiplicative frailties, but assumes a Gamma population distribution.

The modified BUGS code needed to include a fraility term in the Leuk example is shown below

   model
   {
   # Set up data
   for(i in 1 : N) {
      for(j in 1 : T) {
# risk set = 1 if obs.t >= t
         Y[i, j] <- step(obs.t[i] - t[j] + eps)
# counting process jump = 1 if obs.t in [ t[j], t[j+1] )
# i.e. if t[j] <= obs.t < t[j+1]
         dN[i, j] <- Y[i, j ] *step(t[j+1] - obs.t[i] - eps)*fail[i]
      }
      }
   # Model
      for(j in 1 : T) {
         for(i in 1 : N) {
            dN[i, j] ~ dpois(Idt[i, j])
            Idt[i, j] <- Y[i, j] * exp(beta * Z[i]+b[pair[i]]) * dL0[j]
         }
         dL0[j] ~ dgamma(mu[j], c)
         mu[j] <- dL0.star[j] * c # prior mean hazard
   # Survivor function = exp(-Integral{l0(u)du})^exp(beta * z)
         S.treat[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * -0.5))
         S.placebo[j] <- pow(exp(-sum(dL0[1 : j])), exp(beta * 0.5))   
      }
      for(k in 1 : Npairs) {
         b[k] ~ dnorm(0.0, tau);
      }
      tau ~ dgamma(0.001, 0.001)
      sigma <- sqrt(1 / tau)
      c <- 0.001 r <- 0.1
      for (j in 1 : T) {
         dL0.star[j] <- r * (t[j+1]-t[j])
      }
      beta ~ dnorm(0.0,0.000001)
   }

   Results