r - 2-way Repeated Measures ANOVA - different degree of freedom among ezAnova, aov and lme [closed]

Question

I would like to run a simple repeated measures ANOVA with 2 factors (Food & Drink).

Let's have a dummy df:

nSubj = 10; nFood = 5; nDrnk = 2;
dummy <- data.frame(Subj = factor(rep(1:nSubj, each=nFood*nDrnk, time=1)),
                Food = factor(rep(letters[1:nFood],  each=nDrnk, time=nSubj)),
                Drnk = factor(rep(c("water","beer"), each=1, time=nSubj*nFood)),  
                Like = sample(0:10, size=nSubj*nFood*nDrnk, replace=T))

Using ezANOVA:

library(ez)
m1 <- ezANOVA(data=dummy, dv=Like, wid=Subj, within=c(Food, Drnk), detailed=F, type=3); m1

     Effect DFn DFd         F          p p<.05        ges
2      Food   4  36 0.3282034 0.85717626       0.01301648
3      Drnk   1   9 7.2461652 0.02472067     * 0.07354474
4 Food:Drnk   4  36 1.5834181 0.19974760       0.07367973

Using lme:

library(nlme)
m2 <- lme(Like ~ Food*Drnk, random = ~1|Subj/Food/Drnk, data=dummy, method = "ML"); anova(m2)

            numDF denDF   F-value p-value
(Intercept)     1    45 288.81124  <.0001
Food            4    36   0.29673  0.8781
Drnk            1    45   7.14446  0.0104
Food:Drnk       4    45   1.78966  0.1476

Using aov:

m3 <- aov(Like ~ Food*Drnk + Error(Subj/(Food*Drnk)),data=dummy); summary(m3)

Error: Subj
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  9  72.45    8.05               

Error: Subj:Food
          Df Sum Sq Mean Sq F value Pr(>F)
Food       4   10.9   2.725   0.328  0.857
Residuals 36  298.9   8.303               

Error: Subj:Drnk
          Df Sum Sq Mean Sq F value Pr(>F)  
Drnk       1  65.61   65.61   7.246 0.0247 *
Residuals  9  81.49    9.05                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Error: Subj:Food:Drnk
          Df Sum Sq Mean Sq F value Pr(>F)
Food:Drnk  4   65.7   16.43   1.583    0.2
Residuals 36  373.7   10.38  

The denDF in the lme is different: why? is my lme model correct?

EDIT: Since nSubj = 10; nFood = 5; nDrnk = 2 and it is a repeated measures, the dof should be: Food=(4,36); Drnk=(1,9); Food:Drnk=(4,36). Right? So, why does lme have different dof?

EDIT 2: following Sal's suggestion, I run aov with another error:

m4 <- aov(Like ~ Food*Drnk + Error(Subj/Food/Drnk),data=dummy); summary(m4)

Error: Subj
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  9  100.4   11.16               

Error: Subj:Food
          Df Sum Sq Mean Sq F value Pr(>F)
Food       4  10.74   2.685   0.319  0.863
Residuals 36 302.66   8.407               

Error: Subj:Food:Drnk
      Df Sum Sq Mean Sq F value  Pr(>F)   
Drnk       1   54.8   54.76   7.267 0.00984 **
Food:Drnk  4   65.1   16.29   2.161 0.08875 . 
Residuals 45  339.1    7.54     

This is interesting, the results are the same of m2 with the Drnk variable ignored (?).

Next, Changing the order of the variables in the lme's random:

m5 <- lme(Like ~ Food*Drnk, random = ~1|Subj/Drnk/Food, data=dummy, method = "ML"); anova(m5)
            numDF denDF   F-value p-value
(Intercept)     1    72 117.06035  <.0001
Food            4    72   0.61245  0.6550
Drnk            1     9   0.36079  0.5629
Food:Drnk       4    72   0.21603  0.9287

The DenDFs of both Food and interaction are wrong... other combinations of lme's random are either wrong or gives error. Therefore, how can I model a 2-way Repeated Measures ANOVA using lme?

Any help is appreciated

Cheers

Answer 1

Answer mostly from comments.

The model you are seeking with the nlme package would have crossed random effects. Pinheiro and Bates (2000) discuss using the lme function with crossed random effects by invoking the pdBlocked and pdIdent objects. It is not obvious how this is implemented for different designs, and, personally, I've had trouble implementing this approach following Pinheiro and Bates or solutions described online.

One approach is to use the lme4 package, with which it is apparently easier to specify crossed random effects. In the example below, the degrees of freedom from the lmer model are similar to those from the aov function.

However, you should be cautious expecting the results from a mixed effects model to be similar to ordinary least squares (OLS) tests.

References Pinheiro, J.C., and D.M. Bates. 2000. Mixed -Effects Models in S and S-Plus.

set.seed(1234)

nSubj = 10; nFood = 5; nDrnk = 2;
dummy <- data.frame(Subj = factor(rep(1:nSubj, each=nFood*nDrnk, time=1)),
                Food = factor(rep(letters[1:nFood],  each=nDrnk, time=nSubj)),
                Drnk = factor(rep(c("water","beer"), each=1, time=nSubj*nFood)),  
                Like = sample(0:10, size=nSubj*nFood*nDrnk, replace=T))

m4 <- aov(Like ~ Food*Drnk + Error(Subj/(Food*Drnk)),data=dummy)
summary(m4)

   ### Error: Subj
   ###           Df Sum Sq Mean Sq F value Pr(>F)
   ### Residuals  9  39.69    4.41               
   ### 
   ### Error: Subj:Food
   ###           Df Sum Sq Mean Sq F value Pr(>F)
   ### Food       4  29.04   7.260   0.854  0.501
   ### Residuals 36 306.16   8.504               
   ###
   ### Error: Subj:Drnk
   ###           Df Sum Sq Mean Sq F value  Pr(>F)   
   ### Drnk       1  86.49   86.49   18.35 0.00204 **
   ### Residuals  9  42.41    4.71
   ###
   ### Error: Subj:Food:Drnk
   ###           Df Sum Sq Mean Sq F value Pr(>F)
   ### Food:Drnk  4   57.2   14.29   1.412   0.25
   ### Residuals 36  364.4   10.12                     

library(lme4)
library(lmerTest)

m7 = lmer(Like ~ Food*Drnk + (1|Subj) + (1|Food:Subj) + (1|Drnk:Subj), data=dummy)
anova(m7, ddf="Kenward-Roger")

   ### Type III Analysis of Variance Table with Kenward-Roger's method
   ###           Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
   ### Food       29.04    7.26     4    36  0.8681 0.49247  
   ### Drnk       86.49   86.49     1     9 10.3416 0.01056 *
   ### Food:Drnk  57.16   14.29     4    36  1.7086 0.16947

ranova(m7)