3
$\begingroup$

It is possible to find discrepancies between the results of ANOVA and post-hoc tests. The results of post-hoc tests are valid (with the exception of the Fisher LSD test) even if the overall P value for the ANOVA is greater than 0.05 (1, 2, 3, 4).

How should such results be presented and interpreted? There are some questions on this subject, but the answers are more or less theoretical (5, 6, 3, 4). I could not find any papers (biology, ecology and related fields) with such results.

Therefore, I would like to know how to present and interpret such results using a real example: I used a linear mixed model (function lmer, R package lme4) to analyse the effect of cover of two invasive plant species (Solidago spp.) on the diversity of invaded plant communities. The response variable is species richness, the explanatory variables are categorical - invasive species (2 levels: S. canadensis and S. gigantea), cover of invasive species (3 levels: C = ~0%, T = ~30% and D = ~80%) and their interaction.

Data

While the results of the anova (function Anova, type 3 and 2; package car) suggest no significant interaction, visual inspection of the boxplot suggests the opposite, which was confirmed by post-hoc tests (function avg_comparisons from the package marginaleffects). See below.

So how should I present and interpret these results?

(By the way, why are there no values for intercept in the type 2 Anova outputs?)

Boxplot

R script:

> library(lme4)
> library(car)
> library(MuMIn)
> library(marginaleffects)
> m1<-lmer(Hillall_0 ~ solidago * plot + (1 | locality), data=dat)

> Anova(m1, type="III", icontrasts=c("contr.sum", "contr.poly"))
Analysis of Deviance Table (Type III Wald chisquare tests)
Response: Hillall_0
Chisq Df Pr(>Chisq)
(Intercept)   185.7266  1  < 2.2e-16 ***
solidago        8.4379  1  0.0036749 **
plot           16.4066  2  0.0002737 ***
solidago:plot   4.7866  2  0.0913297 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> Anova(m1, type="II")
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: Hillall_0
       Chisq Df Pr(>Chisq) 
solidago      14.3900  1  0.0001486 ***
plot          12.8453  2  0.0016243 ** 
solidago:plot  4.7866  2  0.0913297 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> r.squaredGLMM(m1, pj2014 = F)
       R2m       R2c
[1,] 0.3518038 0.3518038

Post-hoc tests:

> avg_comparisons(m1, variables = "solidago", by = "plot")
Term                          Contrast plot Estimate Std. Error      z Pr(>|z|)   S 2.5 % 97.5 %
solidago mean(gigantea) - mean(canadensis)    C    -12.6       4.34 -2.905  0.00367 8.1 -21.1   -4.1
solidago mean(gigantea) - mean(canadensis)    D     -1.8       4.34 -0.415  0.67816 0.6 -10.3    6.7
solidago mean(gigantea) - mean(canadensis)    T    -14.1       4.34 -3.251  0.00115 9.8 -22.6   -5.6
Columns: term, contrast, plot, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type:  response

> avg_comparisons(m1, variables = "plot", by = "solidago")
Term          Contrast   solidago Estimate Std. Error      z Pr(>|z|)    S 2.5 % 97.5 %
plot mean(D) - mean(C) canadensis    -15.6       4.34 -3.596   <0.001 11.6 -24.1   -7.1
plot mean(D) - mean(C) gigantea       -4.8       4.34 -1.107    0.268  1.9 -13.3    3.7
plot mean(T) - mean(C) canadensis     -0.8       4.34 -0.184    0.854  0.2  -9.3    7.7
plot mean(T) - mean(C) gigantea       -2.3       4.34 -0.530    0.596  0.7 -10.8    6.2
Columns: term, contrast, solidago, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted
Type:  response

> avg_comparisons(m1, variables = "solidago", by = "plot", hypothesis = "pairwise")
Term Estimate Std. Error      z Pr(>|z|)   S   2.5 % 97.5 %
C - D    -10.8       6.13 -1.761   0.0783 3.7 -22.823   1.22
C - T      1.5       6.13  0.245   0.8068 0.3 -10.523  13.52
D - T     12.3       6.13  2.005   0.0450 4.5   0.277  24.32
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type:  response

