Hepatitis: a normal hierarchical model with measurement error

This example is taken from Spiegelhalter et al (1996) (chapter in Markov Chain Monte Carlo in Practice) and concerns 106 children whose post-vaccination anti Hb titre was measured 2 or 3 times. Both measurements and times have been transformed to a log scale. One covariate y0 = log titre at baseline, is available.

The model is essentially a random effects linear growth curve

   Y_ij~ Normal(α_i + β_i (t_ij -t_bar), τ)

   α_i ~ Normal(α_c, τ_α)

   β_i ~ Normal(β_c, τ_β)

where τ 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.

α_c , τ_α , β_c , τ_β , τ are given independent ``noninformative'' priors.

Graphical model for hep example:

BUGS language for hep example:
model
{
for( i in 1 : N ) {
for( j in 1 : T ) {
Y[i , j] ~ dnorm(mu[i , j],tau)
mu[i , j] <- alpha[i] + beta[i] * (t[i,j] - 6.5) +
                           gamma * (y0[i] - mean(y0[]))
}
alpha[i] ~ dnorm(alpha0,tau.alpha)
beta[i] ~ dnorm(beta0,tau.beta)
}
tau ~ dgamma(0.001,0.001)
sigma <- 1 / sqrt(tau)
alpha0 ~ dnorm(0.0,1.0E-6)   
tau.alpha ~ dgamma(0.001,0.001)
beta0 ~ dnorm(0.0,1.0E-6)
tau.beta ~ dgamma(0.001,0.001)
gamma ~ dnorm(0.0,1.0E-6)
}

Note the use of a very flat but conjugate prior for the population effects: a locally uniform prior could also have been used.

   
Results





With measurement error


model
      {
         tau.alpha ~ dgamma(0.001,0.001)
         alpha0 ~ dnorm( 0.0,1.0E-6)
         beta0 ~ dnorm( 0.0,1.0E-6)
         tau.beta ~ dgamma(0.001,0.001)
         for( i in 1 : N ) {
            alpha[i] ~ dnorm(alpha0,tau.alpha)
            beta[i] ~ dnorm(beta0,tau.beta)
            y0[i] ~ dnorm(mu0[i],tau)
            mu0[i] ~ dnorm(theta,psi)
         }
         for( j in 1 : T ) {
            for( i in 1 : N ) {
               Y[i , j] ~ dnorm(mu[i , j],tau)
               mu[i , j] <- alpha[i] + beta[i] * (t[i , j] - 6.5) +
                  gamma * (mu0[i] - mean(y0[]))
            }
         }
         tau ~ dgamma(0.001,0.001)
         sigma <- 1 / sqrt(tau)
         gamma ~ dnorm( 0.0,1.0E-6)
         theta ~ dnorm( 0.0,1.0E-6)
         psi ~ dgamma(0.001,0.001)
      }    
Results