Home forums ANOVA precise estimates from emmeans across lm and aov_ez

• Author
Posts
• #229
klhenry
Participant

I am comparing the results of emmeans across different models with a toy dataset. I am fitting a model in which a continuous outcome is predicted by a 3-level factor and a control variable. I have fit the model in two ways, with a linear regression model using the lm function in R (representing the factor my two dummy codes). And as an ANCOVA using the aov_ez function. The emmeans that I obtain from the lm and ancova model are not exactly the same, although I expected them to be. Likewise, the contrasts from the ANCOVA model that should correspond with the dummy coded effects aren’t precisely the same. Is there a reason these would be slightly different? Do I need to modify my aov_ez model in some way to make the emmeans results and contrasts equivalent? I do get the exact same results if there is no covariate.

Linear model:
lm.1 <- lm(sleep ~ cond.f + anxiety.m, data = slp)
summary(lm.1)

compare <- emmeans(lm.1, ~ cond.f + anxiety.m)
compare

ANCOVA:
anvmod <- aov_ez(id = “id”, dv = “sleep”, data = slp, between = c(“cond.f”),
covariate = c(“anxiety.m”),
factorize = FALSE)
summary(anvmod)

compare <- lsmeans(anvmod, ~ cond.f + anxiety.m)
compare

# make constrasts
cond2_1 <- c(-1, 1, 0)
cond3_1 <- c(-1, 0, 1)
summary(contrast(compare, list(cond2_1, cond3_1)))

data are here: https://pastebin.com/qTajz9ak

• #233
henrik
Keymaster

This is indeed a bug in `afex` with ANCOVAs. It has to do with the way the `aov` object is used internally and I am working on a fix. So long, one can circumvent the use of the `aov` object and simply use the `lm` object. This somewhat non-intuitively is done by setting `model = "multivariate"` in the call to `emmeans`:

``````library(afex)

lm.1 <- lm(sleep ~ cond.f + anxiety.m, data = slp)
emmeans(lm.1, ~ cond.f + anxiety.m)
#  cond.f             anxiety.m   emmean        SE  df lower.CL upper.CL
#  group-based     4.746059e-17 69.27959 0.7289347 596 67.84800 70.71119
#  group + partner 4.746059e-17 76.35960 0.7279758 596 74.92989 77.78931
#  self help       4.746059e-17 60.99081 0.7293691 596 59.55836 62.42325
#
# Confidence level used: 0.95

anvmod <- aov_ez(id = "id", dv = "sleep", data = slp, between = c("cond.f"),
covariate = c("anxiety.m"), factorize = FALSE)
emmeans(anvmod, ~ cond.f + anxiety.m, model = "multivariate")
#  cond.f             anxiety.m   lsmean        SE  df lower.CL upper.CL
#  group-based     4.746059e-17 69.27959 0.7289347 596 67.84800 70.71119
#  group + partner 4.746059e-17 76.35960 0.7279758 596 74.92989 77.78931
#  self help       4.746059e-17 60.99081 0.7293691 596 59.55836 62.42325
#
# Confidence level used: 0.95``````

I will post here again, when I have fixed the bug in `afex`. Thanks a lot for reporting it!

• #239
henrik
Keymaster

The development version of `afex` from `github` now returns the correct results. I will update the post when I submit this version to `CRAN`.

``````devtools::install_github("singmann/afex")
library(afex)
packageVersion("afex")
# [1] ‘0.21.0’

lm.1 <- lm(sleep ~ cond.f + anxiety.m, data = slp)
emmeans(lm.1, ~ cond.f + anxiety.m)
#  cond.f             anxiety.m   emmean        SE  df lower.CL upper.CL
#  group-based     4.746059e-17 69.27959 0.7289347 596 67.84800 70.71119
#  group + partner 4.746059e-17 76.35960 0.7279758 596 74.92989 77.78931
#  self help       4.746059e-17 60.99081 0.7293691 596 59.55836 62.42325
#
# Confidence level used: 0.95

anvmod <- aov_ez(id = "id", dv = "sleep", data = slp, between = c("cond.f"),
covariate = c("anxiety.m"), factorize = FALSE)
emmeans(anvmod, ~ cond.f + anxiety.m)
#  cond.f             anxiety.m   emmean        SE  df lower.CL upper.CL
#  group-based     4.746059e-17 69.27959 0.7289347 596 67.84800 70.71119
#  group + partner 4.746059e-17 76.35960 0.7279758 596 74.92989 77.78931
#  self help       4.746059e-17 60.99081 0.7293691 596 59.55836 62.42325
#
# Confidence level used: 0.95

emmeans(anvmod, ~ cond.f + anxiety.m, model = "multivariate")
#  cond.f             anxiety.m   emmean        SE  df lower.CL upper.CL
#  group-based     4.746059e-17 69.27959 0.7289347 596 67.84800 70.71119
#  group + partner 4.746059e-17 76.35960 0.7279758 596 74.92989 77.78931
#  self help       4.746059e-17 60.99081 0.7293691 596 59.55836 62.42325
#
# Confidence level used: 0.95 ``````

Thanks again for reporting this.

• #272
henrik
Keymaster

A new version of `afex`, version `0.21-2`, that fixes this bug, has just appeared on `CRAN`: https://cran.r-project.org/package=afex

I apologize for the delay!