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)