#### 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!