Introduction to Mixed Modelling: Beyond Regression and Analysis of Variance, 2nd Edition

7.9 Use of R to fit the random-effects model with several factors

The following statements import the data into R, convert ‘location’, ‘sowing_occasion’, ‘sowing_date’ and ‘variety’ to factors and fit Model 7.4:

rm(list = ls())
canola.oil <- read.table(
   "IMM edn 2\\Ch 7\\canola oil gxe.txt",
   header=TRUE)
attach(canola.oil)
flocation <- factor(location)
fsowing_occasion <- factor(sowing_occasion)
fsowing_date <- factor(sowing_date)
fvariety <- factor(variety)
library(lme4)
canola.model1lmer <- lmer(oil ~ fsowing_occasion +
   (1|flocation) + (1|fvariety) +
   (1|fsowing_occasion:flocation) + (1|fsowing_occasion:fvariety) +
   (1|flocation:fvariety),
   data = canola.oil)
summary(canola.model1lmer)
anova(canola.model1lmer)

Note that the three-way interaction term ‘sowing_occasion:location:variety’ is not included explicitly in the model, as it was in the first GenStat model, but is recognized from the outset as the residual term. This is necessary in order to obtain a correct analysis from the function lmer().

The output of the summary() and anova() functions is as follows:

Linear mixed model fit by REML Formula: oil ~ fsowing_occasion + (1 | flocation) + (1 | fvariety) + (1 | fsowing_occasion:flocation) + (1 | fsowing_occasion:fvariety) + (1 | flocation:fvariety) Data: canola.oil AIC BIC logLik deviance REMLdev 378.1 412.1 -177.1 362.6 354.1 Random effects: Groups Name Variance Std.Dev. fsowing_occasion:fvariety (Intercept) 0.11829 0.34393 flocation:fvariety (Intercept) 0.11481 0.33883 ...