Hierachical calibration with prediction of new values

[Roger Schürch]

The Newbie Question(s)

How do you model a hierarchical calibration curve in jags?

I am new to jags (and Bayesian statistics), and I am quite unsure whether I am doing the right thing, so any pointers to improve the model would be most welcome.

A more specific question concerns the prior for mux2 (see model below; which is derived from Hamada et al. (2003), J Quality Technology 35:194-205, on which I started building this model). If I expect the value to be positive, are there better ways to implement it?

Background

I would like to perform a calibration from known x and y of known individuals, and I would like to use this calibration to predict new x from observed y, without any knowledge of what individual is showing y.

I tried to follow Hamada et al. (2003), who do not have random intercepts, and I tried to elaborate from there:

R code

library("rjags")
library("lme4")

# create data set
sim.data <- data.frame(x = c(rep(100, 10), rep(200, 10), rep(300, 10), rep(400, 10)))
sim.data$indid <- c(rep(LETTERS[1:10], 4))
sim.data$ranintr <- c(rep(rnorm(10, 0, 10), 4))
sim.data$resid <- rnorm(40, 0, 30)

a <- 200
b <- 1

sim.data$y <- a + b*sim.data$x + sim.data$ranintr + sim.data$resid

# prepare simulated data for jags
N1 <- length(data.sim$y)
x <- sim.data$x
y <- sim.data$y
indid <- factor(sim.data$indid)
K <- length(unique(indid))

# new data for prediction
N2 <- 1
x2 <- c(NA)
y2 <- c(250)

jags <- jags.model('simModelCalibLinearHierarchical_predict.bug',
                   data = list('x' = x, 'y' = y,
                       'N1' = N1, 'K' = K, 'indid' = indid,
                       'N2' = N2, 'x2' = x2, 'y2' = y2),
                   n.chains = 2,
                   n.adapt = 100)

update(jags, 300000)

samples <- coda.samples(jags, c('alpha0', 'beta', 
                                'sigma2', 'sigma2.a.indid', 'ICC', 'x2'), 
                        100000, thin = 100)

fit.1 <- lmer(sim.data$y ~ sim.data$x + (sim.data$x | sim.data$indid))

summary(samples[,c('alpha0', 'beta', 'x2[1]', 'x2[2]')])
summary(fit.1)

JAGS code

model{

  # calibration
  for(i in 1:N1){       # the calibration observations
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- alpha0 + alpha[indid[i]] + beta * x[i]
  }
  # K should be the number of unique subjects
  for(i in 1:K){
    alpha[i] ~ dnorm(0, tau.a.indid)
  }

  # prediction
  for(i in 1:N2){       # the number of predictions
    y2[i] ~ dnorm(mu2[i], tau)
    mu2[i] <- alpha0 + alpha2[i] + beta * x2[i]
    x2[i] ~ dnorm(mux2, taux2)
    alpha2[i] ~ dnorm(mean.alpha, tau.a.indid)
  }

  # priors for calibration
  tau ~ dgamma(0.01, 0.01)
  alpha0 ~ dnorm(0, 0.0001)
  beta ~ dnorm(0, 0.0001) 
  tau.a.indid ~ dgamma(0.01, 0.01)

  # priors for predictions
  mux2 ~ dnorm(300, 0.0001)  # prior chosen so that (0, 600) contains 99.97% 
                             # of its distribution; this reflects our knowledge 
                             # that values must be positive 
  taux2 ~ dgamma(0.01, 0.01)

  # calculations
  mean.alpha <- mean(alpha)
  sigma2 <- 1 / tau
  sigma2.a.indid <- 1/tau.a.indid
  ICC <- sigma2.a.indid / (sigma2 + sigma2.a.indid)
}

In the end, I would like to do the equivalent to the lmer fit that can be found towards the end of R code block.

Again, any help is greatly appreciated!