plant-growth-no-insects.Rmd

---
title: "Plant growth in the absence of insects"
author: "Matthew A. Barbour"
date: "`r Sys.Date()`"
output: workflowr::wflow_html
editor_options:
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE)

# load required libraries
library(tidyverse)
library(kableExtra)
library(clubSandwich)
library(cowplot)

# set ggplot theme
theme_set(theme_cowplot())
```

# Setup

```{r}
# load data
ChamberNoInsectsDF <- read_csv("data/PreExperimentNoInsectsPlantBiomass.csv") %>%
  mutate(Cage = as.character(Cage),
         Pot = as.character(Pot))

# conduct analyses at cage level
CageLevelBiomass <- ChamberNoInsectsDF %>%
  # sum biomass across both pots
  group_by(Cage, Temperature, Richness, Composition, Col, gsm1, AOP2, AOP2.gsoh) %>%
  summarise_at(vars(Biomass_g), list(sum)) %>%
  # tidy data
  ungroup() %>%
  select(cage = Cage, temp = Temperature, rich = Richness, com = Composition, Col, gsm1, AOP2, AOP2.gsoh, Biomass_g) %>%
  # adjust temp and rich so effect of +1 C is comparable to +1 genotype
  mutate(temp = ifelse(temp == "20 C", 0, 3),
         rich = rich - 1,
         # define orthogonal constrasts to test for above-average allele effects.
        # aop2_vs_AOP2 must be included first
        aop2_vs_AOP2 = Col + gsm1 - AOP2 - AOP2.gsoh,
        mam1_vs_MAM1 = gsm1 - Col, # aop2_vs_AOP2 must be included in model
        gsoh_vs_GSOH = AOP2.gsoh - AOP2)

# source in ANOVA GLM for adjusted F-tests
source('code/glm-ftest.R')
```

# Show equivalence of Analysis of deviance and ANOVA

An Analysis of deviance on a GLM with a gaussian error distribution is equivalent to ANOVA. However, unadjusted *F*-tests are inappropriate because all terms are tested against residual variation rather than the intended error level (e.g. `com` for `rich`). The analysis below is just to prove this equivalence. I'm doing this so I can use the same function `glm.ftest.v2` for the ANOVA in the following section.

```{r}
glm.ftest.v2(
  model = glm(data = CageLevelBiomass,
              family = gaussian(link = "identity"),
              # logging improves residual distribution
              formula = log(Biomass_g) ~ temp + rich + com + temp:rich + temp:com),
 test.formula = list(c("temp","temp:com"),
                     c("rich","com"),
                     c("temp:rich","temp:com")))[[1]]

anova(aov(log(Biomass_g) ~ temp + rich + com + temp:rich + temp:com, CageLevelBiomass))
```

# Reproduce Table S5

```{r}
# fit ANOVA
biomass_noinsects_glmf <- glm.ftest.v2(
  model = glm(data = CageLevelBiomass,
              family = gaussian(link = "identity"),
              # logging improves residual distribution
              formula = log(Biomass_g) ~ temp + rich + aop2_vs_AOP2 + mam1_vs_MAM1 + gsoh_vs_GSOH + com + temp:(rich + aop2_vs_AOP2 + mam1_vs_MAM1 + gsoh_vs_GSOH) + temp:com),
 test.formula = list(c("temp","temp:com"),
                     c("rich","com"),
                     c("aop2_vs_AOP2","com"),
                     c("mam1_vs_MAM1","com"),
                     c("gsoh_vs_GSOH","com"),
                     c("temp:rich","temp:com"),
                     c("temp:aop2_vs_AOP2","temp:com"),
                     c("temp:mam1_vs_MAM1","temp:com"),
                     c("temp:gsoh_vs_GSOH","temp:com")))[[3]] %>%
  # tidy table
  select(Source = treatment, 
         `df (Source)` = num_df, 
         `df (Error)` = den_df,
         Deviance = deviance,
         `Mean Deviance` = mean_deviance,
         F = F, P = P, Error = error)

# reproduce table S5 in Supplementary Materials
biomass_noinsects_glmf %>%
  kable(., caption = "Analysis of variance for plant biomass (log transformed) in the absence of insects.", booktabs = T) %>%
  kable_styling(latex_options = c("striped", "hold_position"))
```

# Reproduce Fig. S6

```{r}
# calculate 95% confidence intervals with `com` as the cluster level
aop2_CI <- conf_int(
  glm(data = CageLevelBiomass,
              family = gaussian(link = "identity"),
              formula = log(Biomass_g) ~ -1 + temp + I(AOP2 + AOP2.gsoh) + I(Col + gsm1)),
  vcov = "CR2", 
  test = "naive-t", 
  coefs = c("I(AOP2 + AOP2.gsoh)","I(Col + gsm1)"), 
  cluster = CageLevelBiomass$com) %>%
  data.frame() %>%
  rownames_to_column(var = "term") %>%
  mutate(allele = c("AOP2","aop2"))
# note that I back transform to original scale for plotting
exp(aop2_CI$beta[2])

# get the effect of each genotype
mean_geno <- conf_int(
  glm(data = CageLevelBiomass,
              family = gaussian(link = "identity"),
              formula = log(Biomass_g) ~ -1 + temp + AOP2 + AOP2.gsoh + Col + gsm1),
  vcov = "CR2",
  test = "naive-t",
  cluster = CageLevelBiomass$com,
  coefs = c("AOP2","AOP2.gsoh","Col","gsm1")
) %>%
  data.frame() %>%
  rownames_to_column(var = "term") %>%
  mutate(allele = c("AOP2","AOP2","aop2","aop2"),
         term = factor(term, levels = c("Col","gsm1","AOP2","AOP2.gsoh"), labels = c("Col","gsm1","AOP2","AOP2/gsoh")))

# plot on original scale
# adding a genotype with an aop2 allele to the population doubles the likelihood of species persistence
ggplot(aop2_CI, aes(x = allele, y = exp(beta))) +
  geom_point(size = 5) +
  geom_point(data = mean_geno, aes(color = term), size = 5, position = position_dodge(width = 0.3)) +
  geom_linerange(aes(ymax = exp(beta + SE), ymin = exp(beta - SE)), size = 1.5) +
  geom_linerange(aes(ymax = exp(CI_U), ymin = exp(CI_L))) +
  scale_x_discrete(labels = c("AOP2\u2013","AOP2+")) +
  scale_y_continuous("Plant biomass (g)") +
  xlab("Allele") + 
  scale_color_manual(values = c("darkgreen","steelblue","darkorange","firebrick1"), name = "") + 
  theme_cowplot(font_size = 18, line_size = 1)
# ggsave(filename = "figures/AOP2-growth-no-insects.pdf", height = 6, width = 8, device=cairo_pdf)
```


# Save analysis

Write out an `.RData` file to use for creating the Supplementary Material Results.

```{r}
# save.image(file = "output/plant-growth-no-insects.RData")
```