# How to run two-way ANOVA on data with neither normality nor equality of variance in R?

I am working on my master thesis at the moment and planned on running the statistics with SigmaPlot. However, after spending some time with my data I came to the conclusion that SigmaPlot might not be fit for my problem (I may be mistaken) so I started my first attempts in R, which did not exactly make it easier.

The plan was to run a simple TWO-WAY-ANOVA on my data which results from 3 different proteins and 8 different treatments on those, so my two factors are proteins and treatments. I tested for normality using both

> shapiro.test(time)


and

> ks.test(time, "norm", mean=mean(time), sd=sqrt(var(time)))


In both cases (maybe not surprising) I ended up with a non-normal distribution.

Which left me with the first questions of which test to use for equality of variances. I came up with

> chisq.test(time)


and the result was, that I don't have equality of variance in my data either.

I tried different data transformations (log, center, standardization), all of which did not solve my problems with the variances.

Now I am at a loss, how to conduct the ANOVA for testing which proteins and which treatments differ significantly from each other. I found something about a Kruskal-Walis-Test, but only for one factor (?). I also found things about ranking or randamization, but not yet how to implement those techniques in R.

Does anyone have a suggestion what I should do?

Edit: thank you for your answers, I am a little overwhelmed by the reading (it just seems getting more and more instead of less), but I will of course keep going.

