Questions about model analysis, residuals for normality, and appropriate statistical test usage

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")

• The QQ plot is fine. Why are you using a whole bunch of tests to test normality? – Glen_b -Reinstate Monica Sep 16 '17 at 8:20
• Thank you. Is there any objective way to know when too many points are off the QQ Line that skewness or heavy tails? Based on the plots above, I was thinking there was an excessive amount of this. I applied the battery of normality tests to serve as a possible indicator regarding if I should be concerned or not (if the tests suggested that it's possible the residuals are normally distributed, I would have been less worried about the QQ Plot) and just to play around with functions in R. – rls Sep 16 '17 at 15:07
• So many things to say -- but these are all discussed on site 1. Real data aren't normal, it's a model; there's little point in testing (e.g. see Harvey's answer) - a test of normality answers the wrong question. 2. You can see how much variation you'd expect to see if the data were normal by using the approach here. 3. However, that still doesn't really answer the question we're interested in -- ... ctd – Glen_b -Reinstate Monica Sep 16 '17 at 22:30
• ctd... which is not "are the data drawn from a normal distribution?" (we know they aren't) but instead "how much difference does it make?" (i.e. is it close enough that I would happily use this approach anyway?). That's a rather subtler question and depends on a host of things - to what extent and in what manner you have non-normality; the sensitivity of the specific analyses you're doing to that; how much you care about the size of the impact (in some situations a given effect on significance level, power, CI-coverage or whatever may not matter; other times it might matter a lot) – Glen_b -Reinstate Monica Sep 16 '17 at 22:35
• In short, the fact that a goodness of fit test rejects or doesn't reject is hardly relevant to a question of how much it matters for your analysis. A Q-Q plot can tell you a little about how much it might matter, but you need to get used to reading them and not over-react to small deviations in the tail -- especially if it's showing a slightly shorter tail. Simulation of data sets with similar characteristics to your own can give you some sense of how much difference those kinds of deviation might make to your inference. – Glen_b -Reinstate Monica Sep 16 '17 at 22:40