Analysis of variance with not normally distributed residuals : how important is normality?

I am using gls and anova to analyse my data. I use gls to aply weights. I have one factor (tree genotype) and I analyse its influence on soil content.

Here is an exemple of my data with one variable (potassium) :

> my_data <- structure(list(Sample_ID = c("08C1", "08C2", "08C3", "09C1",
+                                         "09C2", "09C3", "10C1", "10C2", "10C3", "13C1", "13C2", "13C3",
+                                         "16C1", "16C2", "21C1", "21C2", "21C3", "23C1", "23C2", "23C3",
+                                         "25C1", "25C2", "25C3", "29C1", "29C2", "29C3", "33C1", "33C2",
+                                         "33C3"), Genotype = c("PB08", "PB08", "PB08", "PB09", "PB09",
+                                                               "PB09", "PB10", "PB10", "PB10", "PB13", "PB13", "PB13", "PB16",
+                                                               "PB16", "PB21", "PB21", "PB21", "PB23", "PB23", "PB23", "PB25",
+                                                               "PB25", "PB25", "PB29", "PB29", "PB29", "PB33", "PB33", "PB33"
+                                         ), K = c(0.85299241359094, 0.465133054980756, 0.00690267582675159,
+                                                  1.01673068757058, 1.030275834371, 0.834524714326693, 1.03316209040447,
+                                                  1.10298412805863, 1.18155841491109, 0.838705999451807, 0.883919392948746,
+                                                  0.901960535626635, 0.75594651348742, 0.756023068853717, 0.958311561443041,
+                                                  0.673810456005313, 1.08746537546638, 0.960314352517763, 1.0648664379158,
+                                                  1.00746055772745, 0.847776805354443, 1.03812953182887, 0.919747454011172,
+                                                  0.632561226408033, 0.664781889579858, 0.833083943531637, 0.92164628717724,
+                                                  0.622415179403636, 1.03783101178149)), row.names = c(NA, 29L), class = "data.frame")


If you look at homogeneity of variance for this data here's what you get :

model1 <- gls(K~Genotype, data=my_data)
layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE))
plot(resid(model1, type = "normalized") ~ fitted(model1, type = "normalized"), col = "blue", xlab="Fitted values", ylab="Standardized residuals")
boxplot(resid(model1, type = "normalized") ~ my_data$Genotype, col = "blue", xlab="Fitted values", ylab="Standardized residuals") hist((resid(model1, type = "normalized") - mean(resid(model1, type = "normalized"))) / sd(resid(model1, type = "normalized")), freq = FALSE, xlab="Standardized residuals", main=""); curve(dnorm, add = TRUE, col = "red") qqnorm((resid(model1, type = "normalized") - mean(resid(model1, type = "normalized"))) / sd(resid(model1, type = "normalized")), col="red", main=""); abline(0,1) bartlett.test(resid(model1, type = "normalized") ~ fitted(model1, type = "normalized")) shapiro.test(resid(model1, type = "normalized"))  > bartlett.test(resid(model1, type = "normalized") ~ fitted(model1, type = "normalized")) Bartlett test of homogeneity of variances data: resid(model1, type = "normalized") by fitted(model1, type = "normalized") Bartlett's K-squared = 27.957, df = 9, p-value = 0.0009697 > shapiro.test(resid(model1, type = "normalized")) Shapiro-Wilk normality test data: resid(model1, type = "normalized") W = 0.92237, p-value = 0.03503  And if you add weights in the gls you solve the problem for homogeneity but it makes it worse for normality : model2 <- gls(K~Genotype, weights=varIdent(form=~1|Genotype), data=my_data) layout(matrix(c(1,2,3,4), 2, 2, byrow = TRUE)) plot(resid(model2, type = "normalized") ~ fitted(model2, type = "normalized"), col = "blue", xlab="Fitted values", ylab="Standardized residuals") boxplot(resid(model2, type = "normalized") ~ my_data$Genotype, col = "blue", xlab="Fitted values", ylab="Standardized residuals")
hist((resid(model2, type = "normalized") - mean(resid(model2, type = "normalized"))) / sd(resid(model2, type = "normalized")), freq = FALSE, xlab="Standardized residuals", main=""); curve(dnorm, add = TRUE, col = "red")
qqnorm((resid(model2, type = "normalized") - mean(resid(model2, type = "normalized"))) / sd(resid(model2, type = "normalized")), col="red", main=""); abline(0,1)
bartlett.test(resid(model2, type = "normalized") ~ fitted(model2, type = "normalized"))
shapiro.test(resid(model2, type = "normalized"))


