library(here)
library(tidyverse)
library(cowplot)
library(lme4)
library(car) # for Anova()
library(broom) # for tidy()
library(ggformula) # for stat_qqline()
options(scipen = 999)
library(bbmle) # for AICtab()
library(latex2exp) # for TeX
#write_bib(c("lme4", "car"), "modeling.bib")
library(emmeans)
Read in data
growth.change <- read_rds(here("cleaned_data", "growth_change.rds"))
growth.change
contains measurements of tea plant size/growth collected at initial and final time points.
Treatment
: this is what rainfall treatment the plot received.
plot
: which plot at BACE the plants were in
unique.id
: a plant ID number within the plot
H_i
and H_f
: initial and final height from the ground to the highest leaf measured in cm
d_height
: \(H_f - H_i\)
RGR_height
: relative growth rate. Unused in analyses
L_i
and L_f
: Initial and final leaf count
d_leaf
: \(L_f - L_i\)
RGR_leaf
: relative growth rate for leaf number. Unused in analyses
shoots_f
: the number of shoots at the final measurement
Growth (change in height)
I tried using RGR, but it created problems due to non-normality of residuals, and difficulty transforming data due to negative and 0 values. All analyses use d_height
instead.
Model selection
height.lmer <- lmer(d_height ~ Treatment + (1|plot), data = growth.change)
height.lm <- lm(d_height ~ Treatment, data = growth.change)
AICtab(height.lmer, height.lm)
dAIC df
height.lmer 0.0 5
height.lm 2.1 4
Simple lm wins because of fewer df and dAIC of only 2.
Check residuals
augment(height.lmer) %>%
ggplot(aes(sample = .resid)) +
geom_qq() +
stat_qqline() +
coord_flip()

Looks fine.
ANOVA
Anova(height.lm)
Anova Table (Type II tests)
Response: d_height
Sum Sq Df F value Pr(>F)
Treatment 127.4 2 3.7723 0.02685 *
Residuals 1469.1 87
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey’s Tests
(height.posthoc <- emmeans(height.lm, pairwise ~ Treatment))
$emmeans
Treatment emmean SE df lower.CL upper.CL
100% 5.833333 0.750249 87 4.342132 7.324534
75% 3.533333 0.750249 87 2.042132 5.024534
50% 3.133333 0.750249 87 1.642132 4.624534
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
100% - 75% 2.3 1.061012 87 2.168 0.0825
100% - 50% 2.7 1.061012 87 2.545 0.0336
75% - 50% 0.4 1.061012 87 0.377 0.9247
P value adjustment: tukey method for comparing a family of 3 estimates
height.letters <- growth.change %>%
group_by(Treatment) %>%
summarise(y = max(d_height)) %>%
add_column(letters = c("a", "b", "b"))
height.errbar <- height.posthoc$emmeans %>%
tidy()
Leaf Count
Model selection
leaf.d.lmer <- lmer(d_leaf ~ Treatment + (1|plot), data = growth.change)
leaf.d.lm <- lm(d_leaf ~ Treatment, data = growth.change)
AICtab(leaf.d.lm, leaf.d.lmer)
dAIC df
leaf.d.lmer 0.0 5
leaf.d.lm 9.8 4
lmer wins.
Check Residuals
augment(leaf.d.lmer) %>%
ggplot(aes(sample = .resid)) +
geom_qq()+
stat_qqline() +
coord_flip()

