# “Mixed effect” ANOVA in R with JAGS/BUGS

I would like to be able to code the following linear mixed model into a bayesian framework. In the traditional/frequentist view, the model would be coded like this:

# Generation of simulated data
set.seed(123)
varY <- rnorm(100, 0, 1)
facX <- gl(n = 4, k = 25, labels = c("A", "B", "C", "D"))
block <- factor(rep(x = paste("B",1:25, sep = ""), times = 4))
df <- data.frame(varY, facX, block)

# Frequentist analysis
library(lme4)
model.freq <- lmer(varY ~ facX + (1|block), data = df)


where varY is a continuous response variable, facX is a categorical predictor with 4 levels and block is a random effect. This is similar to a one-way ANOVA with a fixed effect and a random block effect.

I can work out the code for the fixed effect one-way ANOVA:

# Data reshaping
matY <- matrix(data = varY, nrow = 25, ncol = 4, byrow = FALSE)

# Model specification
model.bayes <- function(){
# Likelihood
for(j in 1:Nlev){
for(i in 1:N){
varY[i,j] ~ dnorm(mu + theta[j], 1/(sig*sig))
}
}
theta <- 0
for(j in 2:Nlev){
theta[j] ~ dnorm(0, 0.001)
}
# Priors
mu ~ dnorm(0, 0.001)
sig ~ dunif(0, 1000)
}

# Bayesian analysis with JAGS
dat <- list(varY = matY, Nlev = ncol(matY), N = nrow(matY))
params <- c("mu", "theta", "sig")
inits <- function(){list(mu = rnorm(1), sig = rlnorm(1))}

library(R2jags)
out <- jags(data = dat, inits = inits, parameters.to.save = params, model.file = model.bayes, n.chains = 3, n.iter = 10000, n.burnin = 1000, n.thin = 1)


However, I don't know how to integrate the random block effect into this model. My question is: how would I integrate this block effect into the model to get a similar analysis as the frequentist one?

• I can't help with jags/bugs but note that your lmer model is not exactly "equivalent to a one-way ANOVA with a fixed effect and a random block effect"... – amoeba says Reinstate Monica Feb 13 '18 at 14:55
• Thanks @amoeba! You're right. I edited the comment to be more accurate. – Cédric_Frenette Feb 13 '18 at 17:31

In order to include a random effect (and potentially other fixed effects) you need to format your data in long (rather than wide) format, and use nested indexing with separate vectors as indicator variables for your fixed and random effects. The easiest way to generate this model is to use the runjags::template.jags function to create the model for you based on your lme4 model syntax:

# Generation of simulated data
set.seed(123)
varY <- rnorm(100, 0, 1)
facX <- gl(n = 4, k = 25, labels = c("A", "B", "C", "D"))
block <- factor(rep(x = paste("B",1:25, sep = ""), times = 4))
df <- data.frame(varY, facX, block)

library('runjags')
template.jags(varY ~ facX + (1|block), data = df)


You will then have a text file within your working directory that contains a fully functional JAGS model. The important parts are within the main loop:

regression_fitted[i] <- intercept + facX_effect[facX[i]] + block_randomeffect[block[i]]


This includes the fixed effect of facX and random effect of block, both indexed using the appropriate explanatory variable. The other part for the random effect is:

for(block_iterator in 1:25){
block_randomeffect[block_iterator] ~ dnorm(0, block_precision)
}


Which ensures that each random effect level is normally distributed with a mean of zero and precision to be estimated, which is equivalent to your frequentist model (although remember that precision is 1/variance!).

Once you have checked/edited the model as appropriate (e.g. adjusting priors) then you can run it using:

results <- run.jags("JAGSmodel.txt")


The help file for ?template.jags contains another example of comparing a frequentist linear model to an equivalent Bayesian model, which may be of interest to you.

• Wow! Thanks a lot @Matt Denwood! I was not aware of the template.jags function. This will certainly help me understand the JAGS syntax better. – Cédric_Frenette Feb 14 '18 at 19:17