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