ANCOVA vs linear mixed model in R, using different functions gives me different results

Question

I have tried different packages, functions, libraries etc in R and my results differ, whether it be p values or F values or DF etc. This is especially true for my within variable (hemi)..

between: group (Old adults, younger adults) unbalanced, ~150 young, 100 old within: brain hemisphere (left vs right) dv= volume (volume of left and right hippocampus brain region) covariate = eitv(whole brain volume)

here is an example:

sub group   etiv    hemi    volume
sub-169 OA  1740000 left_vol    120
sub-169 OA  1740000 right_vol   134
sub-198 YA  1730000 left_vol    183

I dont understand why they differ. Below are the different functions/code I used. I would really appreciate some insight as to what might is most appropriate and why I might be getting different results.

A
cova.aov2 <- anova_test(
  data = anova_long, dv = volume, wid = sub, between = group, within = hemi, covariate = etiv, type=3, effect.size = "ges", detailed=TRUE)
ANOVA Table (type III tests)
   Effect DFn DFd        SSn       SSd         F        p p&lt;.05      ges

1 (Intercept)   1 216 202437.268 1082993.4 40.376000 1.22e-09     * 1.39e-01
2        etiv   1 216   8281.849 1082993.4  1.652000 2.00e-01       7.00e-03
3       group   1 216 127405.844 1082993.4 25.411000 9.79e-07     * 9.20e-02
4        hemi   1 216   4035.160  167664.3  5.198000 2.40e-02     * 3.00e-03
5   etiv:hemi   1 216      0.172  167664.3  0.000222 9.88e-01       1.38e-07
6  group:hemi   1 216    737.333  167664.3  0.950000 3.31e-01       5.89e-04

B

my_aov = aov(data=anova_long, formula = volume~group*hemi+etiv+Error(sub/(hemi)))
summary(my_aov)
Error: sub
           Df  Sum Sq Mean Sq F value  Pr(>F)    
group       1  129464  129464  25.821 8.1e-07 ***
etiv        1    8282    8282   1.652     0.2    
Residuals 216 1082993    5014                    

Error: sub:hemi
            Df Sum Sq Mean Sq F value Pr(>F)    
hemi         1 317501  317501 410.926 <2e-16 ***
group:hemi   1    739     739   0.956  0.329    
Residuals  217 167664     773                   

C

m1<-lme(volume~etiv+group*hemi, random=~1|sub, method="REML",data=ancova_long)
result_m1=anova(m1)
> result_m1
            numDF denDF  F-value p-value
(Intercept)     1   217 4606.777  <.0001
etiv            1   216    2.062  0.1524
group           1   216   25.411  <.0001
hemi            1   217  410.927  <.0001
group:hemi      1   217    0.956  0.3293

D

lmer_m = lmer(volume~etiv+group*hemi+(1|sub), data=ancova_long)
result2=anova(lmer_m)

Type III Analysis of Variance Table with Satterthwaite's method
           Sum Sq Mean Sq NumDF DenDF  F value    Pr(>F)    
etiv         1276    1276     1   216   1.6518    0.2001    
group       19634   19634     1   216  25.4107 9.791e-07 ***
hemi       274814  274814     1   217 355.6785 < 2.2e-16 ***
group:hemi    739     739     1   217   0.9560    0.3293    

Warning messages:
1: Some predictor variables are on very different scales: consider rescaling 
2: Some predictor variables are on very different scales: consider rescaling 
i get same result using: lmeretiv=lmer(volume~etiv+group*hemi +(1|sub:group)+(1:sub:group:hemi), data=ancova_long)

E

X<- aov_car(formula = volume~group*hemi+etiv+Error(sub/(hemi)), factorize = FALSE, observed = c("etiv", "group"),type=3, print.formula=TRUE,data = ancova_long) #anova_table = list(correction = "none", es = "pes", p_adjust_method="holm)) #with holm, only group is sig. 
Anova Table (Type 3 tests)

Response: volume
      Effect     df     MSE         F   ges p.value
