Pinheiro-Bates describes an experiment whereby the productivity of six randomly chosen workers are assessed three times on each of three machines, yielding the 54 observations tabulated below.
library(nlme)
data(Machines)
head(Machines)## Grouped Data: score ~ Machine | Worker
## Worker Machine score
## 1 1 A 52.0
## 2 1 A 52.8
## 3 1 A 53.1
## 4 2 A 51.8
## 5 2 A 52.8
## 6 2 A 53.1Machines %>% group_by(Worker,Machine) %>% tally()## # A tibble: 18 x 3
## # Groups: Worker [6]
## Worker Machine n
## <ord> <fct> <int>
## 1 6 A 3
## 2 6 B 3
## 3 6 C 3
## 4 2 A 3
## 5 2 B 3
## 6 2 C 3
## 7 4 A 3
## 8 4 B 3
## 9 4 C 3
## 10 1 A 3
## 11 1 B 3
## 12 1 C 3
## 13 3 A 3
## 14 3 B 3
## 15 3 C 3
## 16 5 A 3
## 17 5 B 3
## 18 5 C 3(Overall mean score = \(59.65\), mean on machine A = \(52.35\) , mean on machine B = \(60.32\), mean on machine C = \(66.27\))
Machines %>% group_by(Machine) %>% summarize( meanscore= mean(score))## # A tibble: 3 x 2
## Machine meanscore
## <fct> <dbl>
## 1 A 52.4
## 2 B 60.3
## 3 C 66.3The worker' factor is modelled with random effects($u_{i}$), whereas themachine’ factor is modelled with fixed effects (\(\beta_{j}\)).
Due to the repeated nature of the data, interaction effects between these factors are assumed to be extant, and shall be examined accordingly.
The interaction effect in this case (\(\tau_{ij}\)) describes whether the effect of changing from one machine to another is different for each worker.
The productivity score \(y_{ijk}\) is the \(k\)th observation taken on worker \(i\) on machine \(j\), and is formulated as follows;
\[\begin{equation} y_{ijk} = \beta_{j} + u_{i} + \tau_{ij} + \epsilon_{ijk} \end{equation}\] \[\begin{equation*} u_{i} \sim N(0, \sigma^{2}_{u}), \epsilon_{ijk} \sim N(0, \sigma^{2}), \tau_{i} \sim N(0, \sigma^{2}_{\tau}) \end{equation*}\]
The {nlme} package is incorporated into the R programming to perform linear mixed model calculations. For the `Machines’ data, Pinheiro-Bates use the following code, with the hierarchical structure specified in the last argument.
lme(score~Machine, data=Machines, random=~1|Worker/Machine)## Linear mixed-effects model fit by REML
## Data: Machines
## Log-restricted-likelihood: -107.8438
## Fixed: score ~ Machine
## (Intercept) MachineB MachineC
## 52.355556 7.966667 13.916667
##
## Random effects:
## Formula: ~1 | Worker
## (Intercept)
## StdDev: 4.78105
##
## Formula: ~1 | Machine %in% Worker
## (Intercept) Residual
## StdDev: 3.729532 0.9615771
##
## Number of Observations: 54
## Number of Groups:
## Worker Machine %in% Worker
## 6 18The output of the R computation is given below.
The crucial pieces of information given in the programme output are the estimates of the intercepts for each of the three machines.
Machine A, which is treated as a control case, is estimated to have an intercept of 52.35.
The intercept estimates for machines B and C are found to be \(60.32\) and \(66.27\) (by adding the values 7.96 and 13.91 to 52.35 respectively).
Estimatefor the variance components are also given; * \(\sigma^{2}_{u} = (4.78)^{2}\) * \(\sigma^{2}_{\tau} = (3.73)^{2}\) * \(\sigma^{2}_{\epsilon} = (0.96)^{2}\)