trillium.RdThe size of Trillium has been hypothesized to be correlated with the intensity of deer browse occuring on thee population.
trillium
A data frame
site code for different forests and woodlots in NW Pennsylvania near the Pymatuning Laboratory of Ecology
species. All are = T for Trillium
length from soil to leaf whorl
length of longest leaf in cm
Author year.
library(ggplot2) library(ggpubr) if (FALSE) { ## Explore data graphically ### Plot boxplots ggpubr::ggboxplot(data = trillium, y = "length.stem.cm", x = "site", fill = "site") ### Plot histograms ggpubr::gghistogram(data = trillium, x = "length.stem.cm", title = "All data") ggpubr::gghistogram(data = trillium, x = "length.stem.cm", facet.by = "site", fill = "site", title = "Faceted by site") ## Plot means with 95% confidence intervals ggpubr::ggerrorplot(trillium, x = "site", y = "length.stem.cm", desc_stat = "mean_ci", add = "mean") } ## 1-way ANOVA to compare sites ### null model model.null <- lm(length.stem.cm ~ 1, data = trillium) ### model of interest model.alt <- lm(length.stem.cm ~ site, data = trillium) ### compare models anova(model.null, model.alt)#> Analysis of Variance Table #> #> Model 1: length.stem.cm ~ 1 #> Model 2: length.stem.cm ~ site #> Res.Df RSS Df Sum of Sq F Pr(>F) #> 1 550 12639 #> 2 547 10012 3 2626.9 47.839 < 2.2e-16 *** #> --- #> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1## Pairwise comparisons after 1-way ANOVA ### no corrections for multiple comparisons pairwise.t.test(x = trillium$length.stem.cm, g = trillium$site, p.adjust.method = "none")#> #> Pairwise comparisons using t tests with pooled SD #> #> data: trillium$length.stem.cm and trillium$site #> #> CN ER TW #> ER 0.0036 - - #> TW 0.0046 0.4113 - #> WC 2.4e-16 4.0e-08 < 2e-16 #> #> P value adjustment method: none### Bonferonni correction pairwise.t.test(x = trillium$length.stem.cm, g = trillium$site, p.adjust.method = "bonferroni")#> #> Pairwise comparisons using t tests with pooled SD #> #> data: trillium$length.stem.cm and trillium$site #> #> CN ER TW #> ER 0.022 - - #> TW 0.028 1.000 - #> WC 1.5e-15 2.4e-07 < 2e-16 #> #> P value adjustment method: bonferroni## Tukey test ### re-fit model with aov() model.alt.aov <- aov(length.stem.cm ~ site, data = trillium) ### TukeyHSD() on model from aov() TukeyHSD(model.alt.aov)#> Tukey multiple comparisons of means #> 95% family-wise confidence level #> #> Fit: aov(formula = length.stem.cm ~ site, data = trillium) #> #> $site #> diff lwr upr p adj #> ER-CN 2.6061111 0.3061946 4.906028 0.0190462 #> TW-CN 2.0765353 0.1942317 3.958839 0.0239234 #> WC-CN 6.2315741 4.3338279 8.129320 0.0000000 #> TW-ER -0.5295758 -2.1894709 1.130319 0.8439879 #> WC-ER 3.6254630 1.9480765 5.302849 0.0000002 #> WC-TW 4.1550388 3.1220450 5.188033 0.0000000 #>### Plot effect sizes #plotTukeysHSD(TukeyHSD(model.alt.aov))