Rcode: One-Way Random-Effects

Load libraries


Data: 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))

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

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

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)

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.


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


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

var(Dyestuff2$Yield)  ## How is sigma estimated