Load libraries

library(lme4) library(lattice)

**Data: Dyestuff**

data(Dyestuff) str(Dyestuff) ?Dyestuff

Dotplot of the data

dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff, ylab = "Batch", jitter.y = TRUE, pch = 21, xlab = "Yield of dyestuff (grams of standard color)", type = c("p", "a"))

ANOVA table and F-test for random-effects

out.anova = anova(lm(Yield ~ Batch, Dyestuff)) out.anova

Fit a random-effect model using ML and REML (default)

fm01 = lmer(Yield ~ 1 + (1 | Batch), Dyestuff) fm01 fm01ML = update(fm01, REML = FALSE) fm01ML

CIs for fixed parameters:

- Profile CIs (For profiling, we need to work with the ML estimates, so if you feed the command
`profile`

the REML fit, R will refit the model using ML); - For the intercept, we can construct an exact t-based CI (see notes);
- Wald CI. For example, the Wald interval for the intercept: mu-hat +/- 1.96xse(mu-hat), where the unknown parameters (sigma and sigma_A) in the expression of se are replaced by the estimated value. Wald-intervals tend to

underestimate the variability, since it ignores the uncertainty associated with the estimates of the variance components.

## Profle CI pr1 = profile(fm01ML) xyplot(pr1, aspect = 1.3) ## profile zeta plot confint(pr1) ## same as: confint(fm01ML) confint(pr1, level=0.99) pr1 round(pnorm(pr1[,1]), dig=3) ## The exact t-based CI for the intercept tmp = qt(0.975,5)*sqrt(out.anova[1,3]/30) c(mean(g.mean)-tmp, mean(g.mean)+tmp) ## Wald intervals confint(fm01, method="Wald") tmp = qnorm(0.975)*sqrt(out.anova[1,3]/30) c(mean(g.mean)-tmp, mean(g.mean)+tmp)

Assessing the random effects: The BLUP of a random effect is basically its “conditional mean” given the data, with parameters in the conditional mean replaced by their estimates. Since the conditional distribution is normal,

mean is also the mode.

ranef(fm01ML) g.mean = rowMeans(matrix(Dyestuff$Yield, ncol = 5, byrow=TRUE)) cbind(g.mean-mean(g.mean), ranef(fm01)$Batch, ranef(fm01ML)$Batch) dotplot(ranef(fm01ML, condVar=TRUE), strip=FALSE)

**Data: Dyestuff2**

data(Dyestuff2) str(Dyestuff2) ?Dyestuff2 ## Dotplot of the data dotplot(reorder(Batch, Yield) ~ Yield, Dyestuff2, ylab = "Batch", jitter.y = TRUE, pch = 21, xlab = "Yield of dyestuff (grams of standard color)", type = c("p", "a")) fm02 = lmer(Yield ~ 1 + (1 | Batch), Dyestuff2) fm02ML = lmer(Yield ~ 1 + (1 | Batch), Dyestuff2, REML = FALSE) summary(fm02) var(Dyestuff2$Yield) ## How is sigma estimated summary(fm02ML) var(Dyestuff2$Yield)*29/30