3

I would like to fit a one-way between-subject anova that assumes unequal variances between groups.

Reproducible example:

library(emmeans)
library(car)

set.seed(123)
n <- 50
DF <- data.frame(score = c(rnorm(n, sd = 10), rnorm(n, sd = 30), rnorm(n, sd = 40)),
                 treatment = rep(c("A", "B", "C"), each = n),
                 subject = 1:(n*3))


leveneTest(score ~ treatment, DF) # Shows heterogeneity of variance


mdl <- lm(score ~ treatment, data = DF)
emmeans(mdl, ~treatment) # same SE for all the means


treatment  emmean   SE  df lower.CL upper.CL


A           0.344 3.99 147    -7.54     8.23


B           4.392 3.99 147    -3.50    12.28


C         -10.156 3.99 147   -18.04    -2.27


Confidence level used: 0.95

Is there a way to tweak lm (or lmer) to take into account unequal variance?

kjetil b halvorsen
  • 77,844
mat
  • 551

2 Answers2

6

You can do heterogeneous variance in a variety of ways in R. A simple way is through the gls package

library(nlme)
mod = gls(score~treatment, data=DF, 
          weights = varIdent(form = ~1|treatment), 
              method="ML")
emmeans(mod, ~treatment)

although lme4 is more efficient and popular, nlme offers a variety of structures for the residuals. Of course you can always go the Bayesian route if you need even more flexibility

library(brms)
modb <- brm(
     bf(score ~ treatment,
        sigma ~ treatment), 
      family = gaussian,
      data=DF)
emmeans(modb, ~ treatment)
kjetil b halvorsen
  • 77,844
bdeonovic
  • 10,127
  • Thanks! What about lmer, is there an alternative? I am used to that function. – mat Aug 15 '22 at 18:31
  • 1
    An alternative is the R package gamlss For an example see https://stats.stackexchange.com/questions/495811/are-there-better-approaches-than-the-weighted-mean/495845#495845 – kjetil b halvorsen Aug 15 '22 at 18:56
  • 1
    we can read in an lme4 vignette: “The main advantage of nlme relative to lme4 is a user interface for fitting models with structure in the residuals (various forms of heteroscedasticity and autocorrelation) and in the random-effects covariance matrices (e.g., compound symmetric models). – bdeonovic Aug 15 '22 at 20:47
0

Found a way to do it via lmer:

library(lme4)

mdl <- lmer(score ~ treatment + (0 + treatment|subject), data = DF,
              + control = lmerControl(check.nobs.vs.nRE = "ignore", 
                                      check.nobs.vs.nlev = "ignore"))
#Warning messages:


1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :


unable to evaluate scaled gradient


2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :


Model failed to converge: degenerate  Hessian with 3 negative eigenvalues


emmeans(mdl, ~ treatment)


treatment  emmean   SE df lower.CL upper.CL


A           0.344 1.31 49    -2.29     2.98


B           4.392 3.84 49    -3.33    12.11


C         -10.156 5.60 49   -21.40     1.09


Degrees-of-freedom method: kenward-roger


Confidence level used: 0.95
kjetil b halvorsen
  • 77,844
mat
  • 551