Mixed-effects models for repeated-measures ANOVA

TODO
Install required packages
Simulate data for all designs
One-way repeated measures ANOVA (RB-\(p\) design)
Two-way repeated measures ANOVA (RBF-\(pq\) design)
Two-way split-plot-factorial ANOVA (SPF-\(p \cdot q\) design)
Three-way split-plot-factorial ANOVA (SPF-\(pq \cdot r\) design)
Three-way split-plot-factorial ANOVA (SPF-\(p \cdot qr\) design)
- Conventional analysis using aov()
- Mixed-effects analysis
Four-way split-plot-factorial ANOVA (SPF-\(pq \cdot rs\) design)
- Conventional analysis using aov()
- Mixed-effects analysis
Further resources
Detach (automatically) loaded packages (if possible)
Get the article source from GitHub

TODO

RBF-\(pq\): lme() with compound symmetry
SPF-\(p \cdot qr\): lme() with compound symmetry

Install required packages

AICcmodavg, lme4, multcomp, nlme, pbkrtest, lmerTest, performance

wants <- c("AICcmodavg", "lme4", "multcomp", "nlme", "pbkrtest",
           "lmerTest", "performance")
has   <- wants %in% rownames(installed.packages())
if(any(!has)) install.packages(wants[!has])

Simulate data for all designs

Two between-subjects factors, two within-subjects factors.

set.seed(123)
P     <- 2               # Xb1
Q     <- 2               # Xb2
R     <- 3               # Xw1
S     <- 3               # Xw2
Njklm <- 20              # obs per cell
Njk   <- Njklm*P*Q       # number of subjects
N     <- Njklm*P*Q*R*S   # number of observations
id    <- gl(Njk,         R*S, N, labels=c(paste("s", 1:Njk, sep="")))
Xb1   <- gl(P,   Njklm*Q*R*S, N, labels=c("CG", "T"))
Xb2   <- gl(Q,   Njklm  *R*S, N, labels=c("f", "m"))
Xw1   <- gl(R,             S, N, labels=c("A", "B", "C"))
Xw2   <- gl(S,   1,           N, labels=c("-", "o", "+"))

Theoretical main effects and interactions

mu      <- 100
eB1     <- c(-5, 5)
eB2     <- c(-5, 5)
eW1     <- c(-5, 0, 5)
eW2     <- c(-5, 0, 5)
eB1B2   <- c(-5, 5, 5, -5)
eB1W1   <- c(-5, 5, 2, -2, 3, -3)
eB1W2   <- c(-5, 5, 2, -2, 3, -3)
eB2W1   <- c(-5, 5, 2, -2, 3, -3)
eB2W2   <- c(-5, 5, 2, -2, 3, -3)
eW1W2   <- c(-5, 2, 3, 2, 3, -5, 2, -5, 3)
eB1B2W1 <- c(-5, 5, 5, -5, 2, -2, -2, 2, 3, -3, -3, 3)
eB1B2W2 <- c(-5, 5, 5, -5, 2, -2, -2, 2, 3, -3, -3, 3)
eB1W1W2 <- c(-5, 5, 2, -2, 3, -3, 3, -3, -5, 5, 2, -2, 2, -2, 3, -3, -5, 5)
eB2W1W2 <- c(-5, 5, 2, -2, 3, -3, 3, -3, -5, 5, 2, -2, 2, -2, 3, -3, -5, 5)
# no 3rd-order interaction B1xB2xW1xW2

Name values according to the corresponding cell in the experimental design

names(eB1)     <- levels(Xb1)
names(eB2)     <- levels(Xb2)
names(eW1)     <- levels(Xw1)
names(eW2)     <- levels(Xw2)
names(eB1B2)   <- levels(interaction(Xb1, Xb2))
names(eB1W1)   <- levels(interaction(Xb1, Xw1))
names(eB1W2)   <- levels(interaction(Xb1, Xw2))
names(eB2W1)   <- levels(interaction(Xb2, Xw1))
names(eB2W2)   <- levels(interaction(Xb2, Xw2))
names(eW1W2)   <- levels(interaction(Xw1, Xw2))
names(eB1B2W1) <- levels(interaction(Xb1, Xb2, Xw1))
names(eB1B2W2) <- levels(interaction(Xb1, Xb2, Xw2))
names(eB1W1W2) <- levels(interaction(Xb1, Xw1, Xw2))
names(eB2W1W2) <- levels(interaction(Xb2, Xw1, Xw2))

