2
$\begingroup$

I have traditionally been running my analyses using aov and would like to switch to lmer. My questions is whether or not this is the correct procedure (see below).

Assume a data set which consist of a 3(A1 vs. A2 vs. A3) x 2(B vs. C) x 2(Y vs. Z) mixed design, with the first factor being between-subjects and the two others within-subjects. The two within-subjects factors are crossed with each other. As far as I know, the only random factor is participant's id. I am interested in the fixed effects.

Here is a reproducible example:

library(data.table)
library(tidyr)
library(lme4)
library(lmerTest)

# -- Data in wide format --
set.seed(33)

X <- data.table( id = paste0("P", 1:12),                  # participant unique ID
                 A  = rep(c("A1", "A2", "A3"), each = 4),
                 BY = rnorm(12, 1),
                 BZ = rnorm(12, 3),
                 CY = rnorm(12, 2),
                 CZ = rnorm(12, 0) )

# -- Orthogonal contrast codes --
X[A == "A1", A1vsA2 := +1] # A1 vs. A2
X[A == "A2", A1vsA2 := -1]
X[A == "A3", A1vsA2 :=  0]

X[A == "A1", A3vsA1A2 := +1] # A3 vs. A1&A2
X[A == "A2", A3vsA1A2 := +1]
X[A == "A3", A3vsA1A2 := -2]

# -- Convert to long format --
XL <- data.table( gather(X, KEY, Y, BY, BZ, CY, CZ) )

XL[KEY %in% c("BY", "BZ"), BvsC := +1]
XL[KEY %in% c("CY", "CZ"), BvsC := -1]

XL[KEY %in% c("BY", "CY"), YvsZ := +1]
XL[KEY %in% c("BZ", "CZ"), YvsZ := -1]

# -- Test model --
aovMdl  <-  aov(Y ~ A1vsA2*A3vsA1A2*BvsC*YvsZ + Error(id/(BvsC*YvsZ)), data = XL, REML = FALSE)
lmerMdl1 <- lmer(Y ~ A1vsA2*A3vsA1A2*BvsC*YvsZ + (1|id),                data = XL, REML = FALSE)

>summary(aovMdl)

Error: id
          Df Sum Sq Mean Sq F value Pr(>F)
A1vsA2     1  1.774  1.7737   2.132  0.178
A3vsA1A2   1  0.147  0.1467   0.176  0.684
Residuals  9  7.487  0.8319               

Error: id:BvsC
              Df Sum Sq Mean Sq F value  Pr(>F)   
BvsC           1 10.005  10.005  12.670 0.00612 **
A1vsA2:BvsC    1  0.339   0.339   0.429 0.52895   
A3vsA1A2:BvsC  1  0.132   0.132   0.167 0.69206   
Residuals      9  7.107   0.790                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: id:YvsZ
              Df Sum Sq Mean Sq F value Pr(>F)
YvsZ           1  0.061  0.0606   0.096  0.764
A1vsA2:YvsZ    1  0.793  0.7927   1.250  0.293
A3vsA1A2:YvsZ  1  0.000  0.0000   0.000  0.996
Residuals      9  5.708  0.6342               

Error: id:BvsC:YvsZ
                   Df Sum Sq Mean Sq F value   Pr(>F)    
BvsC:YvsZ           1  64.67   64.67  95.985 4.24e-06 ***
A1vsA2:BvsC:YvsZ    1   1.37    1.37   2.034    0.188    
A3vsA1A2:BvsC:YvsZ  1   0.04    0.04   0.064    0.806    
Residuals           9   6.06    0.67                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

>anova(lmerMdl1)

fixed-effect model matrix is rank deficient so dropping 4 columns / coefficients
Analysis of Variance Table of type III  with  Satterthwaite 
approximation for degrees of freedom
                   Sum Sq Mean Sq NumDF DenDF F.value   Pr(>F)    
A1vsA2              1.491   1.491     1                            
A3vsA1A2            0.123   0.123     1    12   0.235 0.6364995    
BvsC               10.005  10.005     1    36  19.079 0.0001017 ***
YvsZ                0.061   0.061     1    36   0.116 0.7358099    
A1vsA2:BvsC         0.339   0.339     1                            
A3vsA1A2:BvsC       0.132   0.132     1    36   0.252 0.6187426    
A1vsA2:YvsZ         0.793   0.793     1                            
A3vsA1A2:YvsZ       0.000   0.000     1    36   0.000 0.9956444    
BvsC:YvsZ          64.667  64.667     1    36 123.316 3.515e-13 ***
A1vsA2:BvsC:YvsZ    1.370   1.370     1                            
A3vsA1A2:BvsC:YvsZ  0.043   0.043     1    36   0.082 0.7765143    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Questions:

  • Why do the degrees of freedom differ so widely between the two procedures?
  • Why are some values not being displayed in the lmerMdl result?
  • Do I have to care about the warnings (e.g. model matrix rank is deficient)?

UPDATE:

In the comments, @amoeba pointed out that one should use (BvsC + YvsZ | id) instead of (1 | id). This makes the dfs smaller.

Apart from that, I figured out how to display all the values from the anova table by explicitly excluding the interaction between the two orthogonal variables (A1vsA2:A3vsA1A2). This also gets rid of the warning:

lmerMdl2 <- lmer(Y ~ A1vsA2*A3vsA1A2*BvsC*YvsZ + (BvsC + YvsZ | id),   data = XL, REML = FALSE)