> bartlett.test(resid(model2, type = "normalized") ~ fitted(model2, type = "normalized"))

Bartlett test of homogeneity of variances

data:  resid(model2, type = "normalized") by fitted(model2, type = "normalized")
Bartlett's K-squared = 3.5453e-12, df = 9, p-value = 1

> shapiro.test(resid(model2, type = "normalized"))

Shapiro-Wilk normality test

data:  resid(model2, type = "normalized")
W = 0.89024, p-value = 0.005766


Also, if I compare AIC, it proves my model with weights has a much better fit :

> AIC(model1)
[1] 21.60778
> AIC(model2)
[1] 6.010106
> anova(model1,model2)
Model df       AIC      BIC    logLik   Test  L.Ratio p-value
model1     1 11 21.607779 31.99661  0.196111
model2     2 20  6.010106 24.89889 16.994947 1 vs 2 33.59767   1e-04


All this considered, I choose to do my anova with my model 2 as follows :

> anova(model2)
Denom. DF: 19
numDF   F-value p-value
(Intercept)     1 390068270  <.0001
Genotype        9        23  <.0001


Considering normality isn't met, are the results for my last anova trust-worty?

A fellow coworker pointed this out to me and told me I should use PERMANOVA when normality isn't met because this test does not require normality. So I decided to give it a shot :

> adonis(K~Genotype, data=my_data, method="bray", permutations=999)

Call:
adonis(formula = K ~ Genotype, data = my_data, permutations = 999,      method = "bray")

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
Genotype   9   0.45880 0.050977  1.2793 0.37732    0.1
Residuals 19   0.75712 0.039849         0.62268
Total     28   1.21592                  1.00000


So here it seems like I have oposite results from my two analysis. Also, because I'm doing an univariate analysis with adonis, homogeneity should be met, but as previously demonstrated, it is not the case.

I also tried putting a transformation in my gls (I tried sqrt, transformTukey and rank) but none of them gives normality.

What should I do?

EDIT

Here's the results for the betadisper:

> betadisper(d=vegdist(my_data$$K),group=my_data$$Genotype)

Homogeneity of multivariate dispersions

Call: betadisper(d = vegdist(my_data$$K), group = my_data$$Genotype)

No. of Positive Eigenvalues: 6
No. of Negative Eigenvalues: 0

Average distance to median:
PB08      PB09      PB10      PB13      PB16      PB21      PB23      PB25      PB29      PB33
4.070e-01 3.501e-02 2.236e-02 1.212e-02 2.532e-05 7.915e-02 1.722e-02 3.373e-02 4.573e-02 8.436e-02

Eigenvalues for PCoA axes:
PCoA1     PCoA2     PCoA3     PCoA4     PCoA5     PCoA6
9.296e-01 2.774e-01 8.741e-03 1.759e-04 2.239e-06 1.571e-08


I found that to interpret those results (with a p-value) you can either do an anova or a permutest :

> modk <- betadisper(d=vegdist(my_data$$K),group=my_data$$Genotype)
> permutest(modK, permutations = 999)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     9 0.38476 0.042751 2.0854    999  0.039 *
Residuals 19 0.38951 0.020501
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> anova(modK)
Analysis of Variance Table