Simulate data given the effects defined above

muJKLM <- mu +
          eB1[Xb1] + eB2[Xb2] + eW1[Xw1] + eW2[Xw2] +
          eB1B2[interaction(Xb1, Xb2)] +
          eB1W1[interaction(Xb1, Xw1)] +
          eB1W2[interaction(Xb1, Xw2)] +
          eB2W1[interaction(Xb2, Xw1)] +
          eB2W2[interaction(Xb2, Xw2)] +
          eW1W2[interaction(Xw1, Xw2)] +
          eB1B2W1[interaction(Xb1, Xb2, Xw1)] +
          eB1B2W2[interaction(Xb1, Xb2, Xw2)] +
          eB1W1W2[interaction(Xb1, Xw1, Xw2)] +
          eB2W1W2[interaction(Xb2, Xw1, Xw2)]
muId  <- rep(rnorm(Njk, 0, 3), each=R*S)
mus   <- muJKLM + muId
sigma <- 50

Y  <- round(rnorm(N, mus, sigma), 1)
d2 <- data.frame(id, Xb1, Xb2, Xw1, Xw2, Y)

Data frame with just one within-subjects factor (average over levels of Xw2)

d1 <- aggregate(Y ~ id + Xw1 + Xb1 + Xb2, data=d2, FUN=mean)

One-way repeated measures ANOVA (RB-\(p\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xw1 + Error(id/Xw1), data=d1))


Error: id
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 79  75040   949.9               

Error: id:Xw1
           Df Sum Sq Mean Sq F value Pr(>F)  
Xw1         2   5756  2878.2   3.363 0.0371 *
Residuals 158 135211   855.8                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mixed-effects analysis

Using `lme()` from package `nlme`

No explicit assumption of compound symmetry, but random intercept model equivalent to compound symmetry iff all variance components are positive (id > id:Xw1 and IV > id:Xw1)

library(nlme)
anova(lme(Y ~ Xw1, random=~1 | id, method="ML", data=d1))

            numDF denDF   F-value p-value
(Intercept)     1   158 2554.7564  <.0001
Xw1             2   158    3.3633  0.0371

Assume compound symmetry

lmeFit <- lme(Y ~ Xw1, random=~1 | id,
              correlation=corCompSymm(form=~1|id),
              method="ML", data=d1)
anova(lmeFit)

            numDF denDF   F-value p-value
(Intercept)     1   158 2554.7627  <.0001
Xw1             2   158    3.3633  0.0371

Only 1 random effect, here all equivalent (not run):

anova(lme(Y ~ Xw1, random=list(id=pdIdent(~ 1)),       method="REML", data=d1))
anova(lme(Y ~ Xw1, random=list(id=pdDiag(~ 1)),        method="REML", data=d1))
anova(lme(Y ~ Xw1, random=list(id=pdCompSymm(~Xw1-1)), method="REML", data=d1))

Using `lmer()` from package `lme4`

By default, lmer() does not calculate p-values for fixed effects since the package author sees no convincing justification for selecting a particular value for the denominator (error) degrees of freedom.

library(lme4)
fitF <- lmer(Y ~ Xw1 + (1|id), data=d1)
summary(fitF)

Linear mixed model fit by REML ['lmerMod']
Formula: Y ~ Xw1 + (1 | id)
   Data: d1

REML criterion at convergence: 2294.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.09309 -0.69165 -0.03672  0.61179  2.77990 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  31.38    5.602  
 Residual             855.76   29.253  
Number of obs: 240, groups:  id, 80

Fixed effects:
            Estimate Std. Error t value
(Intercept)   93.629      3.330  28.116
Xw1B          10.364      4.625   2.241
Xw1C          10.414      4.625   2.252

Correlation of Fixed Effects:
     (Intr) Xw1B  
Xw1B -0.694       
Xw1C -0.694  0.500

anova(fitF)

Analysis of Variance Table
    npar Sum Sq Mean Sq F value
Xw1    2 5756.4  2878.2  3.3633

Package lmerTest enables Satterthwaite’s apprixmation of degrees of freedom for the summary() method.

library(lmerTest)
## need to re-fit model after lmerTest is loaded
fitF_p <- lmer(Y ~ Xw1 + (1|id), data=d1)
summary(fitF_p)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Y ~ Xw1 + (1 | id)
   Data: d1

REML criterion at convergence: 2294.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-3.09309 -0.69165 -0.03672  0.61179  2.77990 

Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept)  31.38    5.602  
 Residual             855.76   29.253  