Here an example of my data, as suggested (I am very sorry for the format, I couldn't figure out another solution or place to put a file. I am still new to this all.):

protein treatment   time
A   con 2329.0
A   HY  1072.0
A   CL1 4435.0
A   CL2 2971.0
A   CL1-HY sim  823.5
A   CL2-HY sim  491.5
A   CL1+HY mix  2510.5
A   CL2+HY mix  2484.5
A   con 2454.0
A   HY  1180.5
A   CL1 3249.7
A   CL2 2106.7
A   CL1-HY sim  993.0
A   CL2-HY sim  817.5
A   CL1+HY mix  1981.0
A   CL2+HY mix  2687.5
B   con 1482.0
B   HY  2084.7
B   CL1 1498.0
B   CL2 1258.5
B   CL1-HY sim  1795.7
B   CL2-HY sim  1804.5
B   CL1+HY mix  1633.0
B   CL2+HY mix  1416.3
B   con 1339.0
B   HY  2119.0
B   CL1 1093.3
B   CL2 1026.5
B   CL1-HY sim  2315.5
B   CL2-HY sim  2048.5
B   CL1+HY mix  1465.0
B   CL2+HY mix  2334.5
C   con 1614.8
C   HY  1525.5
C   CL1 426.3
C   CL2 1192.0
C   CL1-HY sim  1546.0
C   CL2-HY sim  874.5
C   CL1+HY mix  1386.0
C   CL2+HY mix  364.5
C   con 1907.5
C   HY  1152.5
C   CL1 639.7
C   CL2 1306.5
C   CL1-HY sim  1515.0
C   CL2-HY sim  1251.0
C   CL1+HY mix  1350.5
C   CL2+HY mix  1230.5

-
It would seem you are well on your way with R, but I feel your difficulties may not be R related at all. Perhaps you need to consult a statistician first to get the grips of the problem, only then you would tackle this practically. Your current question is pretty elaborate and may deter people from taking a whack at it. For the technical part, here's a few sites: ats.ucla.edu/stat/R/seminars/Repeated_Measures/… and personality-project.org/R/r.anova.html –  Roman Luštrik May 16 '12 at 13:30
A very good alternative to rank tests is to use permutation ANOVA's (uvm.edu/~dhowell/StatPages/More_Stuff/Permutation%20Anova/…). With that approach heteroscedasticity does not matter. The main reason people use rank tests is because they are computationally much easier. That does not matter any longer. We have R and efficient computers... –  Largh May 16 '12 at 14:01
What is reply? Note that we don't actually care if the response is normally distribution - we want to check if the residuals from the analysis are approximately normal with similar variances. But if there is actually a treatment effect then we wouldn't expect that the response would be normally distributed marginally. –  Dason May 16 '12 at 22:50
I'll note that I just grabbed the data you posted and ran the two way anova and did a shapiro wilks test on the residuals and this gave a p-value of 0.5022 which doesn't imply too much of a concern. –  Dason May 16 '12 at 23:01
@Sabine normality is often not a big concern what it comes to ANOVA, but your samples should come from the same population meaning that equality of variance is the most important, assumption after randomization, of course. You should do something, if your variances are not approximately equal (try ?bartlett.test) –  Largh May 18 '12 at 11:30

This may be more of a comment than an answer, but it won't fit as a comment. We may be able to help you here, but this may take a few iterations; we need more information.

First, what is your response variable?

Second, note that the marginal distribution of your response does not have to be normal, rather the distribution conditional on the model (i.e., the residuals) should be--it is not clear that you have examined your residuals. Furthermore, normality is the least important assumption of a linear model (e.g., an ANOVA); the residuals may not need to be perfectly normal. Tests of normality are not generally worthwhile (see here for a discussion on CV), plots are much better. I would try a qq-plot of your residuals. In R this is done with qqnorm(), or try qqPlot() in the car package. It's also worth considering the manner in which the residuals are non-normal: skewness is more damaging than excess kurtosis, in particular if the skews alternate directions amongst the groups.

If there really is a problem worth worrying about, a transformation is a good strategy. Taking the log of your raw data is one option, but not the only one. Note that centering and standardizing aren't really transformations in this sense. You want to look into the Box & Cox family of power transformations. And remember, the result doesn't have to be perfectly normal, just good enough.

Next, I don't follow your use of the chi-squared test for homogeneity of variance, although it may be perfectly fine. I would suggest you use Levene's test (use leveneTest() in car). Heterogeneity is more damaging than non-normality, but the ANOVA is pretty robust if the heterogeneity is minor. A standard rule of thumb is that the largest group variance can be up to four times the smallest without posing strong problems. A good transformation should also address heterogeneity.

If these strategies are insufficient, I would probably explore robust regression before trying a non-parametric approach.

If you can edit your question and say more about your data, I may be able to update this to provide more specific information.

-
Re my 2nd point, I recently wrote an answer here that may help make this issue clearer. You may want to read it. –  gung Aug 6 '12 at 13:19

(note: this answer was posted before the question was migrated and merged from SO, so details have been added to the question that are not addressed here. Many are addressed in the comments and the answer by @gung).

There are many different approaches, and this question has been covered elsewhere on this site. Here is a list of some approaches, with links to other questions on the site and some references.:

1. The Box-Cox power transformation can normalize residuals that are on a non-linear scale
2. ANOVA on ranked data is very easy but has reduced power and is difficult to interpret. See Conover and Iman, (1981)
3. Proportional Odds ordinal logistic model
4. Permutation Tests (Anderson and ter Braak 2003), implemented in and described by Anderson and as the adonis function in the R Vegan package
5. Bootstrapping
6. Hierarchical Bayesian modeling (Gelman 2005)
-
+1, this is a nice list of some available options to explore. –  gung May 19 '12 at 3:11
Below are comments originally exchanged on Stack Overflow, somehow lost during migration. See meta.stats.stackexchange.com/q/1157/930. –  chl May 19 '12 at 20:34
Do you know if there is an implementation of the ideas of Anderson and Braak? - Henrik May 16 at 15:15 –  chl May 19 '12 at 20:35
@Henrik there is an implementation in FORTRAN Anderson 2005 which is available through the R function adonis in the Vegan R package / David May 16 at 16:20 –  chl May 19 '12 at 20:35
Thanks. adonis seems to work even with univariate dvs. However, I have the feeling that it uses something like type 1 sums of squares as I get the message Terms added sequentially (first to last) whenever running it. Have you used it or can say anything about it? – Henrik May 16 at 17:03 –  chl May 19 '12 at 20:35
show 1 more comment