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.

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, m ij ) i = 1,...,38; j = 1,2
log m ij = a + b age AGE ij + b sex SEX i + b disease1 DISEASE i1 +
b disease2 DISEASE i2 + b disease3 DISEASE i3 + b i
b i ~ Normal(0, t )
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 t are given independent ``non-informative'' priors, namely
b k ~ Normal(0, 0.0001)
t ~ 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:
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]) I(t.cen[i, j], );
log(mu[i,j ]) <- alpha + beta.age * age[i, j]
+ beta.sex *sex[i]
+ beta.dis[disease[i]] + b[i];
culmative.t[i,j] <- culmative(t[i,j], t[i,j])
}
# 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 ( click to open )
Inits ( click to open )
Results
A 1000 update burn in followed by a further 10000 updates gave the parameter estimates