Homogeneity of variances but non normal residuals of a two way ANOVA

I have a problem, actually two. I am doing a bunch two-way ANOVAs and some of them can't fullfil the assumptions, no matter what I transform the data to.

1. Question: Does there exist a suitable non-paremetric alternative. I have looked at the Friedmann test several times, but as I understand it, my data needs to be paired to use it. Which my data is not.

Secondly some of my test only violate the normal distribution of residuals, and not the homogeneity of variances.

1. Question: Is there a way to correct your two-way ANOVA accordingly, and how would you do it in R? I have tried looking it up, but all I get is "What to do when you violate the homogeneity of variance assumption" which is what my data does not do!!

Data description: My data is continous (otherwise I wouldn't do ANOVA test). Both my factors have two levels. I have four treatments all composed of 12 observations except one which only has 11.

I hope someone find the time to help a hopefull individual as me :)

• Try taking the ranks of all your observations. Then doing the 2-way ANOVA on the ranks. Without seeing your data, I can't say whether it will work. May be OK, but also power may be impaired so you'd get no significant results. Ranks aren't normal, but maybe closer to normal than original observations. You'd still need to do diagnostics on residuals. // If you give more details of your design, maybe we can be more helpful: number of levels of each of the two factors, number of replications for each treatment combination would be a good stat. – BruceET Jun 25 at 15:03
• Both my factors have two levels, and there are 12 replicates per treatment except my fourth treatment which only has 11 replicates. Is that enough, I can't share my data from where I am right now. – Peder Nielsen Jun 25 at 15:08
• Can you include any such information in the original post? Please also explain more about your variables in the question; suitable analyses depend on having a better understanding of what your response variable is (among other things). – Glen_b Jun 26 at 1:47
• Are you still there? I have some fake data that might illustrate my Answer and would post that if you want to see it. If you're happy as is, that's fine. – BruceET Jun 26 at 3:31
• It would be awesome if you used some fake data to illustrate your point :) – Peder Nielsen Jun 26 at 7:33

Say, one factor has levels A,B, other has a,b: First, do a Kruskal-Wallis test on four groups Aa, Ab, Ba, Bb to see if anything at all is significant.

If so, try two 2-sample Wilcoxon (rank sum) tests on (i) A vs B, and then, separately, on (ii) a vs. b.

If you have significant results, then you can see whether the ANOVA on rank-transformed data mentioned in my Comment helps clarify your findings. You can hope they are not much different from the ANOVA you've already done.

[A full ANOVA design (original or rank-transformed) is slightly out of balance because of 12-12-12-11 replications you mentioned in your Comment. I assume you know, you can use a general linear model instead of basic ANOVA for a suitable approximate analysis. Neither K-W nor Wilcoxon RS test will have a problem with the slight imbalance.]

How to report findings for non-statisticians? You might give results of the first ANOVA you've already done. Note that data are non-normal, so that the exact P-values may be in doubt, but say to what extent you have been able to confirm that your basic findings (accept or reject) are OK-- by using nonparametric tests. Then in a footnote: give P-values of the Kruskal-Wallis and Wilcoxon tests. If you try to report your whole analytic procedure all at once in the methods section, you're likely to confuse everyone.

Addendum: Computations with simulated data: Simulate data for four groups in R: Aa, Ab, Ba, and Bb, according my notation above, where capital letters A, B are levels of Factor 1 and a, b are levels of Factor 2. Sample sizes are as in your experiment. The model for the simulated data is normal observations that are somewhat 'contaminated' by an exponential component in observation.

set.seed(1234) # for reproducibility
Aa = rnorm(12, 100, 10) + rexp(12, .5)
Ab = rnorm(12, 103, 10) + rexp(12, .5)
Ba = rnorm(12, 100, 10) + rexp(12, .5)
Bb = rnorm(11, 106, 10) + rexp(11, .5)
boxplot(Aa,Ab,Ba,Bb, notch=T, names=c("Aa","Ab","Ba","Bb"),
col="skyblue2", pch=19)


