library(lme4) Dyestuff <- data.frame(Batch = gl(6, 5, lab = LETTERS[1:6]), Yield = c( 1545, 1440, 1440, 1520, 1580, 1540, 1555, 1490, 1560, 1495, 1595, 1550, 1605, 1510, 1560, 1445, 1440, 1595, 1465, 1545, 1595, 1630, 1515, 1635, 1625, 1520, 1455, 1450, 1480, 1445)) str(Dyestuff) dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff, ylab = "Batch", jitter.y = TRUE, aspect = 0.3, type = c("p", "a")) (fm1 <- lmer(Yield ~ 1|Batch, Dyestuff)) set.seed(4376521) mcs1 <- mcmcsamp(fm1, 10000) xyplot(mcs1) HPDinterval(mcs1)