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.

Summarize data

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=