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,                 # regression coefficients for location models
 beta0.uncentred,
 R[M,M],                       # prior guess at magnitude of Sigma2 
 resid[N,M],resid2[N,M],RSS;   # residuals and residual sum of squares
model {
  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 + beta1 * (age[j] - mean(age)); # linear
  }
  beta0.uncentred <- beta0 - beta1 * mean(age);

  beta0 ~ dnorm(0.0, 0.001); 
  beta1 ~ 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] <- resid[i,j]^2;     # squared residuals
     } 
  }
  RSS <- sum(resid2[,]);                    # Residual Sum of Squares
}