Satterthwaite degrees of freedom in lmerTest for a 3x2x2 mixed model vs. mixed ANOVA dfs

Question

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

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)`. — amoeba, Jan 02 '18 at 23:11
@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` — mat, Jan 03 '18 at 14:54
Ah yes, this is because you have only 1 observation per within-subject factor combination. Then `(BvsC + YvsZ | id)`. — amoeba, Jan 03 '18 at 19:13
@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? — mat, Jan 04 '18 at 08:07
Could you update your question with the corresponding output? — amoeba, Jan 04 '18 at 09:32
@amoeba for lmer dfs were not actually 36, but they vary from 13 to 24 (see update) — mat, Jan 04 '18 at 15:28
Looking at your bullet list of questions, the one about warnings is not relevant anymore because you don't get any with `lmerMdl3`, right? The Q about REML is very different from everything else and is covered elsewhere on this site - make a search. This leaves two Qs: one about effect sizes and one about dfs. I suggest you separate them. The Q about effect sizes is not very clear (what effect size in `aov` you want it to be comparable to?). Do you maybe want to center *this* Q on dfs? — amoeba, Jan 05 '18 at 13:17
@amoeba If I had to choose one question to solve, it would be the one about the dfs — mat, Jan 05 '18 at 13:46
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. — amoeba, Jan 05 '18 at 14:19
This https://stats.stackexchange.com/questions/84268 looks relevant. Would be interesting to run `anova(lmerMdl3, ddf="Kenward-Roger")` and see what comes out, but AFAIK this can be much slower. — amoeba, Jan 05 '18 at 17:19
OK, my current understanding based on BenBolker's answer in the linked thread (see also my comments under his answer) is that `lmerTest` for your `lmerMdl3` substantially overestimates some of the within-subject dfs. — amoeba, Jan 07 '18 at 13:47
Update: but do try Kenward-Roger. It performs correctly in the example in the linked thread. — amoeba, Jan 07 '18 at 15:24
OK, I tried K-R myself and it ousputs 9 dfs throughout. I think this settles it. — amoeba, Jan 07 '18 at 15:51
BTW, why did you say that aov returns (1, 11) dfs? Isn't it (1, 9)? — amoeba, Jan 07 '18 at 15:51
@amoeba indeed, but that means it's 1,11 for each factor; I tried the K-R dfs, but the dfs are not 9 (I updated my post) — mat, Jan 08 '18 at 13:51
You should add `ddf = "Kenward-Roger"` to the `anova` call, not to the `lmer` call. — amoeba, Jan 08 '18 at 13:56
Can you explan why 9 in the `aov` output means 11 for each factor? — amoeba, Jan 08 '18 at 13:57
@amoeba My bad, (1, 9) indeed the correct *dfs*! For some reason, when I run the same analysis with my real data base (256 rows in wide format; 1024 rows in long format) the *dfs* get messed up. For example, the double interaction A3vsA1A2:BvsC:YvsZ *dfs*' are (3, 762). The only difference between the simulated data and the real data being that groups (A1 vs. A2 vs. A3) are not exactly balanced. — mat, Jan 09 '18 at 15:04

Satterthwaite degrees of freedom in lmerTest for a 3x2x2 mixed model vs. mixed ANOVA dfs

0 Answers0

Linked