# Anova with unequal variances via regression modeling

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?

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)

• Thanks! What about lmer, is there an alternative? I am used to that function.
– mat
Commented Aug 15, 2022 at 18:31
• An alternative is the R package gamlss For an example see stats.stackexchange.com/questions/495811/… Commented Aug 15, 2022 at 18:56
• 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). Commented Aug 15, 2022 at 20:47

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

• I'm not convinced this fits a model for unequal variances. If it does, can you explain? Commented Oct 2, 2022 at 18:26
• Looks like maybe based on something like this: rpubs.com/bbolker/factorvar Commented Jul 1 at 15:04