3
$\begingroup$

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?

$\endgroup$
1

2 Answers 2

6
$\begingroup$

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

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   
$\endgroup$
1
  • $\begingroup$ I'm not convinced this fits a model for unequal variances. If it does, can you explain? $\endgroup$
    – Russ Lenth
    Oct 2, 2022 at 18:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.