Nonparametric tests: In the boxplots below, the notches in the sides of the boxes indicate nonparametric confidence intervals, calibrated for pairwise comparisons. Non-overlapping intervals suggest differences in medians at the 5% level. Thus it seems that Aa & Ba are not significantly different; also Ab & Bb are are not significantly different. So this preliminary look at the data suggests that Factor 1 (with levels A and B) is not significant. By contrast, it seems that Factor 2 (levels a and b) may be significant.

A Kruskal-Wallis test of the four samples shown in the boxplot finds significant differences among them.

kruskal.test(list(Aa,Ab,Ba,Bb))

Kruskal-Wallis rank sum test

data:  list(Aa, Ab, Ba, Bb)
Kruskal-Wallis chi-squared = 13.694, df = 3, p-value = 0.003352


For subsequent analyses, it is convenient to 'stack' the data, putting all 47 observations into vector X and using factor variables g1 and g2 to represent the two main effects.

X = c(Aa, Ab, Ba, Bb)
g1 = as.factor(c(rep(1,24), rep(2,23)))
g2 = as.factor( c( rep(1,12), rep(2,12), rep(1,12), rep(2,11) ) )


In particular, we can now use two two-sample Wilcoxon (rank sum) tests to see that Factor 1 is not significant and Factor 2 is highly significant.

wilcox.test(X ~g1)$$p.val [1] 0.4927196 wilcox.test(X ~g2)$$p.val
[1] 0.0003489386


The nonparametric tests shown above do not assume that data are normal. They do pretty much settle the major research issues of the experiment.

ANOVA tests: Now we look at a traditional two-way ANOVA, which does assume normal data (most profitably checked by looking at residuals). The tests in the following ANOVA table (if valid) confirm what we found via nonparametric tests above: Factor 1 is nowhere near significant and Factor 2 is highly significant. Also, the interaction between the two factors is not significant.

aov.org = aov(X ~ g1*g2)    # ANOVA procedure on original data
summary(aov.out)
Df Sum Sq Mean Sq F value   Pr(>F)
g1           1      5     4.8   0.060 0.807620
g2           1   1154  1153.9  14.432 0.000452 ***
g1:g2        1    145   145.1   1.815 0.184929
Residuals   43   3438    80.0
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


While the residuals pass a Shapiro-Wilk test of normality, we see other indications later that the residuals may not be exactly normal. Even so, it seems to use the results of the ANOVA on the original data--especially because its major conclusions agree with nonparametric tests.

Another nonparametric procedure is to rank-transform the original data, replacing each of the 47 observations with its rank from 1 through 47, and then to do a standard ANOVA procedure on the ranks. Results, shown below, are essentially the same as for the ANOVA on the original data.

aov.rnk = aov(rank(X)~g1*g2)
summary(aov.rnk)
Df Sum Sq Mean Sq F value   Pr(>F)
g1           1     93    92.7   0.656 0.422269
g2           1   2243  2243.4  15.883 0.000256 ***
g1:g2        1    238   238.4   1.688 0.200805
Residuals   43   6073   141.2
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


Below we compare the residuals of the two ANOVAs just above (on original data, at the left, and on ranks). Normal probability plots of residuals seem to show that the residuals from the ranked data may be somewhat more nearly normal than those from the original data. However, both sets of residuals pass a Shapiro-Wilk test (P-values about 0.16 for original data and about 0.56 for ranks).

Of course, for your real data, a rank transformation, might make the difference between extremely non-normal residuals and residuals that are more nearly normal.

 par(mfrow=c(1,2))
aov.org = aov(X ~ g1*g2);  res.org = aov.org$$resid qqnorm(res.org, pch=20); qqline(res.out) aov.rnk = aov(rank(X)~g1*g2); res.rnk = aov.rnk$$resid
qqnorm(res.rnk, pch=20);  qqline(res.rnk)
par(mfrow=c(1,1))