Looks like it would be improved by sqrt transformation, but there are negative numbers. There are also 0’s, so no sure about log transformation. I think it looks good enough, and playing around with sqrt() transformation (which induces NaNs) makes me think it won’t change the result.
ANOVA
Anova(leaf.d.lmer)
Analysis of Deviance Table (Type II Wald chisquare tests)
Response: d_leaf
Chisq Df Pr(>Chisq)
Treatment 10.205 2 0.006081 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Tukey’s Tests
(leaf.posthoc <- emmeans(leaf.d.lmer, pairwise ~ Treatment))
$emmeans
Treatment emmean SE df lower.CL upper.CL
100% 25.50000 3.159055 6 17.770071 33.22993
75% 18.03333 3.159055 6 10.303404 25.76326
50% 11.23333 3.159055 6 3.503404 18.96326
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
100% - 75% 7.466667 4.467579 6 1.671 0.2900
100% - 50% 14.266667 4.467579 6 3.193 0.0429
75% - 50% 6.800000 4.467579 6 1.522 0.3464
P value adjustment: tukey method for comparing a family of 3 estimates
leaf.letters <- growth.change %>%
group_by(Treatment) %>%
summarise(y = max(d_leaf)) %>%
add_column(letters = c("a", "ab", "b"))
leaf.errbar <- leaf.posthoc$emmeans %>%
tidy()
Plots
d_height_plot <- ggplot(growth.change, aes(x = Treatment, y = d_height)) +
geom_jitter(width = 0.2, alpha = 0.5) +
geom_point(aes(y = estimate), shape = 15, size = 3, data = height.errbar) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high, y = estimate), width = 0.2, data = height.errbar) +
geom_text(aes(y = y, label = letters), vjust = -1, data = height.letters) +
geom_label(aes(y = 19, x = 3, label = paste0("p = ", Anova(height.lm) %>% tidy() %>% .[1,5] %>% round(3)))) +
ylim(-9, 19) +
ylab(TeX("$\\Delta$ Height (cm)")) +
xlab("% Ambient Rainfall") +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
d_leaf_plot <- ggplot(growth.change, aes(x = Treatment, y = d_leaf)) +
geom_jitter(width = 0.2, alpha = 0.5) +
geom_point(aes(y = estimate), shape = 15, size = 3, data = leaf.errbar) +
geom_errorbar(aes(ymin = conf.low, ymax = conf.high, y = estimate), width = 0.2, data = leaf.errbar) +
geom_text(aes(y = y, label = letters), vjust = -1, data = leaf.letters) +
geom_label(aes(y = 80, x = 3, label = paste0("p = ", Anova(leaf.d.lmer)$`Pr(>Chisq)` %>% round(3)))) +
ylim(-10, 80) +
ylab(TeX("$\\Delta$ Leaf Number")) +
xlab("% Ambient Rainfall") +
theme(axis.title = element_text(size = 10),
axis.text = element_text(size = 9))
growth_plots <- plot_grid(d_height_plot, d_leaf_plot,
ncol = 2, nrow = 1, align = "h",
labels = "AUTO")
growth_plots
save_plot(here("figs", "growth plots.png"), growth_plots,
ncol = 2, nrow = 1,
base_width = 3)
LS0tCnRpdGxlOiAiZ3Jvd3RoIGFuYWx5c2lzIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCmBgYHtyIHBhY2thZ2VzLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KGhlcmUpCmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGNvd3Bsb3QpCmxpYnJhcnkobG1lNCkKbGlicmFyeShjYXIpICMgZm9yIEFub3ZhKCkKbGlicmFyeShicm9vbSkgIyBmb3IgdGlkeSgpCmxpYnJhcnkoZ2dmb3JtdWxhKSAjIGZvciBzdGF0X3FxbGluZSgpCm9wdGlvbnMoc2NpcGVuID0gOTk5KQpsaWJyYXJ5KGJibWxlKSAjIGZvciBBSUN0YWIoKQpsaWJyYXJ5KGxhdGV4MmV4cCkgIyBmb3IgVGVYCiN3cml0ZV9iaWIoYygibG1lNCIsICJjYXIiKSwgIm1vZGVsaW5nLmJpYiIpCmxpYnJhcnkoZW1tZWFucykKYGBgCgojIFJlYWQgaW4gZGF0YQpgYGB7cn0KZ3Jvd3RoLmNoYW5nZSA8LSByZWFkX3JkcyhoZXJlKCJjbGVhbmVkX2RhdGEiLCAiZ3Jvd3RoX2NoYW5nZS5yZHMiKSkKYGBgCmBncm93dGguY2hhbmdlYCBjb250YWlucyBtZWFzdXJlbWVudHMgb2YgdGVhIHBsYW50IHNpemUvZ3Jvd3RoIGNvbGxlY3RlZCBhdCBpbml0aWFsIGFuZCBmaW5hbCB0aW1lIHBvaW50cy4KCi0gYFRyZWF0bWVudGA6IHRoaXMgaXMgd2hhdCByYWluZmFsbCB0cmVhdG1lbnQgdGhlIHBsb3QgcmVjZWl2ZWQuCi0gYHBsb3RgOiB3aGljaCBwbG90IGF0IEJBQ0UgdGhlIHBsYW50cyB3ZXJlIGluCi0gYHVuaXF1ZS5pZGA6IGEgcGxhbnQgSUQgbnVtYmVyIHdpdGhpbiB0aGUgcGxvdAotIGBIX2lgIGFuZCBgSF9mYDogaW5pdGlhbCBhbmQgZmluYWwgaGVpZ2h0IGZyb20gdGhlIGdyb3VuZCB0byB0aGUgaGlnaGVzdCBsZWFmIG1lYXN1cmVkIGluIGNtCi0gYGRfaGVpZ2h0YDogJEhfZiAtIEhfaSQKLSBgUkdSX2hlaWdodGA6IHJlbGF0aXZlIGdyb3d0aCByYXRlLiAgVW51c2VkIGluIGFuYWx5c2VzCi0gYExfaWAgYW5kIGBMX2ZgOiBJbml0aWFsIGFuZCBmaW5hbCBsZWFmIGNvdW50Ci0gYGRfbGVhZmA6ICRMX2YgLSBMX2kkCi0gYFJHUl9sZWFmYDogcmVsYXRpdmUgZ3Jvd3RoIHJhdGUgZm9yIGxlYWYgbnVtYmVyLiAgVW51c2VkIGluIGFuYWx5c2VzCi0gYHNob290c19mYDogdGhlIG51bWJlciBvZiBzaG9vdHMgYXQgdGhlIGZpbmFsIG1lYXN1cmVtZW50CgojIFN1bW1hcml6ZSBkYXRhCgpgYGB7cn0KZ3Jvd3RoLmNoYW5nZSAlPiUgCiAgZ3JvdXBfYnkoVHJlYXRtZW50KSAlPiUgCiAgc3VtbWFyaXplKEhfaS5tZWFuID0gbWVhbihIX2kpLAogICAgICAgICAgICBIX2kuc2QgPSBzZChIX2kpLAogICAgICAgICAgICBIX2YubWVhbiA9IG1lYW4oSF9mKSwKICAgICAgICAgICAgSF9mLnNkID0gc2QoSF9mKSwKICAgICAgICAgICAgTF9pLm1lYW4gPSBtZWFuKExfaSksCiAgICAgICAgICAgIExfaS5zZCA9IHNkKExfaSksCiAgICAgICAgICAgIExfZi5tZWFuID0gbWVhbihMX2YpLAogICAgICAgICAgICBMX2Yuc2QgPSBzZChMX2YpKSAlPiUgCiAgbXV0YXRlX2lmKGlzLmRvdWJsZSwgfnJvdW5kKC4sMSkpCgpncm93dGguY2hhbmdlICU+JSAKICBzdW1tYXJpemUoc2hvb3QubWVhbiA9IG1lYW4oc2hvb3RzX2YpLAogICAgICAgICAgICBzaG9vdC5zZCA9IHNkKHNob290c19mKSkKCmdyb3d0aC5jaGFuZ2UgJT4lIAogICMgZ3JvdXBfYnkoVHJlYXRtZW50KSAlPiUgCiAgc3VtbWFyaXplKG1lYW5sZWFmID0gbWVhbihMX2YpLAogICAgICAgICAgICBzZGxlYWYgPSBzZChMX2YpKQpgYGAKCgojIEdyb3d0aCAoY2hhbmdlIGluIGhlaWdodCkKSSB0cmllZCB1c2luZyBSR1IsIGJ1dCBpdCBjcmVhdGVkIHByb2JsZW1zIGR1ZSB0byBub24tbm9ybWFsaXR5IG9mIHJlc2lkdWFscywgYW5kIGRpZmZpY3VsdHkgdHJhbnNmb3JtaW5nIGRhdGEgZHVlIHRvIG5lZ2F0aXZlIGFuZCAwIHZhbHVlcy4gIEFsbCBhbmFseXNlcyB1c2UgYGRfaGVpZ2h0YCBpbnN0ZWFkLgoKCiMjIE1vZGVsIHNlbGVjdGlvbgpgYGB7cn0KaGVpZ2h0LmxtZXIgPC0gbG1lcihkX2hlaWdodCB+IFRyZWF0bWVudCArICgxfHBsb3QpLCBkYXRhID0gZ3Jvd3RoLmNoYW5nZSkKaGVpZ2h0LmxtIDwtIGxtKGRfaGVpZ2h0IH4gVHJlYXRtZW50LCBkYXRhID0gZ3Jvd3RoLmNoYW5nZSkKQUlDdGFiKGhlaWdodC5sbWVyLCBoZWlnaHQubG0pCmBgYApTaW1wbGUgbG0gd2lucyBiZWNhdXNlIG9mIGZld2VyIGRmIGFuZCBkQUlDIG9mIG9ubHkgMi4KCiMjIENoZWNrIHJlc2lkdWFscwoKYGBge3J9CmF1Z21lbnQoaGVpZ2h0LmxtZXIpICU+JSAKICBnZ3Bsb3QoYWVzKHNhbXBsZSA9IC5yZXNpZCkpICsgCiAgZ2VvbV9xcSgpICsKICBzdGF0X3FxbGluZSgpICsKICBjb29yZF9mbGlwKCkKYGBgCkxvb2tzIGZpbmUuCgojIyBBTk9WQQpgYGB7cn0KQW5vdmEoaGVpZ2h0LmxtKQpgYGAKCiMjIFR1a2V5J3MgVGVzdHMKYGBge3J9CihoZWlnaHQucG9zdGhvYyA8LSBlbW1lYW5zKGhlaWdodC5sbSwgcGFpcndpc2UgfiBUcmVhdG1lbnQpKQpoZWlnaHQubGV0dGVycyA8LSBncm93dGguY2hhbmdlICU+JQogIGdyb3VwX2J5KFRyZWF0bWVudCkgJT4lCiAgc3VtbWFyaXNlKHkgPSBtYXgoZF9oZWlnaHQpKSAlPiUKICBhZGRfY29sdW1uKGxldHRlcnMgPSBjKCJhIiwgImIiLCAiYiIpKQoKaGVpZ2h0LmVycmJhciA8LSBoZWlnaHQucG9zdGhvYyRlbW1lYW5zICU+JQogIHRpZHkoKQpgYGAKCiMgTGVhZiBDb3VudAojIyBNb2RlbCBzZWxlY3Rpb24KYGBge3J9CmxlYWYuZC5sbWVyIDwtIGxtZXIoZF9sZWFmIH4gVHJlYXRtZW50ICsgKDF8cGxvdCksIGRhdGEgPSBncm93dGguY2hhbmdlKQpsZWFmLmQubG0gPC0gbG0oZF9sZWFmIH4gVHJlYXRtZW50LCBkYXRhID0gZ3Jvd3RoLmNoYW5nZSkKQUlDdGFiKGxlYWYuZC5sbSwgbGVhZi5kLmxtZXIpCmBgYApsbWVyIHdpbnMuCgojIyBDaGVjayBSZXNpZHVhbHMKYGBge3J9CmF1Z21lbnQobGVhZi5kLmxtZXIpICU+JSAKICBnZ3Bsb3QoYWVzKHNhbXBsZSA9IC5yZXNpZCkpICsKICBnZW9tX3FxKCkrCiAgc3RhdF9xcWxpbmUoKSArCiAgY29vcmRfZmxpcCgpCmBgYApMb29rcyBsaWtlIGl0IHdvdWxkIGJlIGltcHJvdmVkIGJ5IHNxcnQgdHJhbnNmb3JtYXRpb24sIGJ1dCB0aGVyZSBhcmUgbmVnYXRpdmUgbnVtYmVycy4gIFRoZXJlIGFyZSBhbHNvIDAncywgc28gbm8gc3VyZSBhYm91dCBsb2cgdHJhbnNmb3JtYXRpb24uICBJIHRoaW5rIGl0IGxvb2tzIGdvb2QgZW5vdWdoLCBhbmQgcGxheWluZyBhcm91bmQgd2l0aCBzcXJ0KCkgdHJhbnNmb3JtYXRpb24gKHdoaWNoIGluZHVjZXMgTmFOcykgbWFrZXMgbWUgdGhpbmsgaXQgd29uJ3QgY2hhbmdlIHRoZSByZXN1bHQuCgojIyBBTk9WQQpgYGB7cn0KQW5vdmEobGVhZi5kLmxtZXIpIApgYGAKIyMgVHVrZXkncyBUZXN0cwpgYGB7cn0KKGxlYWYucG9zdGhvYyA8LSBlbW1lYW5zKGxlYWYuZC5sbWVyLCBwYWlyd2lzZSB+IFRyZWF0bWVudCkpCmxlYWYubGV0dGVycyA8LSBncm93dGguY2hhbmdlICU+JQogIGdyb3VwX2J5KFRyZWF0bWVudCkgJT4lIAogIHN1bW1hcmlzZSh5ID0gbWF4KGRfbGVhZikpICU+JQogIGFkZF9jb2x1bW4obGV0dGVycyA9IGMoImEiLCAiYWIiLCAiYiIpKQoKbGVhZi5lcnJiYXIgPC0gbGVhZi5wb3N0aG9jJGVtbWVhbnMgJT4lIAogIHRpZHkoKQpgYGAKCiMgUGxvdHMKCmBgYHtyfQpkX2hlaWdodF9wbG90IDwtIGdncGxvdChncm93dGguY2hhbmdlLCBhZXMoeCA9IFRyZWF0bWVudCwgeSA9IGRfaGVpZ2h0KSkgKwogIGdlb21faml0dGVyKHdpZHRoID0gMC4yLCBhbHBoYSA9IDAuNSkgKwogIGdlb21fcG9pbnQoYWVzKHkgPSBlc3RpbWF0ZSksIHNoYXBlID0gMTUsIHNpemUgPSAzLCBkYXRhID0gaGVpZ2h0LmVycmJhcikgKwogIGdlb21fZXJyb3JiYXIoYWVzKHltaW4gPSBjb25mLmxvdywgeW1heCA9IGNvbmYuaGlnaCwgeSA9IGVzdGltYXRlKSwgd2lkdGggPSAwLjIsIGRhdGEgPSBoZWlnaHQuZXJyYmFyKSArCiAgZ2VvbV90ZXh0KGFlcyh5ID0geSwgbGFiZWwgPSBsZXR0ZXJzKSwgdmp1c3QgPSAtMSwgZGF0YSA9IGhlaWdodC5sZXR0ZXJzKSArCiAgZ2VvbV9sYWJlbChhZXMoeSA9IDE5LCB4ID0gMywgbGFiZWwgPSBwYXN0ZTAoInAgPSAiLCBBbm92YShoZWlnaHQubG0pICU+JSB0aWR5KCkgJT4lIC5bMSw1XSAlPiUgcm91bmQoMykpKSkgKwogIHlsaW0oLTksIDE5KSArCiAgeWxhYihUZVgoIiRcXERlbHRhJCBIZWlnaHQgKGNtKSIpKSArCiAgeGxhYigiJSBBbWJpZW50IFJhaW5mYWxsIikgKwogIHRoZW1lKGF4aXMudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEwKSwKICAgICAgICBheGlzLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDkpKQogIApkX2xlYWZfcGxvdCA8LSBnZ3Bsb3QoZ3Jvd3RoLmNoYW5nZSwgYWVzKHggPSBUcmVhdG1lbnQsIHkgPSBkX2xlYWYpKSArCiAgZ2VvbV9qaXR0ZXIod2lkdGggPSAwLjIsIGFscGhhID0gMC41KSArIAogIGdlb21fcG9pbnQoYWVzKHkgPSBlc3RpbWF0ZSksIHNoYXBlID0gMTUsIHNpemUgPSAzLCBkYXRhID0gbGVhZi5lcnJiYXIpICsKICBnZW9tX2Vycm9yYmFyKGFlcyh5bWluID0gY29uZi5sb3csIHltYXggPSBjb25mLmhpZ2gsIHkgPSBlc3RpbWF0ZSksIHdpZHRoID0gMC4yLCBkYXRhID0gbGVhZi5lcnJiYXIpICsKICBnZW9tX3RleHQoYWVzKHkgPSB5LCBsYWJlbCA9IGxldHRlcnMpLCB2anVzdCA9IC0xLCBkYXRhID0gbGVhZi5sZXR0ZXJzKSArCiAgZ2VvbV9sYWJlbChhZXMoeSA9IDgwLCB4ID0gMywgbGFiZWwgPSBwYXN0ZTAoInAgPSAiLCBBbm92YShsZWFmLmQubG1lcikkYFByKD5DaGlzcSlgICU+JSByb3VuZCgzKSkpKSArCiAgeWxpbSgtMTAsIDgwKSArCiAgeWxhYihUZVgoIiRcXERlbHRhJCBMZWFmIE51bWJlciIpKSArCiAgeGxhYigiJSBBbWJpZW50IFJhaW5mYWxsIikgKwogIHRoZW1lKGF4aXMudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEwKSwKICAgICAgICBheGlzLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDkpKQoKZ3Jvd3RoX3Bsb3RzIDwtIHBsb3RfZ3JpZChkX2hlaWdodF9wbG90LCBkX2xlYWZfcGxvdCwKICAgICAgICAgICAgICAgICAgICAgIG5jb2wgPSAyLCBucm93ID0gMSwgYWxpZ24gPSAiaCIsCiAgICAgICAgICAgICAgICAgICAgICBsYWJlbHMgPSAiQVVUTyIpCmdyb3d0aF9wbG90cwpzYXZlX3Bsb3QoaGVyZSgiZmlncyIsICJncm93dGggcGxvdHMucG5nIiksIGdyb3d0aF9wbG90cywKICAgICAgICAgIG5jb2wgPSAyLCBucm93ID0gMSwKICAgICAgICAgIGJhc2Vfd2lkdGggPSAzKQpgYGA=