Rcode: Two-way Mixed Model


Machines Data (from package MEMSS)

Data on an experiment to compare three brands of machines
used in an industrial process are presented in Milliken
and Johnson (p. 285, 1992). Six workers were chosen randomly
among the employees of a factory to operate each machine three
times. The response is an overall productivity score taking
into account the number and quality of components produced.

library(lme4)
library(lattice)
library(MEMSS)

?Machines
Machines[1:5,]
str(Machines)

The data are completely crossed, which means that we have three
observations for each combination of a level of “Worker” and a level of “Machine.”

xtabs( ~ Worker + Machine, Machines)

dotplot(reorder(Worker, score) ~ score, Machines, groups = Machine,
        ylab = "Worker", xlab = "Productivity Score",
        type = c("p", "a"), auto.key = list(columns = 3, lines = TRUE))

Two-way Random Effects Model. For illustrative purposes, treat both worker and machine as random.

Model 1: two-way random-effects with interaction.

fr1 = lmer(score ~ 1 + (1|Worker) + (1|Machine) + 
     (1|Worker:Machine), data = Machines)
summary(fr1)

Check how to obtain the REML fit above by using the ANOVA output

machine.anova = anova(lm(score ~ Worker * Machine, Machines))

## REML estimates for the variance components
s2.e=machine.anova[4,3]
s2.WM=(machine.anova[3,3]-s2.e)/3
s2.W=(machine.anova[1,3]-machine.anova[3,3])/(3*3)
s2.M=(machine.anova[2,3]-machine.anova[3,3])/(3*6)

## estimate of the fixed effect and its sd
mean(Machines$score)
sqrt((machine.anova[1,3]+machine.anova[2,3] - machine.anova[3,3])/54)

## F-tests for random effects

## Test the random "Worker" effects
tmp.F = machine.anova[1,3]/machine.anova[3,3]
1 - pf(tmp.F, machine.anova[1,1], machine.anova[3,1])

## Test the random "Machine" effects
tmp.F = machine.anova[2,3]/machine.anova[3,3]
1 - pf(tmp.F, machine.anova[2,1], machine.anova[3,1])

## Test the random interaction effects
tmp.F = machine.anova[3,3]/machine.anova[4,3]
1 - pf(tmp.F, machine.anova[3,1], machine.anova[4,1])

## BLUP of random effects´┐╝
ranef(fr1)
myplot = dotplot(ranef(fr1, condVar=TRUE), strip=FALSE)
print(myplot[[1]], pos = c(0,0,1,0.75), more = TRUE)
print(myplot[[2]], pos = c(0,0.65,1,1))
print(myplot[[3]], pos = c(0,0.85,1,1))

Model 2: two-way random-effects without interaction

fr2 = lmer(score ~ 1 + (1|Worker)+
      (1|Machine), data = Machines)
summary(fr2)
anova(fr2, fr1)

Two-way Mixed Effects Model

  • Machine: fixed effects
  • Worker: random effects

Compare models with and without interaction

fm1 = lmer(score ~ Machine  + (1 | Worker), Machines, REML=FALSE)
fm2 = lmer(score ~ Machine + (1 | Worker) + 
            (1 | Worker:Machine), Machines, REML=FALSE)
fm3 = lmer(score ~ Machine + (Machine | Worker), Machines, REML=FALSE)
anova(fm1, fm2, fm3)

summary(fm1)
summary(fm2)
summary(fm3)

fm2 = lmer(score ~ Machine -1 + (1 | Worker) + 
     (1 | Worker:Machine), Machines)

## The fixed effects are the same as the sample means (for each Machine)
## The variance components and sd of the fixed effects can be computed
## through the ANOVA table
summary(fm2)$coef

mean(Machines$score[Machines$Machine=='A'])
mean(Machines$score[Machines$Machine=='B'])
mean(Machines$score[Machines$Machine=='C'])

summary(fm2)$varcor  ## standard dev instead of var
machine.anova
sqrt(machine.anova[4,3])
sqrt((machine.anova[3,3]-machine.anova[4,3])/3)
sqrt((machine.anova[1,3]-machine.anova[3,3])/9)