I have data from 3 groups of algae biomass ($A$, $B$, $C$) which contain unequal sample sizes ($n_A=15$, $n_B=13$, $n_C=12$) and I would like compare if these groups are from the same population.
One-way ANOVA would definitely be the way to go, however upon conducting normality tests on my data, heteroskedascity seems to the main issue. My raw data, without any transformation, produced a ratio of variances ($F_{\max} = 19.1$) which is very much higher than the critical value ($F_{\rm crit} = 4.16$) and therefore I cannot perform one-way ANOVA.
I also tried transformation to normalise my data. Even after trials of various transformations (log, square root, square), the lowest $F_{\max}$ produced after transformation with a $\log_{10}$ transformation was $7.16$, which was still higher compared to $F_{\rm crit}$.
Can anyone here advise me on where to go from here? I can't think of other methods of transformation to normalise by data. Are there any alternatives to a one-way ANOVA?
P.S.: my raw data are below:
A: 0.178 0.195 0.225 0.294 0.315 0.341 0.36 0.363 0.371 0.398 0.407 0.409 0.432
0.494 0.719
B: 0.11 0.111 0.204 0.416 0.417 0.441 0.492 0.965 1.113 1.19 1.233 1.505 1.897
C: 0.106 0.114 0.143 0.435 0.448 0.51 0.576 0.588 0.608 0.64 0.658 0.788 0.958
Best Answer
There are a number of options available when dealing with heteroscedastic data. Unfortunately, none of them is guaranteed to always work. Here are some options I'm familiar with:
Update: Here is a demonstration in
Rof some ways of fitting a linear model (i.e., an ANOVA or a regression) when you have heteroscedasticity / heterogeneity of variance.Let's start by taking a look at your data. For convenience, I have them loaded into two data frames called
my.data(which is structured like above with one column per group) andstacked.data(which has two columns:valueswith the numbers andindwith the group indicator).We can formally test for heteroscedasticity with Levene's test:
Sure enough, you have heteroscedasticity. We'll check to see what the variances of the groups are. A rule of thumb is that linear models are fairly robust to heterogeneity of variance so long as the maximum variance is no more than $4\!\times$ greater than the minimum variance, so we'll find that ratio as well:
Your variances differ substantially, with the largest,
B, being $19\!\times$ the smallest,A. This is a problematic level of heteroscedsaticity.You had thought to use transformations such as the log or square root to stabilize the variance. That will work in some cases, but Box-Cox type transformations stabilize variance by squeezing the data asymmetrically, either squeezing them downwards with the highest data squeezed the most, or squeezing them upwards with the lowest data squeezed the most. Thus, you need the variance of your data to change with the mean for this to work optimally. Your data have a huge difference in variance, but a relatively small difference amongst the means and medians, i.e., the distributions mostly overlap. As a teaching exercise, we can create some
parallel.universe.databy adding $2.7$ to allBvalues and $.7$ toC's to show how it would work:Using the square root transformation stabilizes those data quite well. You can see the improvement for the parallel universe data here:
Rather than just trying different transformations, a more systematic approach is to optimize the Box-Cox parameter $\lambda$ (although it is usually recommended to round that to the nearest interpretable transformation). In your case either the square root, $\lambda = .5$, or the log, $\lambda = 0$, are acceptable, though neither actually works. For the parallel universe data, the square root is best:
Since this case is an ANOVA (i.e., no continuous variables), one way to deal with heterogeneity is to use the Welch correction to the denominator degrees of freedom in the $F$-test (n.b.,
df = 19.445, a fractional value, rather thandf = 38):A more general approach is to use weighted least squares. Since some groups (
B) spread out more, the data in those groups provide less information about the location of the mean than the data in other groups. We can let the model incorporate this by providing a weight with each data point. A common system is to use the reciprocal of the group variance as the weight:This yields slightly different $F$ and $p$-values than the unweighted ANOVA (
4.5089,0.01749), but it has addressed the heterogeneity well:Weighted least squares is not a panacea, however. One uncomfortable fact is that it is only just right if the weights are just right, meaning, among other things, that they are known a-priori. It does not address non-normality (such as skew) or outliers, either. Using weights estimated from your data will often work fine, though, particularly if you have enough data to estimate the variance with reasonable precision (this is analogous to the idea of using a $z$-table instead of a $t$-table when you have $50$ or $100$ degrees of freedom), your data are sufficiently normal, and you don't appear to have any outliers. Unfortunately, you have relatively few data (13 or 15 per group), some skew and possibly some outliers. I'm not sure that these are bad enough to make a big deal out of, but you could mix weighted least squares with robust methods. Instead of using the variance as your measure of spread (which is sensitive to outliers, especially with low $N$), you could use the reciprocal of the inter-quartile range (which is unaffected by up to 50% outliers in each group). These weights could then be combined with robust regression using a different loss function like Tukey's bisquare:
The weights here aren't as extreme. The predicted group means differ slightly (
A: WLS0.36673, robust0.35722;B: WLS0.77646, robust0.70433;C: WLS0.50554, robust0.51845), with the means ofBandCbeing less pulled by extreme values.In econometrics the Huber-White ("sandwich") standard error is very popular. Like the Welch correction, this does not require you to know the variances a-priori and doesn't require you to estimate weights from your data and/or contingent on a model that may not be correct. On the other hand, I don't know how to incorporate this with an ANOVA, meaning that you only get them for the tests of individual dummy codes, which strikes me as less helpful in this case, but I'll demonstrate them anyway:
The function
vcovHCcalculates a heteroscedasticicy consistent variance-covariance matrix for your betas (your dummy codes), which is what the letters in the function call stand for. To get standard errors, you extract the main diagonal and take the square roots. To get $t$-tests for your betas, you divide your coefficient estimates by the SEs and compare the results to the appropriate $t$-distribution (namely, the $t$-distribution with your residual degrees of freedom).For
Rusers specifically, @TomWenseleers notes in the comments below that the ?Anova function in thecarpackage can accept awhite.adjustargument to get a $p$-value for the factor using heteroscedasticity consistent errors.You can try to get an empirical estimate of what the actual sampling distribution of your test statistic looks like by bootstrapping. First, you create a true null by making all group means exactly equal. Then you resample with replacement and calculate your test statistic ($F$) on each bootsample to get an empirical estimate of the sampling distribution of $F$ under the null with your data whatever their status with regard to normality or homogeneity. The proportion of that sampling distribution that is as extreme or more extreme than your observed test statistic is the $p$-value:
In some ways, bootstrapping is the ultimate reduced assumption approach to conducting an analysis of the parameters (e.g., means), but it does assume that your data are a good representation of the population, meaning you have a reasonable sample size. Since your $n$'s are small, it may be less trustworthy. Probably the ultimate protection against non-normality and heterogeneity is to use a non-parametric test. The basic non-parametric version of an ANOVA is the Kruskal-Wallis test:
Although the Kruskal-Wallis test is definitely the best protection against type I errors, it can only be used with a single categorical variable (i.e., no continuous predictors or factorial designs) and it has the least power of all strategies discussed. Another non-parametric approach is to use ordinal logistic regression. This seems odd to a lot of people, but you only need to assume that your response data contain legitimate ordinal information, which they surely do or else every other strategy above is invalid as well:
It may not be clear from the output, but the test of the model as a whole, which in this case is the test of your groups, is the
chi2underDiscrimination Indexes. Two versions are listed, a likelihood ratio test and a score test. The likelihood ratio test is typically considered the best. It yields a $p$-value of0.0363.