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