Number of obs: 240, groups:  id, 80

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   93.629      3.330 236.408  28.116   <2e-16 ***
Xw1B          10.364      4.625 158.007   2.241   0.0264 *  
Xw1C          10.414      4.625 158.007   2.252   0.0257 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
     (Intr) Xw1B  
Xw1B -0.694       
Xw1C -0.694  0.500

Model comparison \(F\)-test with \(p\)-value with Kenward-Roger corrected degrees of freedom from package pbkrtest. Needs a full model and a restricted model with the effect of interest.

# restricted model
fitR <- lme4::lmer(Y ~ 1 + (1|id), data=d1)
library(pbkrtest)
KRmodcomp(fitF, fitR)

large : Y ~ Xw1 + (1 | id)
small : Y ~ 1 + (1 | id)
          stat      ndf      ddf F.scaling p.value  
Ftest   3.3633   2.0000 158.0000         1 0.03712 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

AIC comparison table using package `AICcmodavg`

Re-fit using ML instead of REML.

fitF_ML <- lme4::lmer(Y ~ Xw1 + (1|id), REML=FALSE, data=d1)
fitR_ML <- lme4::lmer(Y ~ 1   + (1|id), REML=FALSE, data=d1)

library(AICcmodavg)
AICc(fitF)

[1] 2304.445

aictab(cand.set=list(fitR_ML, fitF_ML),
       modnames=c("restricted", "full"),
       sort=FALSE, second.ord=FALSE)


Model selection based on AIC:

           K     AIC Delta_AIC AICWt       LL
restricted 3 2319.57      2.67  0.21 -1156.78
full       5 2316.90      0.00  0.79 -1153.45

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

library(performance)
r2_nakagawa(fitF)

# R2 for Mixed Models

  Conditional R2: 0.061
     Marginal R2: 0.026

icc(fitF)

# Intraclass Correlation Coefficient

     Adjusted ICC: 0.035
  Conditional ICC: 0.034

Multiple comparisons using `glht()` from package `multcomp`

library(multcomp)
contr <- glht(lmeFit, linfct=mcp(Xw1="Tukey"))
summary(contr)


     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = Y ~ Xw1, data = d1, random = ~1 | id, correlation = corCompSymm(form = ~1 | 
    id), method = "ML")

Linear Hypotheses:
           Estimate Std. Error z value Pr(>|z|)  
B - A == 0 10.36375    4.59638   2.255   0.0624 .
C - A == 0 10.41417    4.59638   2.266   0.0607 .
C - B == 0  0.05042    4.59638   0.011   0.9999  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)

confint(contr)


     Simultaneous Confidence Intervals

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = Y ~ Xw1, data = d1, random = ~1 | id, correlation = corCompSymm(form = ~1 | 
    id), method = "ML")

Quantile = 2.3448
95% family-wise confidence level
 

Linear Hypotheses:
           Estimate  lwr       upr      
B - A == 0  10.36375  -0.41370  21.14120
C - A == 0  10.41417  -0.36329  21.19162
C - B == 0   0.05042 -10.72704  10.82787

Two-way repeated measures ANOVA (RBF-\(pq\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))


Error: id
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 79 225120    2850               

Error: id:Xw1
           Df Sum Sq Mean Sq F value Pr(>F)  
Xw1         2  17269    8635   3.363 0.0371 *
Residuals 158 405633    2567                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw2
           Df Sum Sq Mean Sq F value Pr(>F)
Xw2         2   9859    4929   1.616  0.202
Residuals 158 481938    3050               

Error: id:Xw1:Xw2
           Df Sum Sq Mean Sq F value Pr(>F)
Xw1:Xw2     4   6118    1530   0.603  0.661
Residuals 316 802069    2538

Mixed-effects analysis

Using `lme()` from package `nlme`

anova(lme(Y ~ Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdIdent(~Xw1-1),
                                        pdIdent(~Xw2-1)))),
          method="ML", data=d2))

            numDF denDF   F-value p-value
(Intercept)     1   632 2440.2286  <.0001
Xw1             2   632    3.3889  0.0344
Xw2             2   632    1.6522  0.1924
Xw1:Xw2         4   632    0.6003  0.6626

Assume compound symmetry

anova(lme(Y ~ Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdCompSymm(~Xw1-1),
                                        pdCompSymm(~Xw2-1)))),
          method="ML", data=d2))

            numDF denDF   F-value p-value
