Roger Schürch - 2013-04-17

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,

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!