18

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
whuber
  • 281,159
  • 54
  • 637
  • 1,101
Sabine
  • 181
  • 1
  • 1
  • 4
  • 2
    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: http://www.ats.ucla.edu/stat/R/seminars/Repeated_Measures/repeated_measures.htm and http://www.personality-project.org/R/r.anova.html – Roman Luštrik May 16 '12 at 13:30
  • I'm guessing your dependent variable is "time", can you say more about that? Is this time to event (a.k.a. [survival analysis](http://en.wikipedia.org/wiki/Survival_analysis))? Can you explain what your experiment is & what you are trying to understand (in layman's terms)? – gung - Reinstate Monica May 16 '12 at 22:15
  • Just a general comment about the levels for treatment. Would it be appropriate to have separate independent variables for whatever CY and HL represent? – RobJackson28 May 17 '12 at 02:36
  • I treated the proteins with the enzymes and then conducted an experiment where all samples started out at the same value and the measured variable decreased with time. The variable "time" in the example is the time it took the samples to decrease to half the original value. I had been thinking about treating CL and HY as different factors as well, but then I would also have to decide between CL+HX mix and CL-HY sim. I figured they are all different treatments and I could look at them as that without making factors out of it. Whereas, if it is easy to do with R, why not? I'll think about it. – Sabine May 17 '12 at 06:59
  • 4
    A very good alternative to rank tests is to use permutation ANOVA's (http://www.uvm.edu/~dhowell/StatPages/More_Stuff/Permutation%20Anova/PermTestsAnova.html). 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... – Mikko May 16 '12 at 14:01
  • 1
    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
  • 1
    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
  • Thank you very much for your replies. I can follow that so far and made some residual plots on the data myself by now and this really doesn't look too bad in terms of normality. But are you implying, that equality of variance would not matter that much either? – Sabine May 17 '12 at 06:55
  • You can also try some robust methods such as trimmed mean/median base methods with bootstrap-t or percentile-bootstrap modified Patel-Hoel method. –  May 16 '12 at 18:40
  • 1
    @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`) – Mikko May 18 '12 at 11:30
  • @Largh. I agree but I think nonnormality can be a problem if there is high skewness or a heavytailed distribution for the noise component. Least squares is very sensitive to outliers and heavytailed distributions. In addition to the constant variable assumption people often neglect the assumption that the X variables are observed without noise. If there is a significant noise component to the Xs an error in variables regression model might give a far different solution than least squares. – Michael R. Chernick May 18 '12 at 11:42

2 Answers2

12

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.

Community
  • 1
gung - Reinstate Monica
  • 132,789
  • 81
  • 357
  • 650
  • Re my 2nd point, I recently wrote an answer [here](http://stats.stackexchange.com/questions/12262/what-if-residuals-are-normally-distributed-but-y-is-not/33320#33320) that may help make this issue clearer. You may want to read it. – gung - Reinstate Monica Aug 06 '12 at 13:19
9

(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)
Community
  • 1
David LeBauer
  • 7,060
  • 6
  • 44
  • 89
  • 1
    +1, this is a nice list of some available options to explore. – gung - Reinstate Monica May 19 '12 at 03:11
  • Below are comments originally exchanged on Stack Overflow, somehow lost during migration. See http://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](http://stats.stackexchange.com/users/442/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](http://stats.stackexchange.com/users/1381/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](http://stats.stackexchange.com/users/442/henrik) May 16 at 17:03 – chl May 19 '12 at 20:35
  • @Henrik nope. I only read some of the the Anderson 2005 documentation of the [implementation of PERMANOVA in FORTRAN](http://www.stat.auckland.ac.nz/~mja/prog/PERMANOVA_UserNotes.pdf), and then found through searching that it was implemented in the R package. – [David](http://stats.stackexchange.com/users/1381/david) May 16 at 18:48 – chl May 19 '12 at 20:35