1. An industrial psychologist was interested in decreasing the time required to assemble an electronic component. Three types of assembly fixtures, including the one currently in use (control), were evaluated. For this study, six factories were selected randomly from a large population. The six factories were assigned randomly to the three assembly fixtures, with the constraint that each type of fixture was tested in two factories. Five workers in each factory were randomly selected from the entire workforce to participate in the study. After a three-week familiarization period, the experimenters measured the assembly rate (units assembled per hour) for each participant. In this experiment, fixture is fixed and factory is random, and factory is nested within fixture. The data are stored in the data frame electronics which can be loaded and inspected with the following commands:

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/electronics.rda") )
sapply(electronics,class)
##   fixture   factory    worker      rate 
##  "factor"  "factor"  "factor" "integer"
xtabs(~fixture + factory,data=electronics) # factory nested in fixture
##          factory
## fixture   f1 f2 f3 f4 f5 f6
##   control  5  5  0  0  0  0
##   design1  0  0  5  5  0  0
##   design2  0  0  0  0  5  5

1. Evaluate the statistical significance of the fixed effect of fixture and the random effect of factory.

# within anova:
electronics.aov.01 <- aov(rate~fixture+Error(factory),data=electronics)
summary(electronics.aov.01)
## 
## Error: factory
##           Df Sum Sq Mean Sq F value Pr(>F)  
## fixture    2   1460   730.0   9.955 0.0474 *
## Residuals  3    220    73.3                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 24    300    12.5
( F.factory <- 73.33/12.5 )
## [1] 5.8664
( p.factory <- 1-pf(F.factory,3,24) )
## [1] 0.003743147
# alternative anova:
electronics.aov.02 <- aov(rate~fixture+factory,data=electronics)
summary(electronics.aov.02)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## fixture      2   1460   730.0  58.400 6.02e-10 ***
## factory      3    220    73.3   5.867  0.00374 ** 
## Residuals   24    300    12.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
( F.fixture <- 730/73.33 )
## [1] 9.954998
( p.fixture <- 1-pf(F.fixture,2,3) )
## [1] 0.04738538

# within lmer:
library(lmerTest)
electronics.lme.01 <- lmer(rate~fixture+(1|factory),data=electronics)
anova(electronics.lme.01)

ranova(electronics.lme.01)

2. Compute association strength for fixture and factory.

library(effectsize)
# omega-squared for fixed effect:
omega_squared(electronics.lme.01) # correct

omega_squared(electronics.aov.01) # correct

omega_squared(electronics.aov.02) # NOT CORRECT!

# Last result is incorrect because the F stat for treatment is not
# correct in this model (we had to recalculate it)
F_to_omega2(f=9.95,df=2,df_error=3) # this gives correct value

# ICC for random effect
library(performance)
icc(electronics.lme.01,by_group=T)

library(VCA)
# Variance components are listed in column "VC"
# ICC is listed in column "%Total"
anovaMM(rate~fixture+(factory),Data=electronics)
## 
## 
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
## 
##  [Fixed Effects]
## 
##            int fixturecontrol fixturedesign1 fixturedesign2 
##             32            -17            -10              0 
## 
## 
##  [Variance Components]
## 
##   Name    DF       SS  MS        VC        %Total    SD       CV[%]    
## 1 total   8.019526               24.666667 100       4.966555 21.593717
## 2 factory 3        220 73.333333 12.166667 49.324324 3.488075 15.165543
## 3 error   24       300 12.5      12.5      50.675676 3.535534 15.371887
## 
## Mean: 23 (N = 30) 
## 
## Experimental Design: unbalanced  |  Method: ANOVA

2. An experiment analyzed the glycogen content of rat livers. There were three experimental treatment groups: a control group, a group treated with compound 217, and a group treated with compound 217 plus sugar. Six rats were selected randomly from a large population, and two were assigned randomly to each treatment. The experimenters took three samples of liver (i.e., three pieces of tissue) from each rat, and then two measures of glycogen were made on each sample for a total of 36 measures. The data are stored in the data frame liver which can be loaded with the following command:

load(url("http://pnb.mcmaster.ca/bennett/psy710/datasets/liver.rda") )
sapply(liver,class)
##  glycogen treatment       rat    sample   measure 
## "integer"  "factor"  "factor"  "factor"  "factor"
# summary(liver)
# xtabs(~treatment + rat,data=liver) # rat nested in treatment
# xtabs(~rat+sample,data=liver) # sample nested in rat