(Intercept)     1   632 2554.7545  <.0001
Xw1             2   632    3.3633  0.0352
Xw2             2   632    1.6160  0.1995
Xw1:Xw2         4   632    0.6026  0.6609

Using `lmer()` from package `lme4`

fitF <- lme4::lmer(Y ~ Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2)
anova(fitF)

Analysis of Variance Table
        npar  Sum Sq Mean Sq F value
Xw1        2 17269.2  8634.6  3.3889
Xw2        2  8419.6  4209.8  1.6523
Xw1:Xw2    4  6118.0  1529.5  0.6003

anova(lmerTest::lmer(Y ~ Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))

Type III Analysis of Variance Table with Satterthwaite's method
         Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
Xw1     17269.2  8634.6     2   474  3.3889 0.03457 *
Xw2      8419.6  4209.8     2   237  1.6523 0.19382  
Xw1:Xw2  6118.0  1529.5     4   474  0.6003 0.66260  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

library(performance)
r2_nakagawa(fitF)

Random effect variances not available. Returned R2 does not account for random effects.

# R2 for Mixed Models

  Conditional R2: NA
     Marginal R2: 0.018

icc(fitF)

[1] NA

Two-way split-plot-factorial ANOVA (SPF-\(p \cdot q\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xb1*Xw1 + Error(id/Xw1), data=d1))


Error: id
          Df Sum Sq Mean Sq F value Pr(>F)  
Xb1        1   5335    5335    5.97 0.0168 *
Residuals 78  69705     894                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw1
           Df Sum Sq Mean Sq F value  Pr(>F)   
Xw1         2   5756    2878   3.541 0.03133 * 
Xb1:Xw1     2   8414    4207   5.176 0.00666 **
Residuals 156 126797     813                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mixed-effects

Using `lme()` from package `nlme`

Random intercept model equivalent to compound symmetry iff all variance components are positive

No explicit assumption of compound symmetry

anova(lme(Y ~ Xb1*Xw1, random=~1 | id, method="ML", data=d1))

            numDF denDF   F-value p-value
(Intercept)     1   156 2715.4696  <.0001
Xb1             1    78    5.9697  0.0168
Xw1             2   156    3.5411  0.0313
Xb1:Xw1         2   156    5.1762  0.0067

Assume compound symmetry

anova(lme(Y ~ Xb1*Xw1, random=~1 | id,
          correlation=corCompSymm(form=~1|id),
          method="ML", data=d1))

            numDF denDF   F-value p-value
(Intercept)     1   156 2715.4765  <.0001
Xb1             1    78    5.9697  0.0168
Xw1             2   156    3.5411  0.0313
Xb1:Xw1         2   156    5.1762  0.0067

Only 1 random effect, here all equivalent (not run):

anova(lme(Y ~ Xb1*Xw1, random=list(id=pdIdent(~ 1)),       method="REML", data=d1))
anova(lme(Y ~ Xb1*Xw1, random=list(id=pdDiag(~ 1)),        method="REML", data=d1))
anova(lme(Y ~ Xb1*Xw1, random=list(id=pdCompSymm(~Xw1-1)), method="REML", data=d1))

Using `lmer()` from package `lme4`

fitF <- lme4::lmer(Y ~ Xb1*Xw1 + (1|id), data=d1)
anova(fitF)

Analysis of Variance Table
        npar Sum Sq Mean Sq F value
Xb1        1 4852.2  4852.2  5.9697
Xw1        2 5756.4  2878.2  3.5411
Xb1:Xw1    2 8414.4  4207.2  5.1762

anova(lmerTest::lmer(Y ~ Xb1*Xw1 + (1|id), data=d1))

Type III Analysis of Variance Table with Satterthwaite's method
        Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)   
Xb1     4852.2  4852.2     1    78  5.9697 0.01681 * 
Xw1     5756.4  2878.2     2   156  3.5411 0.03133 * 
Xb1:Xw1 8414.4  4207.2     2   156  5.1762 0.00666 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

library(performance)
r2_nakagawa(fitF)

# R2 for Mixed Models

  Conditional R2: 0.118
     Marginal R2: 0.089

icc(fitF)

# Intraclass Correlation Coefficient

     Adjusted ICC: 0.032
  Conditional ICC: 0.029

Three-way split-plot-factorial ANOVA (SPF-\(pq \cdot r\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xb1*Xb2*Xw1 + Error(id/Xw1), data=d1))