1      group 1, 216 5013.86 25.41 ***  .092   <.001
2       etiv 1, 216 5013.86      1.65  .006    .200
3       hemi 1, 216  776.22    5.20 *  .003    .024
4 group:hemi 1, 216  776.22      0.95 <.001    .331
5  etiv:hemi 1, 216  776.22      0.00 <.001    .988
(i can reproduce this using aov_4 and aov_ex of course

F

mixy<-mixed(volume~etiv+group*hemi+(1|sub), anova_longfact, method = "S")
> mixy
Mixed Model Anova Table (Type 3 tests, S-method)

Model: volume ~ etiv + group * hemi + (1 | sub)
Data: anova_longfact
      Effect        df          F p.value
1       etiv 1, 216.00       1.65    .200
2      group 1, 216.00  25.41 ***   <.001
3       hemi 1, 217.00 355.68 ***   <.001
4 group:hemi 1, 217.00       0.96    .329

(note sometimes i refer to a different data sheet, but they are the same.)

Answer 1

To summarise: You've discovered that different packages & functions may implement different methods, so you shouldn't expect to get the same results when you throw your data at a sundry collection of R functions.

A more efficient approach is to start reading the documentation pages for the R methods you want to use and understand. In R and RStudio you can access the documentation with the help() function or the ? operator. For example, type help("lmer") or ?lmer to learn more about how you can use lmer to fit linear mixed-effects models. If you come across terms or theory that you don't understand, look them up. Or check here on Cross Validated. The difference from your current approach is that you'll have a very specific question and it's more likely you'll get a satisfactory answer.

Now with this caveat aside, let's look at why C and D give different results. It requires delving into a couple of advanced topics.

The first important point is about statistics: to do ANOVA analysis and compare models, you should fit the models using maximum likelihood (ML), not restricted maximum likelihood (REML).

library("nlme")
library("lme4")
Use REML to fit two mixed-effects models.
m1 <- lme(distance ~ age * Sex, random = ~ 1 | Subject, method = "ML", data = Orthodont)
m2 <- lmer(distance ~ age * Sex + (1 | Subject), REML = FALSE, data = Orthodont)
anova(m1)
#>             numDF denDF  F-value p-value
#> (Intercept)     1    79 4288.082  <.0001
#> age             1    79  120.900  <.0001
#> Sex             1    25    9.664  0.0046
#> age:Sex         1    79    6.223  0.0147
anova(m2)
#> Analysis of Variance Table
#>         npar  Sum Sq Mean Sq  F value
#> age        1 235.356 235.356 125.5485
#> Sex        1  18.813  18.813  10.0359
#> age:Sex    1  12.114  12.114   6.4622

So our first ANOVA tables are close but not quite the same.

This brings us to the second advanced topic, which is about R. anova is a generic function: it calls other functions which implement ANOVA analysis for specific model types. This process is called "method dispatch".

So even though we typed anova(m1) and anova(m2), we ended up calling two different anova implementations.

?anova.lme # Does ANOVA on `nlme::lme` models.
?anova.merMod # Does ANOVA on `lme4` models, incl. `lmer` and `glmer`

By reading the docs I learned that by default anova.lme does some adjustment to the residual standard errors. We want to turn the adjustment off for our comparisons of ANOVAs.

anova(m1, adjustSigma = FALSE)
#>             numDF denDF  F-value p-value
#> (Intercept)     1    79 4453.008  <.0001
#> age             1    79  125.550  <.0001
#> Sex             1    25   10.035   0.004
#> age:Sex         1    79    6.462   0.013
anova(m2)
#> Analysis of Variance Table
#>         npar  Sum Sq Mean Sq  F value
#> age        1 235.356 235.356 125.5485
#> Sex        1  18.813  18.813  10.0359
#> age:Sex    1  12.114  12.114   6.4622

Finally we have (almost) the same ANOVAs, from the different methods to fit mixed-effects models. I think the remaining numeric differences are due to the details of the implementations behind lme4 and nlme.

Lesson learned: Always start with the documentation. It will save you lots of time.

Answer 2

This is mostly an extended comment on the answer from @dipetkov (+1), which I suggest that you accept.

First, note that you specified two different sets of interaction terms in your models. Some added an etiv:hemi interaction to the group:hemi interaction that all models included. You might not have thought you were asking for it in those cases, but the way you formulated your models did. Those models should differ from the others. They seem to be the models in which the p-value for hemi on its own differs from the others (0.02 with that extra interaction, versus < 0.0001 without).

Second, consider why you are requesting Type III sums of squares for your evaluations. See this page for discussion of different types of sums of squares.The R car package warns in the help page for its Anova() function: "Be very careful in formulating the model for type-III tests, or the hypotheses tested will not make sense."

Third, there is not universal agreement on how to estimate degrees of freedom for mixed models. See this page and its links for some discussion, and why lme4 as a default does not provide p values.

Fourth, most results of your models don't seem to be that different. For example, the p-value for the group:hemi interaction seems to be $0.33 \pm 0.001$ for all models.