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

# 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[1] <- 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)