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
Λ0(t) i = 1,...,42;
pairi = 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)
}
Data
list(N = 42, T = 17, eps = 0.00001, Npairs = 21,
t = c(1,2,3,4,5,6,7,8,10,11,12,13,15,16,17,22,23,35),
obs.t = c(1,1,2,2,3,4,4,5,5,8,8,8,8,11,11,12,12,15,17,22,23,
6,6,6,6,7,9,10,10,11,13,16,17,19,20,22,23,25,32,32,34,35),
pair = c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,
19,18,8,1,20,6,2,10,3,14,4,11,7,9,12,16,17,5,13,15,21),
fail = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,
1,1,1,0,1,0,1,0,0,1,1,0,0,0,1,1,0,0,0,0,0),
Z = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5,
0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5,
-0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5)
)
Inits for chain 1
list( beta = 0.0,
dL0 = c(1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,
1.0,1.0,1.0,1.0,1.0,1.0, 1.0,1.0),
tau=1
)
Inits for chain 2
list( beta = 1.0,
dL0 = c(2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,
2.0,2.0,2.0,2.0,2.0,2.0, 2.0,2.0),
tau = 0.1)
Results