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?

Improve this question
$\endgroup$
1

2 Answers 2

Reset to default
7
$\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)
Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks! What about lmer, is there an alternative? I am used to that function. $\endgroup$
    – mat
    Commented 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$
    – kjetil b halvorsen ♦
    Commented 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
    Commented 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
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ I'm not convinced this fits a model for unequal variances. If it does, can you explain? $\endgroup$
    – Russ Lenth
    Commented Oct 2, 2022 at 18:26
  • $\begingroup$ Looks like maybe based on something like this: rpubs.com/bbolker/factorvar $\endgroup$
    – bdeonovic
    Commented Jul 1 at 15:04

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.