Rats: a normal hierarchical model

This example is taken from section 6 of Gelfand et al (1990), and concerns 30 young rats whose weights were measured weekly for five weeks. Part of the data is shown below, where Y_ij is the weight of the ith rat measured at age x_j.

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

The model is essentially a random effects linear growth curve

   Y_ij ~ Normal(α_i + β_i(x_j - x_bar), τ_c)

   α_i ~ Normal(α_c, τ_α)

   β_i ~ Normal(β_c, τ_β)

where x_bar = 22, and τ represents the precision (1/variance) of a normal distribution. We note the absence of a parameter representing correlation between α_i and β_i unlike in Gelfand et al 1990. However, see the Birats example in Volume 2 which does explicitly model the covariance between α_i and β_i. For now, we standardise the x_j's around their mean to reduce dependence between α_i and β_i 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).

α_c , τ_α , β_c , τ_β , τ_c are given independent "noninformative'' priors. Interest particularly focuses on the intercept at zero time (birth), denoted α₀ = α_c - β_c x_bar.
Graphical model for rats example:

BUGS language for rats example:
   model
   {
      for( i in 1 : N ) {
         for( j in 1 : T ) {
            Y[i , j] ~ dnorm(mu[i , j],tau.c)
            mu[i , j] <- alpha[i] + beta[i] * (x[j] - xbar)
         }
         alpha[i] ~ dnorm(alpha.c,alpha.tau)
         beta[i] ~ dnorm(beta.c,beta.tau)
      }
      tau.c ~ dgamma(0.001,0.001)
      sigma <- 1 / sqrt(tau.c)
      alpha.c ~ dnorm(0.0,1.0E-6)   
      alpha.tau ~ dgamma(0.001,0.001)
      beta.c ~ dnorm(0.0,1.0E-6)
      beta.tau ~ dgamma(0.001,0.001)
      alpha0 <- alpha.c - xbar * beta.c   
   }
Note the use of a very flat but conjugate prior for the population effects: a locally uniform prior could also have been used.

   

(Note: the response data (Y) for the rats example can also be found in the file ratsy in rectangular format. The covariate data (x) can be found in S-Plus format in file ratsx. To load data from each of these files, focus the window containing the open data file before clicking on "load data" from the "Specification" dialog.)
Results

These results may be compared with Figure 5 of Gelfand et al 1990 --- we note that the mean gradient of independent fitted straight lines is 6.19.

Gelfand et al 1990 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 obtained by simply replacing data values by NA (see below). The model specification is unchanged, since the distinction between observed and unobserved quantities is made in the data file and not the model specification.


Gelfand et al 1990 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 estimates were obtained:


We note that our estimate 6.58 of β_c is substantially greater than that shown in Figure 6 of Gelfand et al 1990. 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.