precise estimates from emmeans across lm and aov_ez

Home forums ANOVA precise estimates from emmeans across lm and aov_ez

This topic contains 3 replies, has 2 voices, and was last updated by  henrik 3 weeks ago.

  • 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)
    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!

  • #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’
    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.

  • #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!

You must be logged in to reply to this topic.

 

Author: klhenry