Error: id
          Df Sum Sq Mean Sq F value  Pr(>F)   
Xb1        1   5335    5335   7.468 0.00781 **
Xb2        1   7246    7246  10.144 0.00210 **
Xb1:Xb2    1   8169    8169  11.436 0.00114 **
Residuals 76  54290     714                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw1
             Df Sum Sq Mean Sq F value   Pr(>F)    
Xw1           2   5756    2878   4.116 0.018166 *  
Xb1:Xw1       2   8414    4207   6.016 0.003058 ** 
Xb2:Xw1       2  11336    5668   8.105 0.000452 ***
Xb1:Xb2:Xw1   2   9167    4583   6.554 0.001861 ** 
Residuals   152 106294     699                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mixed-effects analysis

Using `lme()` from package `nlme`

anova(lme(Y ~ Xb1*Xb2*Xw1, random=~1 | id, method="ML", data=d1))

            numDF denDF  F-value p-value
(Intercept)     1   152 3397.096  <.0001
Xb1             1    76    7.468  0.0078
Xb2             1    76   10.144  0.0021
Xw1             2   152    4.116  0.0182
Xb1:Xb2         1    76   11.436  0.0011
Xb1:Xw1         2   152    6.016  0.0031
Xb2:Xw1         2   152    8.105  0.0005
Xb1:Xb2:Xw1     2   152    6.554  0.0019

Assume compound symmetry

anova(lme(Y ~ Xb1*Xb2*Xw1, random=~1 | id,
          correlation=corCompSymm(form=~1 | id),
          method="ML", data=d1))

            numDF denDF  F-value p-value
(Intercept)     1   152 3397.098  <.0001
Xb1             1    76    7.468  0.0078
Xb2             1    76   10.144  0.0021
Xw1             2   152    4.116  0.0182
Xb1:Xb2         1    76   11.436  0.0011
Xb1:Xw1         2   152    6.016  0.0031
Xb2:Xw1         2   152    8.105  0.0005
Xb1:Xb2:Xw1     2   152    6.554  0.0019

Only 1 random effect, here all equivalent (not run):

anova(lme(Y ~ Xb1*Xb2*Xw1, random=list(id=pdIdent(~ 1)),        method="ML", data=d1))
anova(lme(Y ~ Xb1*Xb2*Xw1, random=list(id=pdDiag(~ 1)),         method="ML", data=d1))
anova(lme(Y ~ Xb1*Xb2*Xw1, random=list(id=pdCompSymm(~Xw1-1)),  method="ML", data=d1))

Using `lmer()` from package `lme4`

fitF <- lme4::lmer(Y ~ Xb1*Xb2*Xw1 + (1|id), data=d1)
anova(fitF)

Analysis of Variance Table
            npar  Sum Sq Mean Sq F value
Xb1            1  5222.5  5222.5  7.4682
Xb2            1  7093.5  7093.5 10.1437
Xw1            2  5756.4  2878.2  4.1158
Xb1:Xb2        1  7997.0  7997.0 11.4357
Xb1:Xw1        2  8414.4  4207.2  6.0163
Xb2:Xw1        2 11335.9  5668.0  8.1052
Xb1:Xb2:Xw1    2  9166.6  4583.3  6.5541

anova(lmerTest::lmer(Y ~ Xb1*Xb2*Xw1 + (1|id), data=d1))

Type III Analysis of Variance Table with Satterthwaite's method
             Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
Xb1          5222.5  5222.5     1    76  7.4682 0.0078066 ** 
Xb2          7093.5  7093.5     1    76 10.1437 0.0021003 ** 
Xw1          5756.4  2878.2     2   152  4.1158 0.0181655 *  
Xb1:Xb2      7997.0  7997.0     1    76 11.4357 0.0011408 ** 
Xb1:Xw1      8414.4  4207.2     2   152  6.0163 0.0030580 ** 
Xb2:Xw1     11335.9  5668.0     2   152  8.1052 0.0004522 ***
Xb1:Xb2:Xw1  9166.6  4583.3     2   152  6.5541 0.0018608 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

library(performance)
r2_nakagawa(fitF)

# R2 for Mixed Models

  Conditional R2: 0.253
     Marginal R2: 0.248

icc(fitF)

# Intraclass Correlation Coefficient

     Adjusted ICC: 0.007
  Conditional ICC: 0.005

Three-way split-plot-factorial ANOVA (SPF-\(p \cdot qr\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xb1*Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))


