Rats: Normal hierarchical models with missing data

This example is taken from section 6 of gelfand:etal:90, and concerns 30 young rats whose weights were measured weekly for five weeks. Part of the data is shown below, where is the weight of the ith rat measured at age .

A plot of the 30 growth curves suggests some evidence of downward curvature.

The model is essentially a random effects linear growth curve

where , and represents the precision (1/variance) of a normal distribution. We note the absence of a parameter representing correlation between and unlike in gelfand:etal:90. However, see the birats example in Volume 2 which does explicitly model the covariance between and . For now, we standardise the 's around their mean to reduce dependence between and in their likelihood: in fact for the full balanced data, complete independence is achieved. (Note that, in general, prior independence does not force the posterior distributions to be independent).

are given independent ``noninformative'' priors. Interest particularly focusses on the intercept at zero time (birth), denoted . The appropriate graphical model is shown in Figure 1.

 
Figure 1:   Graphical model for rats example

Rats: model specification in BUGS

model rats;
const
   N = 30,  # number of rats
   T = 5;   # number of time points
var
   tau.c, alpha0, alpha.c, beta.c, x[T],
   mu[N,T], Y[N,T], alpha[N], beta[N],
   tau.alpha, tau.beta, sigma, x.bar;

data Y in "ratsy.dat", x in "ratsx.dat";
inits in "rats.in";

{
    for (i in 1:N) {
       for (j in 1:T) {
          mu[i,j] <- alpha[i] + beta[i]*(x[j] - x.bar);
          Y[i,j]   ~ dnorm(mu[i,j],tau.c)
       }
       alpha[i] ~ dnorm(alpha.c,tau.alpha);
       beta[i]  ~ dnorm(beta.c,tau.beta);
    }
    alpha.c   ~ dnorm(0,1.0E-4);
    beta.c    ~ dnorm(0,1.0E-4);
    tau.c     ~ dgamma(1.0E-3,1.0E-3);
    tau.alpha ~ dgamma(1.0E-3,1.0E-3);
    tau.beta  ~ dgamma(1.0E-3,1.0E-3);
    sigma    <- 1.0/sqrt(tau.c);
    x.bar    <- mean(x[]);
    alpha0   <- alpha.c - beta.c*x.bar;
}

Note the use of a very flat but conjugate prior for the population effects: a locally uniform prior could also have been used.

If the data are input in rectangular format, 2 files are required. The response data y are in file ratsy.dat:

 151 199 246 283 320
 145 199 249 293 354
 147 214 263 312 328
 155 200 237 272 297
 .....
 .....
 157 205 248 289 316
 137 180 219 258 291
 153 200 244 286 324

and the measurements times x are in ratsx.dat:

8.0
15.0
22.0
29.0
36.0

Alternatively, the data may be input in S format, as in file ratsS.dat. In this case, both y and x may be included in the same file:

list(Y =  c(151.0,199.0,246.0,283.0,320.0,
            145.0,199.0,249.0,293.0,354.0,
            .....
            153.0,200.0,244.0,286.0,324.0),
      x = c(8.0,15.0,22.0,29.0,36.0))

and the data statement on line 10 of the rats.bug file must be changed to

data in "ratsS.dat";

Analysis

A naive run, using no diagnostics for convergence, gave the following results for the population intercept at time 0 and the population gradient .

Bugs>update(500)   500    updates took  00:00:02 
Bugs>monitor(alpha0)
Bugs>monitor(beta.c)
Bugs>update(1000)   1000   updates took  00:00:04 
Bugs>stats(alpha0)
                  mean        sd       2.5% :  97.5%  CI    median      sample
                1.063E+2   3.590E+0   9.968E+1   1.132E+2   1.061E+2     1000  
Bugs>stats(beta.c)
                  mean        sd       2.5% :  97.5%  CI    median      sample
                6.183E+0   1.095E-1   5.968E+0   6.393E+0   6.179E+0     1000

These results may be compared with Figure 5 of gelfand:etal:90 -- we note that the mean gradient of independent fitted straight lines is 6.19.

gelfand:etal:90 also consider the problem of missing data, and delete the last observation of cases 6-10, the last two from 11-20, the last 3 from 21-25 and the last 4 from 26-30. The appropriate data file is called ratsmiss.dat, and is obtained by simply replacing data values by NA (see below). The rats.bug file only has to change the data declaration to data Y in "ratsmiss.dat"; we note that this is the only change necessary, since the distinction between observed and unobserved quantities is made in the data file and not the model specification.

Data file "ratsmiss.dat"

 151 199 246 283 320
 145 199 249 293 354
 .....
 153  NA  NA  NA  NA

gelfand:etal:90 focus on the parameter estimates and the predictions for the final 4 observations on rat 26. These predictions are obtained automatically in BUGS by monitoring the relevant Y[] nodes. The following is a sample run.

Bugs>update(500)   500    updates took  00:00:02 
Bugs>monitor(beta.c)
Bugs>monitor(Y[26,])
Bugs>update(1000)   1000   updates took  00:00:04 
Bugs>stats(beta.c)
                  mean        sd       2.5% :  97.5%  CI    median      sample
                6.537E+0   1.411E-1   6.260E+0   6.811E+0   6.533E+0     1000  
Bugs>stats(Y[26,])
                  mean        sd       2.5% :  97.5%  CI    median      sample
[26,2]          2.044E+2   8.937E+0   1.865E+2   2.212E+2   2.046E+2     1000  
[26,3]          2.497E+2   1.076E+1   2.294E+2   2.706E+2   2.493E+2     1000  
[26,4]          2.952E+2   1.280E+1   2.700E+2   3.216E+2   2.949E+2     1000  
[26,5]          3.413E+2   1.603E+1   3.115E+2   3.742E+2   3.401E+2     1000

We note that our estimate 6.54 of is substantially greater than that shown in Figure 6 of gelfand:etal:90. However, plotting the growth curves indicates some curvature with steeper gradients at the beginning: the mean of the estimated gradients of the reduced data is 6.66, compared to 6.19 for the full data. Hence we are inclined to believe our analysis. The observed weights for rat 26 were 207, 257, 303 and 345, compared to our predictions of 204, 250, 295 and 341.