BiRats

       Birats: a bivariate normal
      hierarchical model

We return to the Rats example, and illustrate the use of a multivariate Normal (MVN) population distribution for the regression coefficients of the growth curve for each rat. This is the model adopted by Gelfand etal (1990) for these data, and assumes a priori that the intercept and slope parameters for each rat are correlated. For example, positive correlation would imply that initially heavy rats (high intercept) tend to gain weight more rapidly (steeper slope) than lighter rats. The model is as follows

   Y _ij ~ Normal( m _ij, t _c)
    m _ij = b _1i + b _2i x _j
    b _i ~ MVN( m _b, W )

where Y _ij is the weight of the ith rat measured at age x _j, and b _i denotes the vector (b _1i, b _2i). We assume 'non-informative' independent univariate Normal priors for the separate components m _b₁ and m _b₂. A Wishart(R, r ) prior was specified for W , the population precision matrix of the regression coefficients. To represent vague prior knowledge, we chose the the degrees of freedom r for this distribution to be as small as possible (i.e. 2, the rank of W ). The scale matrix was specified as

This represents our prior guess at the order of magnitude of the covariance matrix W ^-1 for b _i (see Classic BUGS manual (version 0.5) section on Multivariate normal models), and is equivalent to the prior specification used by Gelfand et al. Finally, a non-informative Gamma(0.001, 0.001) prior was assumed for the measurement precision t _c.

[birats2]

   model
   {
       for( i in 1 : N ) {
          beta[i , 1 : 2] ~ dmnorm(mu.beta[], R[ , ])
          for( j in 1 : T ) {
             Y[i, j] ~ dnorm(mu[i , j], tauC)
             mu[i, j] <- beta[i, 1] + beta[i, 2] * x[j]
          }
       }

       mu.beta[1 : 2] ~ dmnorm(mean[], prec[ , ])
       R[1 : 2 , 1 : 2] ~ dwish(Omega[ , ], 2)
       tauC ~ dgamma(0.001, 0.001)
       sigma <- 1 / sqrt(tauC)
   }

Data    ( click to open )

Inits for chain 1        Inits for chain 2    ( click to open )

Results

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

[birats3]