Name of Contributed Example

Description of problem

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 , τ_a )

β_i ~ Normal( β_c, τ_b )

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 &beta_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 τ_a β_c, τ_b and τ_c are given independent "noninformative" priors. Interest particularly focuses on the intercept at zero time (birth), denoted α₀ = α_c - β_c x_bar.

BUGS Code

        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)
                                culmative.Y[i , j] <- culmative(Y[i , j], Y[i , j])
                                post.pv.Y[i , j] <- post.p.value(Y[i , j])
                                prior.pv.Y[i , j] <- prior.p.value(Y[i , j])
                                replicate.post.Y[i , j] <- replicate.post(Y[i , j])
                                pv.post.Y[i , j] <- step(Y[i , j] - replicate.post.Y[i , j])
                                replicate.prior.Y[i , j] <- replicate.prior(Y[i , j])
                                pv.prior.Y[i , j] <- step(Y[i , j] - replicate.prior.Y[i , j])
                        }
                        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.

Data

Data

Inits

Initial values

Results

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

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.