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

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<.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.)

• Why are you trying different packages, functions, libraries etc in R? What are trying to achieve? It would be more productive to explain what analysis you are aiming for and to show us your one best attempt at implementing this analysis. Jul 22 at 5:51
• I second @dipetkov 's point about needing to understand your goal here. It would also be helpful if you could explain what you do understand about differences between these methods/functions.
– mkt
Jul 22 at 6:45
• this is my first time running ancova/anovas in R and i am new to it. I wanted to double check my work by trying a different library to see if results are the same. I found a lot of resources doing things differently and I wasnt sure what is most appropriate and if i am doing something wrong to get different results. I dont fully grasp the math either, sorry. I am trying to do the same thing in all of them. run an ancovoa betwwen group, within hemi, and etiv as covar.
– matt
Jul 22 at 19:56
• i THINK i understand that lme/lmer fits it into an anova after (not sure what that means). I read also that it might be best as my data is unbalanced. But not sure. and still, i dont understand then why the results between C and D would differ
– matt
Jul 22 at 19:58

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.

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.