> avg_comparisons(m1, variables = "plot", by = "solidago", hypothesis = "pairwise")
Term Estimate Std. Error      z Pr(>|z|)    S 2.5 % 97.5 %
(mean(D) - mean(C), canadensis) - (mean(D) - mean(C), gigantea)      -10.8       6.13 -1.761   0.0783  3.7 -22.8   1.22
(mean(D) - mean(C), canadensis) - (mean(T) - mean(C), canadensis)    -14.8       4.34 -3.412   <0.001 10.6 -23.3  -6.30
(mean(D) - mean(C), canadensis) - (mean(T) - mean(C), gigantea)      -13.3       6.13 -2.168   0.0301  5.1 -25.3  -1.28
(mean(D) - mean(C), gigantea) - (mean(T) - mean(C), canadensis)       -4.0       6.13 -0.652   0.5144  1.0 -16.0   8.02
(mean(D) - mean(C), gigantea) - (mean(T) - mean(C), gigantea)         -2.5       4.34 -0.576   0.5644  0.8 -11.0   6.00
(mean(T) - mean(C), canadensis) - (mean(T) - mean(C), gigantea)        1.5       6.13  0.245   0.8068  0.3 -10.5  13.52
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high
Type:  response
$\endgroup$

1 Answer 1

3
+50
$\begingroup$

There are two remarks with the comparison using the boxplot

Medians vs Means

Your boxplots are plotting the medians and it looks like the case of D cover has an opposite effect.

But this is not the case when you compare the means. Then you have effect sizes in the same direction.

So the plot with the medians make the differences in effect sizes appear a bit stronger.

Contrast                                 Estimate  Pr(>|z|)   
mean(gigantea) - mean(canadensis)    C  -12.6      0.00367 
mean(gigantea) - mean(canadensis)    D  -1.8       0.67816 
mean(gigantea) - mean(canadensis)    T  -14.1      0.00115 

Multiple comparisons

Still you have one effect size that is a lot different, why wouldn't this result into a significant interaction? The reason is that you have three levels which increases the odds that one of them happens to be different. Your case is a bit borderline. The result is not very significant, but neither is it insignificant. Remember that absence of statistical significance is not proof that there is no effect.

If you would compute a t-statistic out of the pairwise differences between those three effects. For example, with one of those comparisons you get an effect size difference of $$(-12.6) - (-1.8) = -10.8$$ then you get the comparisons in the lines from the final table

Term                                                                 Estimate  z        Pr(>|z|)
(mean(D) - mean(C), canadensis) - (mean(D) - mean(C), gigantea)     -10.8     -1.761    0.0783
(mean(D) - mean(C), canadensis) - (mean(T) - mean(C), gigantea)     -13.3     -2.168    0.0301 
(mean(T) - mean(C), canadensis) - (mean(T) - mean(C), gigantea)       1.5      0.245    0.8068

The z/t-scores of 0.245, -2.168 and -1.761 are not very large, or at the border of typically applied confidence levels. Compare with this figure from The theory behind Tukey's HSD test.

comparison

The individual comparisons can be significant at a 5% level. Like the individual z/t-score of -2.168 is described with a p-value of 0.0301, and this is what you intuitively see in the boxplots. However, when you look at multiple comparisons and wish to compensate for that, then you need to be stricter and use the blue or the red borders from the image above.

$\endgroup$
4
  • 1
    $\begingroup$ Halfway the post, I have been writing z/t-score instead of plain z-score or t-score. The table from the avg_comparison function confuses me. I believe that it should be a t-statistic that is computed and not a z-score. Having little motivation to figure out what that function exactly does I leave the question open by writing remaining ambiguous about it. For the idea of the answer it is not very important. $\endgroup$ Commented Oct 17, 2023 at 7:35
  • $\begingroup$ (1) Thank you for the answer. I understand that the non-intuitive appearing differences between the ANOVA and post/hoc test results are due to the different "recipes" behind them, and due to problems of p-values and using cutoff values like 0.05. $\endgroup$ Commented Oct 17, 2023 at 8:46
  • $\begingroup$ (2) I remain confused when it comes to the interpretation and presentation of these results. I could not find any papers (ecology and related fields) that present these kinds of results. Based on the graphical representation of the data and the results of the post-hoc tests, I am inclined to conclude that the effect of the two invasive species is different at low cover, and that this difference disappeared at high cover. But I have no idea how to present this in an article... $\endgroup$ Commented Oct 17, 2023 at 8:46
  • 1
    $\begingroup$ @benjaminjarcuska the data and image do indeed indicate that you have this interaction, but the effect is not very statistically significant. You can present it in the way that you would normally present an observed effect with a p-value of 0.0913297 $\endgroup$ Commented Oct 17, 2023 at 10:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.