[kidney0]     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

ij ~ Weibull(r, m ij ) i = 1,...,38; j = 1,2
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
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   

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:    [kidney2]

BUGS language for kidney example

      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];
            cumulative.t[i,j] <- cumulative(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 for chain 1        Inits for chain 2    ( click to open )


A 1000 update burn in followed by a further 10000 updates gave the parameter estimates