Response: Distances
Df  Sum Sq  Mean Sq F value Pr(>F)
Groups     9 0.38476 0.042751  2.0854  0.085 .
Residuals 19 0.38951 0.020501
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1


It seems to me like the permutest tells me that I don't have homogeneity of variances and the anova sais the opposite...

SECOND EDIT

Since I'm not sure about homogeneity I tried transforming my data and it seemed to have solve the problem :

> my_data <- cbind(my_data, K_sqrt=sqrt(my_data$$K)) > bd <- betadisper(d=vegdist(my_data$$K_sqrt),group=my_data$Genotype) > permutest(bd, permutations = 999) Permutation test for homogeneity of multivariate dispersions Permutation: free Number of permutations: 999 Response: Distances Df Sum Sq Mean Sq F N.Perm Pr(>F) Groups 9 0.23187 0.025764 1.5013 999 0.118 Residuals 19 0.32606 0.017161 > anova(bd) Analysis of Variance Table Response: Distances Df Sum Sq Mean Sq F value Pr(>F) Groups 9 0.23187 0.025764 1.5013 0.2176 Residuals 19 0.32606 0.017161  Is this ok? Still I'd like to know wich one is better : anova or permutest? And what's the difference? • The vegan::adonis function is deprecated - you should be using vegan:adonis2. – André.B May 27 at 20:56 1 Answer Before changing weights the normality assumption was fine - ANOVAs are quite robust to deviations and that was a nice unimodal peak. However, equality of variance (an arguably more important assumption) was not met. Ignore Shapiro-Wilk in the first instance. It is much more sensitive to deviations from normality than the actual test you are using is. The gls weighted analysis corrected the heteroscedasticity but violated the normality assumption. This is not a huge problem if you have a decent sample size (i.e. n > 30), as the central limit theorem applies, but it does mean the model is unreliable for predictions. In saying that, it looks like you have quite a small sample size. Akaike's information criterion (AIC) is completely inappropriate in this case. A) your sample size is much too small. If anything you should be using AICc (AIC's second order derivative), and; B) AIC (or AICc) cannot be used to compare two completely different methods that are maximising two different likelihoods. PERMANOVA is probably the most sensible approach if you want to see how tree genotype is influencing multiple soil parameters, as it avoids some type II error issues. However, as mentioned in the comments, vegan::adonis is deprecated, you should be using vegan::adonis2. Moreover, PERMANOVA does require that the multivariate spread of your groups is approximately equal. This is similar to the homogeneity of variance assumption for an ANOVA. You can check this with the vegan::betadisper function. • I used betadisper to calculate multivariante dispersions : mod<-betadisper(d=vegdist(my_data$K),group=my_data\$Genotype) now to check if the homogeneity of variance assumption, should I do anova(mod) or permutest(mod,pairwise=TRUE,permutations=99) ? – Karelle Rheault May 28 at 14:45
• Multivariate dispersions is the multivariate equivalent of the ANOVA's homogeneity assumption so you don't need to do anything more. Hopefully, betadisper returned favourable results... – André.B May 28 at 21:40
• In this case I don't understand how to interpret the results of betadisper... I'm use to look at a p-value but here I have Average distance to median and Eigenvalues for PCoA axes... – Karelle Rheault May 29 at 14:24
• Could you update your question with the values? I am happy to help you interpret them. – André.B May 29 at 21:07
• There are methods for permutational ANOVAs, but PERMANOVA is specifically multivariate. If you do the univariate version then you would just check for heteroscedasticity as you normally would (i.e. residuals v.s. fitted values plot - *Bartlett's and Levene's tests are much to sensitive). If you do multiple univariate models though you need to correct for multiple comparisons - PERMANOVA avoids this - you can do this with the predictmeans::predictmeans() function in R. Here is a blog about univariate permutational tests: statmethods.wordpress.com/2012/05/21/permutation-tests-in-r – André.B May 30 at 21:10