welcome: please sign in
location: RatsExample

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 Yij is the weight of the ith rat measured at age xj .

Weight Yij of rat i on day xj

xj = 8

15

22

29

36

Rat 1

151

199

246

283

320

Rat 2

145

199

249

293

354

...

Rat 1

153

200

244

286

324

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

The model is essentially a random effects linear growth curve

Yij ~ Normal( αi + βi (xj - xbar ), τc )

αi ~ Normal( αc , τa )

βi ~ Normal( βc, τb )

where xbar = 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 xj's around their mean to reduce dependence between αi and &betai 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 α0 = αc - βc xbar.

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:

mean

sd

MC_error

val2.5pc

median

val97.5pc

start

sample

alpha0

106.6

3.655

0.04079

99.44

106.5

113.8

1001

10000

beta.c

6.185

0.1061

0.00132

5.975

6.185

6.394

1001

10000

sigma

6.086

0.4606

0.007398

5.255

6.061

7.049

1001

10000

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.

RatsExample (last edited 2009-10-26 15:59:15 by BobOH)