Kidney: Weibull regression with random efects
McGilchrist and Aisbett (1991) analyse time to first and second recurrence of infection in kidney patients on dialysis using a Cox model with a multiplicative frailty parameter for each individual. The risk variables considered are age, sex and underlying disease (coded other, GN, AN and PKD). A portion of the data are shown below.
![[kidney1]](kidney1.bmp)
We have analysed the same data assuming a parametric Weibull distribution for the survivor function, and including an additive random effect b
i for each patient in the exponent of the hazard model as follows
t
ij ~ Weibull(r,
μij) i = 1,...,38; j = 1,2
log
μij =
α +
βageAGE
ij +
βsexSEX
i +
βdisease1DISEASE
i1 +
βdisease2DISEASE
i2 +
βdisease3DISEASE
i3 + b
i b
i ~ Normal(0,
τ)
where AGE
ij is a continuous covariate, SEX
i is a 2-level factor and DISEASE
ik (k = 1,2,3) are dummy variables representing the 4-level factor for underlying disease. Note that the the survival distribution is a truncated Weibull for censored observations as discussed in the mice example. The regression coefficients and the precision of the random effects
τ are given independent ``non-informative'' priors, namely
b
k ~ Normal(0, 0.0001)
τ ~ Gamma(0.0001, 0.0001)
The shape parameter of the survival distribution r is given a Gamma(1, 0.0001) prior which is slowly decreasing on the positive real line.
The graphical model and
BUGS language are given below.
Graphical model for kidney example: ![[kidney2]](kidney2.bmp)
BUGS language for kidney example
model
{
for (i in 1 : N) {
for (j in 1 : M) {
# Survival times bounded below by censoring times:
t[i,j] ~ dweib(r, mu[i,j])C(t.cen[i, j], );
log(mu[i,j ]) <- alpha + beta.age * age[i, j]
+ beta.sex *sex[i]
+ beta.dis[disease[i]] + b[i];
}
# Random effects:
b[i] ~ dnorm(0.0, tau)
}
# Priors:
alpha ~ dnorm(0.0, 0.0001);
beta.age ~ dnorm(0.0, 0.0001);
beta.sex ~ dnorm(0.0, 0.0001);
# beta.dis[1] <- 0; # corner-point constraint
for(k in 2 : 4) {
beta.dis[k] ~ dnorm(0.0, 0.0001);
}
tau ~ dgamma(1.0E-3, 1.0E-3);
r ~ dgamma(1.0, 1.0E-3);
sigma <- 1 / sqrt(tau); # s.d. of random effects
}
Data
list(N = 38, M = 2,
t = structure(
.Data = c( 8, 16,
23, NA,
22, 28,
447, 318,
30, 12,
24, 245,
7, 9,
511, 30,
53, 196,
15, 154,
7, 333,
141, NA,
96, 38,
NA, NA,
536, NA,
17, NA,
185, 177,
292, 114,
NA, NA,
15, NA,
152, 562,
402, NA,
13, 66,
39, NA,
12, 40,
NA, 201,
132, 156,
34, 30,
2, 25,
130, 26,
27, 58,
NA, 43,
152, 30,
190, NA,
119, 8,
NA, NA,
NA, 78,
63, NA), .Dim = c(38, 2)),
t.cen = structure(
.Data = c( 0, 0,
0, 13,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 8,
0, 0,
149, 70,
0, 25,
0, 4,
0, 0,
0, 0,
22, 159,
0, 108,
0, 0,
0, 24,
0, 0,
0, 46,
0, 0,
113, 0,
0, 0,
0, 0,
0, 0,
0, 0,
0, 0,
5, 0,
0, 0,
0, 5,
0, 0,
54, 16,
6, 0,
0, 8), .Dim = c(38, 2)),
age = structure(
.Data = c(28, 28,
48, 48,
32, 32,
31, 32,
10, 10,
16, 17,
51, 51,
55, 56,
69, 69,
51, 52,
44, 44,
34, 34,
35, 35,
42, 42,
17, 17,
60, 60,
60, 60,
43, 44,
53, 53,
44, 44,
46, 47,
30, 30,
62, 63,
42, 43,
43, 43,
57, 58,
10, 10,
52, 52,
53, 53,
54, 54,
56, 56,
50, 51,
57, 57,
44, 45,
22, 22,
42, 42,
52, 52,
60, 60), .Dim = c(38, 2)),
beta.dis = c(0, NA, NA, NA),
sex = c(0, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1,
0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0),
disease = c(1, 2, 1, 1, 1, 1, 2, 2, 3, 2, 3, 1, 3, 3, 1, 3, 1, 1,2, 1, 4, 1, 3, 3,
3, 3, 2, 3, 2, 2, 3, 3, 4, 2, 1, 1, 4, 4))
Inits for chain 1
list(beta.age = 0, beta.sex = 0, beta.dis=c(NA,0,0,0),
alpha = 0, r=1, tau=0.3)
Inits for chain 2
list(beta.age = -1, beta.sex = 1, beta.dis=c(NA,1,1,1),
alpha = 1, r=1.5, tau=1)
Results