I've conducted an experiment on whether the presence of an insect that uses animal carcasses for food and reproduction impacts the soil in the immediate area and am now attempting to analyze the data.
During my analysis, I've started to question the validity of my statistical approach/test, primarily because I am not able to achieve normal residuals, but also because of my lack of experience in statistical analysis (I have only taken one graduate course in regression).
My data is fairly simple: three treatment levels (soil alone, soil with carcass, soil with insect and carcass) with two measurements (organic matter and soil ph). I have additional data measurements, but I'm starting the analysis with this set and hoping to learn from it.
My model is straightforward:
model <- organic_matter ~ trt + soil_ph + trt:soil_ph
aovOrganicMatter <- aov(model, data=df, na.action=na.omit)
The output looks about as one might expect:
> summary.aov(aovOrganicMatter)
Df Sum Sq Mean Sq F value Pr(>F)
trt 2 49.16 24.580 6.925 0.00292 **
soil_ph 1 19.27 19.266 5.428 0.02571 *
trt:soil_ph 2 30.35 15.177 4.276 0.02180 *
Residuals 35 124.22 3.549
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
3 observations deleted due to missingness
Similarly, the summary seems straightforward:
> summary.lm(aovOrganicMatter)
Call:
aov(formula = model, data = df, na.action = na.omit)
Residuals:
Min 1Q Median 3Q Max
-4.3060 -1.2696 0.3315 1.1608 3.7315
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 108.056 4.952 21.820 < 2e-16 ***
trtInsect -22.997 9.800 -2.347 0.024730 *
trtNoInsect -11.040 22.029 -0.501 0.619394
soil_ph -4.770 1.286 -3.710 0.000715 ***
trtInsect:soil_ph 5.228 1.793 2.916 0.006153 **
trtNoInsect:soil_ph 3.623 4.268 0.849 0.401714
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.884 on 35 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.4429, Adjusted R-squared: 0.3634
F-statistic: 5.566 on 5 and 35 DF, p-value: 0.000711
The resulting plots start to give the impression that the residuals might not be normal, though:
In particular, the QQ Plot seems to show that the tails might be a bit heavy and that there is some skewness. I've run a number of tests on the residuals (Anderson-Darling, Shapiro-Wilk, Anscombe-Glynn, and Jarque-Bera) that also suggest I can't conclude that the distribution is normal.
I've looked over most posts on Cross Validated to learn more about normal QQ Plots, but I still don't think my eye/inuition is that strong. I've also tried transforming the data (using square root and log2), but the resulting residual QQ Plots never look much better to me, and the tests above do not change in their results (that is, they still suggest non-normality).
So, I appeal to anyone out there who has more experience interpreting QQ Plots and designing models than myself. Does my model seem appropriate for the analysis? Does the QQ Plot of residuals look "normal enough," and can I therefore continue to use the Anova? If it is not normal, how would you advise transforming the data, or should I change my method of analysis entirely? Lastly, is there a more suitable version of Analysis of Variance in R, for this instance, than aov?
Here is a dput of the data (note that I've included the full set of data here, but the above analysis removes three outliers. I'm also questioning this: one (the 72.5 value) I'd agree is an outlier, the other two seem reasonable, though admittedly high):
structure(list(id = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L,
11L, 12L, 13L, 14L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L,
25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L,
38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L), sample = c("nbom1",
"nbom3", "nbom4", "nbom5", "nbom6",
"nbom7", "nbom8", "nbom9", "nbom10",
"nbom11", "nbom12", "nbom13",
"nbom14", "m199", "m175", "m266", "m155", "m156",
"m164", "m166", "m173", "m174", "m176", "m185", "m186", "m187",
"m188", "m197", "m198", "m200", "m224", "m225", "m227", "m360",
"nb01", "nb02", "nb03", "nb04", "nb05",
"nb06", "nb07", "nb08", "nb09", "nb10"
), trt = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L), .Label = c("Soil", "Insect", "NoInsect"), class = "factor"),
organic_matter = c(89.4, 92.2, 90.6, 90.4, 91.8, 90.5, 92.9,
89.9, 91.8, 85.8, 88.2, 87.7, 86, 85.7, 89.5, 87.3, 89.5,
88.2, 88.7, 85.5, 87.4, 90.7, 92, 90.3, 89.9, 88.1, 88.6,
84.1, 86, 85.7, 85.5, 89.2, 90.5, 88.9, 89.2, 92.8, 90.5,
89.1, 94.2, 95.5, 72.5, 99.9, 89.6, 91.4), soil_ph = c(3.8,
3.5, 3.5, 3.8, 3.4, 3.4, 3.4, 3.6, 3.8, 4.4, 4.3, 4.5, 4.4,
6.7, 6.9, 6.7, 6.8, 6, 7.2, 6.3, 7, 6.7, 7, 7, 6.8, 6.9,
7, 7.3, 7, 6.4, 6.7, 6.5, 6.9, 6.1, 5.2, 5.1, 5.5, 5.5, 5.3,
5.3, 5.2, 5.1, 5, 5.3)), .Names = c("id", "sample", "trt",
"organic_matter", "soil_ph"), row.names = c(1L, 2L, 3L, 4L, 5L,
6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 16L, 17L, 18L, 19L,
20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L,
33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 45L
), class = "data.frame")