Error: id
          Df Sum Sq Mean Sq F value Pr(>F)  
Xb1        1  16005   16005    5.97 0.0168 *
Residuals 78 209116    2681                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw1
           Df Sum Sq Mean Sq F value  Pr(>F)   
Xw1         2  17269    8635   3.541 0.03133 * 
Xb1:Xw1     2  25243   12622   5.176 0.00666 **
Residuals 156 380390    2438                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw2
           Df Sum Sq Mean Sq F value Pr(>F)
Xw2         2   9859    4929   1.604  0.204
Xb1:Xw2     2   2462    1231   0.400  0.671
Residuals 156 479476    3074               

Error: id:Xw1:Xw2
             Df Sum Sq Mean Sq F value Pr(>F)
Xw1:Xw2       4   6118    1530   0.601  0.662
Xb1:Xw1:Xw2   4   7609    1902   0.747  0.561
Residuals   312 794460    2546

Mixed-effects analysis

Using `lme()` from package `nlme`

anova(lme(Y ~ Xb1*Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdIdent(~Xw1-1),
                                        pdIdent(~Xw2-1)))),
          method="ML", data=d2))

            numDF denDF   F-value p-value
(Intercept)     1   624 2473.9526  <.0001
Xb1             1    78    5.4387  0.0223
Xw1             2   624    3.4396  0.0327
Xw2             2   624    1.6751  0.1881
Xb1:Xw1         2   624    5.0278  0.0068
Xb1:Xw2         2   624    0.4183  0.6584
Xw1:Xw2         4   624    0.6093  0.6561
Xb1:Xw1:Xw2     4   624    0.7577  0.5531

Assume compound symmetry

anova(lme(Y ~ Xb1*Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdCompSymm(~Xw1-1),
                                        pdCompSymm(~Xw2-1)))),
          method="ML", data=d2))

            numDF denDF   F-value p-value
(Intercept)     1   624 2715.3729  <.0001
Xb1             1    78    5.9695  0.0168
Xw1             2   624    3.4396  0.0327
Xw2             2   624    1.6038  0.2020
Xb1:Xw1         2   624    5.0278  0.0068
Xb1:Xw2         2   624    0.4005  0.6702
Xw1:Xw2         4   624    0.6093  0.6561
Xb1:Xw1:Xw2     4   624    0.7577  0.5531

Using `lmer()` from package `lme4`

anova(lme4::lmer(Y ~ Xb1*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))

Analysis of Variance Table
            npar  Sum Sq Mean Sq F value
Xb1            1 13654.3 13654.3  5.4390
Xw1            2 17269.2  8634.6  3.4395
Xw2            2  8410.8  4205.4  1.6752
Xb1:Xw1        2 25243.2 12621.6  5.0277
Xb1:Xw2        2  2100.3  1050.2  0.4183
Xw1:Xw2        4  6118.0  1529.5  0.6093
Xb1:Xw1:Xw2    4  7608.8  1902.2  0.7577

anova(lmerTest::lmer(Y ~ Xb1*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))

Type III Analysis of Variance Table with Satterthwaite's method
             Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
Xb1         13654.3 13654.3     1 233.98  5.4390 0.020541 * 
Xw1         17269.2  8634.6     2 467.96  3.4395 0.032894 * 
Xw2          8410.8  4205.4     2 233.98  1.6752 0.189514   
Xb1:Xw1     25243.2 12621.6     2 467.96  5.0277 0.006913 **
Xb1:Xw2      2100.3  1050.2     2 233.98  0.4183 0.658644   
Xw1:Xw2      6118.0  1529.5     4 467.96  0.6093 0.656152   
Xb1:Xw1:Xw2  7608.8  1902.2     4 467.96  0.7577 0.553229   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Four-way split-plot-factorial ANOVA (SPF-\(pq \cdot rs\) design)

Conventional analysis using `aov()`

summary(aov(Y ~ Xb1*Xb2*Xw1*Xw2 + Error(id/(Xw1*Xw2)), data=d2))


Error: id
          Df Sum Sq Mean Sq F value  Pr(>F)   
Xb1        1  16005   16005   7.468 0.00781 **
Xb2        1  21738   21738  10.144 0.00210 **
Xb1:Xb2    1  24507   24507  11.436 0.00114 **
Residuals 76 162871    2143                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw1
             Df Sum Sq Mean Sq F value   Pr(>F)    
