elston:grizzle:62 present repeated measurements of ramus (jaw) bone height on a cohort of 20 boys over an 18 month period:
Interest focuses on describing the average growth curve of the ramus bone. The 4 measurements 
= { 
} for each child i are assumed to be correlated and follow a multivariate normal (MVN) distribution with unknown population mean vector 
and precision matrix 
. That is 
The following location models for the population mean were fitted in turn:
where 
= age at jth measurement. Non-informative independent normal priors were specified for the regression coefficients 
, 
and 
. The population precision matrix was assumed to follow a Wishart 
distribution. To represent vague prior knowledge, we chose the the degrees of freedom 
for this distribution to be as small as possible (i.e. 4, the rank of ). The scale matrix R was specified as 
, which represents an assessment of the order of magnitude of the covariance matrix 
for (see subsection on the use of the Wishart distribution in the ``Multivariate normal nodes'' section of the BUGS manual 0.50). Note that except for cases with very few individuals, the choice of R has little effect on the posterior estimate of [Lindley1970].
Comparision of the fit of the 3 location models may be assessed by calculating the deviance. This is given by -twice the sum of the log-likelihood contributions for each boy i:
where M=4, the number of measurements per boy. The change in (minimum) deviance between the constant and linear models or the linear and quadratic models may be compared to a 
distribution on 1 degree of freedom. In addition, we may compute the residual sum of squares RSS 
for each model.
The graph of the quadratic model is shown in Figure 7, and the BUGS code is given below.
model jaw;
const
M = 4, # number of time points
N = 20, # number of boys
PI = 3.141593;
var
Y[N,M], age[M], # jaw bone length in mm and age
mu[M],Omega[M,M],Sigma2[M,M], # mean, precision & covariance matrix for Y
beta0, beta1, beta2, # regression coefficients for location models
R[M,M], # prior guess at magnitude of Sigma2
resid[N,M],resid2[N,M],RSS, # residuals and residual sum of squares
L1[N,M],llike[N],deviance; # log-likelihood terms and deviance
data Y in "jawy.dat", age, R in "jawcov.dat";
inits in "jaw.in";
{
for (i in 1:N) {
Y[i,] ~ dmnorm(mu[], Omega[,]); # The 4 measurements for each
} # boy are multivariate normal
for(j in 1:M) { # location model for mean bone length at each age
# mu[j] <- beta0; # constant
mu[j] <- beta0 + beta1*age[j]; # linear
# mu[j] <- beta0 + beta1*age[j] + beta2*pow(age[j],2); # quadratic
}
beta0 ~ dnorm(0.0, 0.001);
beta1 ~ dnorm(0.0, 0.001);
beta2 ~ dnorm(0.0, 0.001);
Omega[,] ~ dwish(R[,], 4); # between-child variance in length at each age
Sigma2[,] <- inverse(Omega[,]);
for (i in 1:N) {
for (j in 1:M) {
resid[i,j] <- Y[i,j] - mu[j]; # residuals
resid2[i,j] <-pow(resid[i,j], 2); # squared residuals
L1[i,j] <- inprod(Omega[j,], resid[i,]);
}
llike[i] <- -0.5*(M*log(2*PI) - logdet(Omega[,])
+ inprod(resid[i,], L1[i,]));
}
RSS <- sum(resid2[,]); # Residual Sum of Squares
deviance <- -2 * sum(llike[]); # Deviance
}
Note that at present, BUGS is unable to perform matrix multiplication. Hence to compute
we use the inprod function to first compute 
, and then to compute 
. Also note the use of the inverse function to compute Sigma2.
Figure 7:
Graphical model for jaw example.
Analysis
After a 500 iteration burn-in, 1000 iterations took between 12 and 45 seconds for the 3 models and produced the following output {
MLE = Maximum likelihood estimates obtainied by goldstein:79 (p. 95), and prosser:rasbash:goldstein (p. 106) using the ML3 software.
Examination of the RSS clearly indicates that the linear model is a superior fit to the constant model, but that a quadratic term is unecessary. This is confirmed by the change in (minimum) deviance between successive models: the linear versus constant model yields a log likelihood ratio statistic of 32.4 on 1 d.f. (p<0.000001). The quadratic versus linear model yield virtually identical deviances, giving a non-significant log likelihood ratio statistic.