Birats: a bivariate Normal hierarchical model

We return to the rats example in Volume 1, 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:90 for this 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

where is the weight of the ith rat measured at age , and denotes the vector . We assume `non-informative' independent univariate Normal priors for the separate components and . A Wishart prior was specified for , the population precision matrix of the regression coefficients. To represent vague prior knowledge, we chose the the degrees of freedom for this distribution to be as small as possible (i.e. 2, the rank of ). The scale matrix R = . This represents our prior guess at the order of magnitude of the covariance matrix for (see BUGS manual 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 .

The appropriate graphical model is shown in Figure 8, and the BUGS code is given below. Note the use of the inverse function to compute Sigma2.beta, the population variance-covariance matrix for the regression coefficients. This matrix is then used to compute r, the correlation between the population mean intercept and slope parameters.

We note that in the original rats analysis we not only assumed were a priori independent, but also centred the covariates around their mean, which ensured that the likelihood for each pair factorised. By not centering the covariates and using a multivariate normal prior for , we have therefore introduced two additional forms of dependency.

We can investigate the influence of allowing such prior and likelihood dependence by fitting the range of models listed below:

• Likelihood for independent (centred covariates)
  • Prior independence of assumed: this is the original rats analysis, and implies that and retain independence even conditional on the data.
  • Prior dependence of allowed: this leads to dependence in the posterior of and to be introduced, essentially due to the estimates of 's and 's being correlated (in fact, fitting completely independent growth curves to the 30 rats leads to an empirical correlation of .50 between the estimated slopes and estimated intercepts, even though each specific pair of estimates are independent).

• Likelihood for dependent (uncentred covariates)
  • Prior independence of assumed: in this formulation dependence between is only introduced through the likelihood: when learning about , say, only current values of 's are taken into account.
  • Prior dependence of allowed: this is the full multivariate analysis described above.

The BUGS code for all these combinations is given below: relevant models may be fitted by changing the commented lines appropriately. We note that the definition of , and depends on whether the likelihood is independent (covariates centred) or not.

Birats: model specification in BUGS

model birats;

const
   N = 30,  # number of rats
   T = 5;   # number of time points

var
   x[T],mu[N,T],Y[N,T],beta[N,2],mu.beta[2],Omega.beta[2,2],
   Sigma2.beta[2,2],sigma.beta[2],tau.c,sigma,R[2,2],r,alpha0,
   tau.beta[2];

data Y in "biratsy.dat", x in "biratsx.dat";
inits in "birats.in"; 

{
    for (i in 1:N) {
      for (j in 1:T) {
        Y[i,j] ~ dnorm(mu[i,j],tau.c);  #
        mu[i,j] <- beta[i,1] + beta[i,2]*x[j];             # uncentred
#       mu[i,j] <- beta[i,1] + beta[i,2]*(x[j]-mean(x[])); # centred
     }
     beta[i,] ~ dmnorm(mu.beta[],Omega.beta[,]);  # bivariate Normal
#    beta[i,1] ~ dnorm(mu.beta[1],tau.beta[1]);  # independent Normals
#    beta[i,2] ~ dnorm(mu.beta[2],tau.beta[2]);  # independent Normals
     }
 
# intercept at zero for centred model
#  alpha0 <- mu.beta[1] - mu.beta[2]* mean(x[]); 

# intercept at mean(x[]) for uncentred model
   alpha0 <- mu.beta[1] + mu.beta[2]* mean(x[]); 
# prior for sampling precision
  tau.c  ~ dgamma(1.0E-3,1.0E-3);  sigma <- 1.0/sqrt(tau.c);



# parameters considered MVN

 Omega.beta[,]  ~ dwish(R[,],2);    # Wishart prior on precision matrix
  R[1,1] <- 200.0; R[1,2] <- 0.0;   # R = prior guess at order of covariance
  R[2,1] <- 0.0; R[2,2] <- 0.2;     # matrix for beta[i,]  
  Sigma2.beta[,] <-inverse(Omega.beta[,]);
  sigma.beta[1]<-sqrt(Sigma2.beta[1,1]);
  sigma.beta[2]<-sqrt(Sigma2.beta[2,2]);
  r <- Sigma2.beta[1,2] / (sqrt(Sigma2.beta[1,1])
                          *sqrt(Sigma2.beta[2,2]));  # correlation

# parameters considered independent
# for (k in 1:2){
#  tau.beta[k]  ~ dgamma(1.0E-3,1.0E-3);
#  sigma.beta[k]<-1/sqrt(tau.beta[k]);
# }
 
  mu.beta[1] ~ dnorm(0,.00001);   # `flat' univariate Normal prior on mean
  mu.beta[2] ~ dnorm(0,.00001);   # `flat' univariate Normal prior on mean 
}

 
Figure 8:   Graphical model for the birats example

Analysis

Each run used a 500 iteration burn-in and a further 5000 iterations, which took around 90 seconds and yielded the following results.

We can make a number of observations from these results.

1. The estimated correlation r between slopes and intercepts is high for the centred data and almost zero for the uncentred, which suggests a `fan' pattern in which the speed of growth is unrelated to initial weight at birth, but as the lines separate the intercept at is higher for those with faster growth.
2. Independence assumptions make no difference to the main location parameter estimates, only their precision (and convergence properties).
3. For centred data, the population intercept at 0 is substantially less precise with the prior independence model than allowing dependence, since it is a function of and and the dependence introduced between them from the highly correlated 's has been ignored.
4. For uncentred data, the population intercept at is only slightly less precise with the prior independence model than allowing dependence, since it is a function of and and the very small dependence introduced between them from the correlated 's has been ignored.

In general we would advise using the multivariate normal model for multiple random effects, particularly when interest lies in functions of the population coefficients, such as predictions.

It is also technically possible to assume a multivariate sampling model for the 's, in which correlated residuals around the growth curve are permitted. However, we have found these models have very poor convergence properties: initial poor estimates lead to the variance-covariance matrix of the observation vector having very large entries, and once this has occurred the sampler tends to `stick' with these values since the fitted curves have little incentive to approach the data. If such models are attempted, extremely good starting values would be necessary, perhaps derived from the type of analysis above.