precise estimates from emmeans across lm and aov_ez

[klhenry]

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

[henrik]

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)
slp <- read.csv("https://pastebin.com/raw/qTajz9ak")

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!

[henrik]

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’
slp <- read.csv("https://pastebin.com/raw/qTajz9ak")

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.

[henrik]

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!

[Author]

Posts