It’s actually not very difficult to handle heteroscedasticity in simple linear models (e.g., one- or two-way ANOVA-like models).
Robustness of ANOVA
First, as others have note, the ANOVA is amazingly robust to deviations from the assumption of equal variances, especially if you have approximately balanced data (equal number of observations in each group). Preliminary tests on equal variances, other the other hand, are not (though Levene’s test is much better than the F-test commonly taught in textbooks). As George Box put it:
To make the preliminary test on variances is rather like putting to sea in a rowing boat to find out whether conditions are sufficiently calm for an ocean liner to leave port!
Even though the ANOVA is very robust, as it’s very easy to take heteroscedaticity into account, there’s little reason not to do so.
If you’re really interested in differences in means, the non-parametric tests (e.g., the Kruskal–Wallis test) are really not of any use. They do test differences between groups, but they do not in general test differences in means.
Let’s generate a simple example of data where one would like to use ANOVA, but where the assumption of equal variances is not true.
pop = data.frame(group=factor(c("A","B","C")),
d = do.call(rbind, rep(list(pop),13))
d$x = rnorm(nrow(d), d$mean, d$sd)
We have three groups, with (clear) differences in both means and variances:
stripchart(x ~ group, data=d)
Not surprisingly, a normal ANOVA handles this quite well:
> mod.aov = aov(x ~ group, data=d)
Df Sum Sq Mean Sq F value Pr(>F)
group 2 199.4 99.69 13.01 5.6e-05 ***
Residuals 36 275.9 7.66
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
So, which groups differ? Let’s use Tukey’s HSD method:
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = x ~ group, data = d)
diff lwr upr p adj
B-A 1.736692 -0.9173128 4.390698 0.2589215
C-A 5.422838 2.7688327 8.076843 0.0000447
C-B 3.686146 1.0321403 6.340151 0.0046867
With a P-value of 0.26, we can’t claim any difference (in means) between group A and B. And even if we didn’t take into account that we did three comparisons, we wouldn’t get a low P-value (P = 0.12):
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5098 0.7678 0.664 0.511
groupB 1.7367 1.0858 1.599 0.118
groupC 5.4228 1.0858 4.994 0.0000153 ***
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 2.768 on 36 degrees of freedom
Why is that? Based on the plot, there is a pretty clear difference. The reason is that ANOVA assumes equal variances in each group, and estimates a common standard deviation of 2.77 (shown as ‘Residual standard error’ in the
summary.lm table, or you can get it by taking the square root of the residual mean square (7.66) in the ANOVA table).
But group A has a (population) standard deviation of 1, and this overestimate of 2.77 makes it (needlessly) difficult to get statistically significant results, i.e., we have a test with (too) low power.
‘ANOVA’ with unequal variances
So, how to fit a proper model, one that takes into account the differences in variances? It’s easy in R:
> oneway.test(x ~ group, data=d, var.equal=FALSE)
One-way analysis of means (not assuming equal variances)
data: x and group
F = 12.7127, num df = 2.000, denom df = 19.055, p-value = 0.0003107
So, if you want to run a simple one-way ‘ANOVA’ in R without assuming equal variances, use this function. It’s basically an extension of the (Welch)
t.test() for two samples with unequal variances.
Unfortunately, it doesn’t work with
TukeyHSD() (or most other functions you use on
aov objects), so even if we’re pretty sure there are group differences, we don’t know where they are.
Modelling the heteroscedasticity
The best solution is to model the variances explicitly. And it’s very easy in R:
> mod.gls = gls(x ~ group, data=d,
weights=varIdent(form= ~ 1 | group))
Denom. DF: 36
numDF F-value p-value
(Intercept) 1 16.57316 0.0002
group 2 13.15743 0.0001
Still significant differences, of course. But now the differences between group A and B have also become statically significant (P = 0.025):
Generalized least squares fit by REML
Model: x ~ group
Structure: Different standard
deviations per stratum
Formula: ~1 | group
A B C
1.000000 2.444532 3.913382
Value Std.Error t-value p-value
(Intercept) 0.509768 0.2816667 1.809829 0.0787
groupB 1.736692 0.7439273 2.334492 0.0253
groupC 5.422838 1.1376880 4.766542 0.0000
Residual standard error: 1.015564
Degrees of freedom: 39 total; 36 residual
So using an appropriate model helps! Also note that we get estimates of the (relative) standard deviations. The estimated standard deviation for group A can be found at the bottom of the, results, 1.02. The estimated standard deviation of group B is 2.44 times this, or 2.48, and the estimated standard deviation of group C is similarly 3.97 (type
intervals(mod.gls) to get confidence intervals for the relative standard deviations of groups B and C).
Correcting for multiple testing
However, we really should correct for multiple testing. This is easy using the ‘multcomp’ library. Unfortunately, it does’t have built-in support for ‘gls’ objects, so we’ll have to add some helper functions first:
model.matrix.gls <- function(object, ...)
model.matrix(terms(object), data = getData(object), ...)
model.frame.gls <- function(object, ...)
model.frame(formula(object), data = getData(object), ...)
terms.gls <- function(object, ...)
Now let’s get to work:
> mod.gls.mc = glht(mod.gls, linfct = mcp(group = "Tukey"))
Estimate Std. Error z value Pr(>|z|)
B - A == 0 1.7367 0.7439 2.334 0.0480 *
C - A == 0 5.4228 1.1377 4.767 <0.001 ***
C - B == 0 3.6861 1.2996 2.836 0.0118 *
Still statistically significant difference between group A and group B! ☺ And we can even get (simultaneous) confidence intervals for the differences between group means:
Estimate lwr upr
B - A == 0 1.73669 0.01014 3.46324
C - A == 0 5.42284 2.78242 8.06325
C - B == 0 3.68615 0.66984 6.70245
Using an approximately (here exactly) correct model, we can trust these results!
Note that for this simple example, the data for group C doesn’t really add any information on the differences between group A and B, since we model both separate means and standard deviations for each group. We could have just used pairwise t-tests corrected for multiple comparisons:
> pairwise.t.test(d$x, d$group, pool.sd=FALSE)
Pairwise comparisons using t tests with non-pooled SD
data: d$x and d$group
B 0.03301 -
C 0.00098 0.02032
P value adjustment method: holm
However, for more complicated models, e.g., two-way models, or linear models with many predictors, using GLS (generalised least squares) and explicitly modelling the variance functions is the best solution.
And the variance function need not simply be a different constant in each group; we can impose structure on it. For example, we can model the variance as a power of the mean of each group (and thus only need to estimate one parameter, the exponent), or perhaps as the logarithm of one of the predictors in the model. All this is very easy with GLS (and
gls() in R).
Generalised least squares is IMHO a very underused statistical modelling technique. Instead of worrying about deviations from the model assumptions, model those deviations!