Xw1           2  17269    8635   4.116 0.018166 *  
Xb1:Xw1       2  25243   12622   6.016 0.003058 ** 
Xb2:Xw1       2  34008   17004   8.105 0.000452 ***
Xb1:Xb2:Xw1   2  27500   13750   6.554 0.001861 ** 
Residuals   152 318882    2098                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw2
             Df Sum Sq Mean Sq F value  Pr(>F)   
Xw2           2   9859    4929   1.728 0.18110   
Xb1:Xw2       2   2462    1231   0.432 0.65031   
Xb2:Xw2       2  11822    5911   2.072 0.12945   
Xb1:Xb2:Xw2   2  34080   17040   5.974 0.00318 **
Residuals   152 433574    2852                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Error: id:Xw1:Xw2
                 Df Sum Sq Mean Sq F value Pr(>F)  
Xw1:Xw2           4   6118    1529   0.610 0.6558  
Xb1:Xw1:Xw2       4   7609    1902   0.759 0.5530  
Xb2:Xw1:Xw2       4  24545    6136   2.447 0.0465 *
Xb1:Xb2:Xw1:Xw2   4   7595    1899   0.757 0.5539  
Residuals       304 762320    2508                 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mixed-effects analysis

Using `lme()` from package `nlme`

No explicit assumption of compound symmetry

anova(lme(Y ~ Xb1*Xb2*Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdIdent(~Xw1-1),
                                        pdIdent(~Xw2-1)))),
          method="ML", data=d2))

                numDF denDF   F-value p-value
(Intercept)         1   608 2782.9276  <.0001
Xb1                 1    76    6.1180  0.0156
Xb2                 1    76    8.3098  0.0051
Xw1                 2   608    3.6417  0.0268
Xw2                 2   608    1.8843  0.1528
Xb1:Xb2             1    76    9.3682  0.0031
Xb1:Xw1             2   608    5.3232  0.0051
Xb2:Xw1             2   608    7.1714  0.0008
Xb1:Xw2             2   608    0.4705  0.6249
Xb2:Xw2             2   608    2.2595  0.1053
Xw1:Xw2             4   608    0.6451  0.6305
Xb1:Xb2:Xw1         2   608    5.7991  0.0032
Xb1:Xb2:Xw2         2   608    6.5138  0.0016
Xb1:Xw1:Xw2         4   608    0.8023  0.5240
Xb2:Xw1:Xw2         4   608    2.5880  0.0359
Xb1:Xb2:Xw1:Xw2     4   608    0.8008  0.5249

Assume compound symmetry

anova(lme(Y ~ Xb1*Xb2*Xw1*Xw2,
          random=list(id=pdBlocked(list(pdIdent(~1),
                                        pdCompSymm(~Xw1-1),
                                        pdCompSymm(~Xw2-1)))),
          method="ML", data=d2))

                numDF denDF   F-value p-value
(Intercept)         1   608 3113.1727  <.0001
Xb1                 1    76    6.8440  0.0107
Xb2                 1    76    9.2959  0.0032
Xw1                 2   608    3.6924  0.0255
Xw2                 2   608    1.7281  0.1785
Xb1:Xb2             1    76   10.4799  0.0018
Xb1:Xw1             2   608    5.3973  0.0047
Xb2:Xw1             2   608    7.2713  0.0008
Xb1:Xw2             2   608    0.4315  0.6497
Xb2:Xw2             2   608    2.0722  0.1268
Xw1:Xw2             4   608    0.6541  0.6242
Xb1:Xb2:Xw1         2   608    5.8798  0.0030
Xb1:Xb2:Xw2         2   608    5.9737  0.0027
Xb1:Xw1:Xw2         4   608    0.8134  0.5168
Xb2:Xw1:Xw2         4   608    2.6240  0.0339
Xb1:Xb2:Xw1:Xw2     4   608    0.8119  0.5178

Using `lmer()` from package `lme4`

anova(lme4::lmer(Y ~ Xb1*Xb2*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))

Analysis of Variance Table
                npar Sum Sq Mean Sq F value
Xb1                1  14506 14506.2  6.1180
Xb2                1  19703 19703.0  8.3098
Xw1                2  17269  8634.6  3.6417
Xw2                2   8936  4467.8  1.8843
Xb1:Xb2            1  22213 22212.6  9.3682
Xb1:Xw1            2  25243 12621.6  5.3232
Xb2:Xw1            2  34008 17003.9  7.1714
Xb1:Xw2            2   2231  1115.7  0.4705
Xb2:Xw2            2  10715  5357.5  2.2595
Xw1:Xw2            4   6118  1529.5  0.6451
Xb1:Xb2:Xw1        2  27500 13749.9  5.7991
Xb1:Xb2:Xw2        2  30889 15444.6  6.5138
Xb1:Xw1:Xw2        4   7609  1902.2  0.8023
Xb2:Xw1:Xw2        4  24545  6136.2  2.5880
Xb1:Xb2:Xw1:Xw2    4   7595  1898.7  0.8008