lmerMdl3 <- lmer(Y ~ A1vsA2 + A3vsA1A2 + BvsC + YvsZ +
                     A1vsA2:BvsC + A3vsA1A2:BvsC + A1vsA2:YvsZ + A3vsA1A2:YvsZ + BvsC:YvsZ +
                     A1vsA2:BvsC:YvsZ + A3vsA1A2:BvsC:YvsZ + (BvsC + YvsZ | id), data = XL, REML = FALSE)

> anova(lmerMdl2)
                   Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
A1vsA2              0.894   0.894     1                             
A3vsA1A2            0.082   0.082     1 13.358   0.217  0.648891    
BvsC                6.342   6.342     1 12.043  16.869  0.001443 ** 
YvsZ                0.041   0.041     1 15.093   0.110  0.744851    
A1vsA2:BvsC         0.417   0.417     1                             
A3vsA1A2:BvsC       0.084   0.084     1 12.043   0.223  0.645364    
A1vsA2:YvsZ         0.540   0.540     1                             
A3vsA1A2:YvsZ       0.000   0.000     1 15.093   0.000  0.995795    
BvsC:YvsZ          64.667  64.667     1 24.000 171.995 1.945e-12 ***
A1vsA2:BvsC:YvsZ    1.370   1.370     1                             
A3vsA1A2:BvsC:YvsZ  0.043   0.043     1 24.000   0.114  0.738474   

> anova(lmerMdl3)
                   Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
A1vsA2              0.986   0.986     1 13.358   2.623  0.128673    
A3vsA1A2            0.082   0.082     1 13.358   0.217  0.648891    
BvsC                6.342   6.342     1 12.043  16.869  0.001443 ** 
YvsZ                0.041   0.041     1 15.093   0.110  0.744851    
A1vsA2:BvsC         0.215   0.215     1 12.043   0.571  0.464417    
A3vsA1A2:BvsC       0.084   0.084     1 12.043   0.223  0.645364    
A1vsA2:YvsZ         0.540   0.540     1 15.093   1.436  0.249203    
A3vsA1A2:YvsZ       0.000   0.000     1 15.093   0.000  0.995795    
BvsC:YvsZ          64.667  64.667     1 24.000 171.995 1.945e-12 ***
A1vsA2:BvsC:YvsZ    1.370   1.370     1 24.000   3.644  0.068306 .  
A3vsA1A2:BvsC:YvsZ  0.043   0.043     1 24.000   0.114  0.738474

Still, I don't understand why the lmer function returns from (1, 13) to (1, 24) dfs for the within-subject factors and (1, 12) for the between-subject whereas aov always returns (1, 11) for each factor no matter if it is within or between?

UPDATE 2

It seems to work correctly using Kenward-Roger approximation instead of the Satterthwaite one:

> anova(lmerMdl3, ddf = "Kenward-Roger")
                   Sum Sq Mean Sq NumDF  DenDF F.value    Pr(>F)    
A1vsA2              0.740   0.740     1     9   1.967  0.194269    
A3vsA1A2            0.061   0.061     1     9   0.163  0.696099    
BvsC                4.757   4.757     1     9  12.652  0.006149 ** 
YvsZ                0.031   0.031     1     9   0.082  0.780568    
A1vsA2:BvsC         0.161   0.161     1     9   0.428  0.529244    
A3vsA1A2:BvsC       0.063   0.063     1     9   0.167  0.692268    
A1vsA2:YvsZ         0.405   0.405     1     9   1.077  0.326382    
A3vsA1A2:YvsZ       0.000   0.000     1     9   0.000  0.996399    
BvsC:YvsZ          48.500  48.500     1     9 128.996 1.229e-06 ***
A1vsA2:BvsC:YvsZ    1.028   1.028     1     9   2.733  0.132693    
A3vsA1A2:BvsC:YvsZ  0.032   0.032     1     9   0.086  0.776523
Improve this question
$\endgroup$
19
  • $\begingroup$ Your lmer call is very different from your aov call because you are not specifying the within-subject factors. The standard way to do it would be to use (BvsC*YvsZ | id) instead of (1 | id). $\endgroup$
    – amoeba
    Commented Jan 2, 2018 at 23:11
  • $\begingroup$ @amoeba It returns an error: Error: number of observations (=48) <= number of random effects (=48) for term (BvsC * YvsZ | id); the random-effects parameters and the residual variance (or scale parameter) are probably unidentifiable $\endgroup$
    – mat
    Commented Jan 3, 2018 at 14:54
  • $\begingroup$ Ah yes, this is because you have only 1 observation per within-subject factor combination. Then (BvsC + YvsZ | id). $\endgroup$
    – amoeba
    Commented Jan 3, 2018 at 19:13
  • $\begingroup$ @amoeba The analysis runs without error. However I don't understand why the lmer functions returns (1, 36) dfs for the within-subject factors and (1, 12) for the between-subject whereas aov always returns (1, 11) for each factor no matter if it is within or between? $\endgroup$
    – mat
    Commented Jan 4, 2018 at 8:07
  • 1
    $\begingroup$ You don't have to, but it is my recommendation to have 1 specific focused question per thread. If you have several distinct questions, you will increase the chances of getting good answers by asking them separately. $\endgroup$
    – amoeba
    Commented Jan 5, 2018 at 14:19

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.