In this experiment, glycogen is the dependent variable, treatment is a fixed factor and rat and sample are random factors. Note that rat is nested within treatment and sample is nested within rat. The factor sample contains 18 unique levels, with each level representing one of the three samples taken from the six rats. The data frame also contains the factor measure, which you can think of as the ID for each glycogen measure (similar to the way subjects often are identified by an ID).

1. Evaluate the effect of treatment on glycogen while ignoring the random effects and nested structure of the data.

liver.aov.00 <- aov(glycogen~treatment,data=liver)
summary(liver.aov.00)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment    2   1558   778.8    14.5 3.03e-05 ***
## Residuals   33   1773    53.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. Evaluate the statistical significance of the fixed effect of treatment on glycogen using a model that includes the random/nested factors. Does the effect of treatment differ from the one estimated in the previous question? Why or why not?

Answer: The simplest way to perform this analysis is to use lmer.

# using lmer
require(lmerTest)
liver.lmer.01 <- lmer(glycogen~treatment+(1|rat)+(1|sample),data=liver)
anova(liver.lmer.01)

You also can use aov.

liver.aov.01 <- aov(glycogen ~ treatment + rat + sample, data=liver)
summary(liver.aov.01)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment    2 1557.6   778.8  36.793 4.38e-07 ***
## rat          3  797.7   265.9  12.562 0.000114 ***
## sample      12  594.0    49.5   2.339 0.050291 .  
## Residuals   18  381.0    21.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
( F.treatment <- 779/266 )
## [1] 2.928571
( p.treatment <- 1-pf(F.treatment,2,3))
## [1] 0.1971249

The effect of treatment was significant when we ignored the effects of the random and nested factors, but was not significant when we included the random and nested factors. The results differ because some of the variation across the levels of treatment actually are due to variations across levels of the nested factors rat and sample, and therefore comparing MS-glycogen to MS-residuals is a biased \(F\) test.

3. Evaluate the statistical significance of the random effects.

# using lmer
ranova(liver.lmer.01)

# using anova
summary(liver.aov.01)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## treatment    2 1557.6   778.8  36.793 4.38e-07 ***
## rat          3  797.7   265.9  12.562 0.000114 ***
## sample      12  594.0    49.5   2.339 0.050291 .  
## Residuals   18  381.0    21.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
( F.rat <- 265.9/49.5) # demoninator from sample
## [1] 5.371717
( p.rat <- 1-pf(F.rat,3,12))
## [1] 0.01410715
( F.sample <- 49.5/21.2)
## [1] 2.334906
( p.sample <- 1-pf(F.sample,12,18))
## [1] 0.0505995

Answer: Both lmer and aov find that the effect of rat is statistically significant, and the effect of sample is not significant. Note that the test of sample using lmer is almost significant whereas the one using aov is not even close. This difference is due to the fact that the null hypothesis is being evaluated with different statistics (i.e., chi-square vs. \(F\)) and the chi-square test generally is more conservative than the \(F\) test (Baayen, Davidson, and Bates 2008; Pinheiro and Bates 2000)

4. Estimate the variance components.

# variance components from VCA:
library(VCA)
liver.vca <- anovaMM(glycogen ~ treatment + (rat) + (sample),Data=liver)
print(liver.vca,digits=3)
## 
## 
## ANOVA-Type Estimation of Mixed Model:
## --------------------------------------
## 
##  [Fixed Effects]
## 
##                 int   treatmentaControl      treatmentbC217 treatmentcC217sugar 
##             135.167               5.333              15.833               0.000 
## 
## 
##  [Variance Components]
## 
##   Name   DF    SS      MS      VC     %Total SD    CV[%]
## 1 total  7.458                 71.398 100    8.45  5.941
## 2 rat    3     797.667 265.889 36.065 50.512 6.005 4.223
## 3 sample 12    594     49.5    14.167 19.842 3.764 2.646
## 4 error  18    381     21.167  21.167 29.646 4.601 3.235
## 
## Mean: 142.222 (N = 36) 
## 
## Experimental Design: unbalanced  |  Method: ANOVA

# variance components from anova table:
var.comp.error <- 21
var.comp.sample <- (50-21)/(2)
var.comp.rat <- (266-50)/(2*3) 
c(var.comp.rat,var.comp.sample,var.comp.error)
## [1] 36.0 14.5 21.0
sqrt(c(var.comp.rat,var.comp.sample,var.comp.error))
## [1] 6.000000 3.807887 4.582576

# variance components using lmer
vc <- VarCorr((liver.lmer.01))
print(vc,comp=c("Variance","Std.Dev."),digits=3)
##  Groups   Name        Variance Std.Dev.
##  sample   (Intercept) 14.2     3.76    
##  rat      (Intercept) 36.1     6.01    
##  Residual             21.2     4.60