anova(lmerTest::lmer(Y ~ Xb1*Xb2*Xw1*Xw2 + (1|id) + (1|Xw1:id) + (1|Xw2:id), data=d2))

Type III Analysis of Variance Table with Satterthwaite's method
                Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)    
Xb1              14506 14506.2     1   228  6.1180 0.014111 *  
Xb2              19703 19703.1     1   228  8.3098 0.004320 ** 
Xw1              17269  8634.6     2   456  3.6416 0.026974 *  
Xw2               8936  4467.8     2   228  1.8843 0.154296    
Xb1:Xb2          22213 22212.6     1   228  9.3682 0.002473 ** 
Xb1:Xw1          25243 12621.6     2   456  5.3232 0.005185 ** 
Xb2:Xw1          34008 17003.9     2   456  7.1714 0.000858 ***
Xb1:Xw2           2231  1115.7     2   228  0.4705 0.625271    
Xb2:Xw2          10715  5357.5     2   228  2.2595 0.106731    
Xw1:Xw2           6118  1529.5     4   456  0.6451 0.630614    
Xb1:Xb2:Xw1      27500 13749.9     2   456  5.7991 0.003258 ** 
Xb1:Xb2:Xw2      30889 15444.6     2   228  6.5138 0.001774 ** 
Xb1:Xw1:Xw2       7609  1902.2     4   456  0.8023 0.524161    
Xb2:Xw1:Xw2      24545  6136.2     4   456  2.5880 0.036291 *  
Xb1:Xb2:Xw1:Xw2   7595  1898.7     4   456  0.8008 0.525103    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Further resources

For an alternative approach using generalized estimating equations (GEE), see package geepack.

Detach (automatically) loaded packages (if possible)

try(detach(package:multcomp))
try(detach(package:mvtnorm))
try(detach(package:TH.data))
try(detach(package:survival))
try(detach(package:pbkrtest))
try(detach(package:lmerTest))
try(detach(package:lme4))
try(detach(package:nlme))
try(detach(package:Matrix))
try(detach(package:AICcmodavg))
try(detach(package:MASS))
try(detach(package:performance))

Get the article source from GitHub

R markdown - markdown - R code - all posts

Mixed-effects models for repeated-measures ANOVA

TODO

Install required packages

Simulate data for all designs

One-way repeated measures ANOVA (RB-\(p\) design)

Conventional analysis using aov()

Mixed-effects analysis

Using lme() from package nlme

Using lmer() from package lme4

AIC comparison table using package AICcmodavg

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

Multiple comparisons using glht() from package multcomp

Two-way repeated measures ANOVA (RBF-\(pq\) design)

Conventional analysis using aov()

Mixed-effects analysis

Using lme() from package nlme

Using lmer() from package lme4

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

Two-way split-plot-factorial ANOVA (SPF-\(p \cdot q\) design)

Conventional analysis using aov()

Mixed-effects

Using lme() from package nlme

Using lmer() from package lme4

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

Three-way split-plot-factorial ANOVA (SPF-\(pq \cdot r\) design)

Conventional analysis using aov()

Mixed-effects analysis

Using lme() from package nlme

Using lmer() from package lme4

\(R^{2}\) (Nakagawa & Schielzeth) & ICC

Three-way split-plot-factorial ANOVA (SPF-\(p \cdot qr\) design)

Conventional analysis using aov()

Mixed-effects analysis

Using lme() from package nlme

Using lmer() from package lme4

Four-way split-plot-factorial ANOVA (SPF-\(pq \cdot rs\) design)

Conventional analysis using aov()

Mixed-effects analysis

Using lme() from package nlme

Using lmer() from package lme4

Further resources

Detach (automatically) loaded packages (if possible)

Get the article source from GitHub

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`

AIC comparison table using package `AICcmodavg`

Multiple comparisons using `glht()` from package `multcomp`

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`

Conventional analysis using `aov()`

Using `lme()` from package `nlme`

Using `lmer()` from package `lme4`