---
title: 'Code: Nonrandom sampling measures the occurrence but not the strength of a
  textbook trophic cascade'
author: "D.R. MacNulty, E.M. Brice, E.J. Larsen"
date: "2023-11-09"
output:
  pdf_document: default
  html_document: default
---

The following code and data can be used to replicate the model displayed in Figure 1 of "Nonrandom sampling measures the occurrence but not the strength of a textbook trophic cascade," which is a Technical Note in response to Painter et al.'s 2023 comment on "Sampling bias exaggerates a textbook example of a trophic cascade" (Brice et al., Ecology Letters, 2022).

# R Code  

## Install necessary packages  
```{r, eval = F}
## Install required packages (if not already installed)
install.packages("tidyverse") # for data manipulation
install.packages("mgcv") # for building the GAMM
install.packages("emmeans") # for model predictions

## Load required packages
library(tidyverse)
library(margins)
library(mgcv)
```

## Read and format data  
The dataset "CAG_data.csv" has 14,095 records of individual young aspen (plants $\ge$ 1 year-old & $\le$ 600 cm) measured from 2007-2017.

The dataframe has the following columns:  

  1) Plot: individual identifier for each of 113 plots distributed randomly across the study area. Each plot was a 1 × 20 m belt transect located randomly within an aspen stand 
  2) Year: year in which aspen was sampled
  3) Height: height (cm) of the leader stem of each individual aspen
  4) Browse: denotes the browsing status (browsed = 1,unbrowsed = 0) of the leader (tallest) stem. A leader was ‘browsed’ if its growth from the previous growing season had been eaten
  5) CAG: current annual growth of the leader stem (cm)


The following code reads and formats the aspen data for the subsequent analysis. 

```{r, eval = F}
## Set working directory
setwd()

## Read aspen data 
cag.df <- read.csv("CAG_data.csv")

# Make year a factor
cag.df$Year.f <- as.factor(cag.df$Year)

# Make plot a factor
cag.df$Plot <- as.factor(cag.df$Plot)
```

## Figure 1: Height Model  
We used a GAMM to test for an effect of leader length (current annual growth) and browsing on young aspen height. We treated the leader as the unit of analysis and used a GAMM with a gamma distribution and a log link to analyze height, which took only non-negative values that were strongly right skewed. Leader length was fit with a cubic regression smoothing function, and browse was a binary variable. We included crossed random effects for stand identity (N = 113 stands) and year to account for (i) correlation between measurements taken on the same stand in multiple years and on multiple stands in the same year, and (ii) unmeasured stand- and year-related effects. 

```{r, eval = F}
## Run model
mod <- gam(Height ~ s(CAG,k = 5,bs = "cr") + Browsed + s(Year.f,bs = "re") + s(Plot,bs = "re"),
           data=cag.df, family = Gamma(link = "log"))

# Get a summary of the model
summary(mod)

# Get coefficients for splines
coef(mod)

# To get the standard error of the splines
sqrt(diag(vcov(mod)))

# To see knot placement:
mod[['smooth']][[1]]$xp

## Plot model predictions
emmeans(mod, ~Browsed*CAG, at = list(CAG = seq(0,100,1),
                                      Browsed = c(0,1)), type = "response") %>% 
  as.data.frame(.) %>% 
  mutate(Browsed = as.factor(Browsed)) %>% 
  ggplot(aes(x = CAG))+
  geom_rug(data = cag.df %>% 
             filter(Browsed == 0),
           aes(x = CAG, y = Height),sides = "t",color = "#0072B2")+
  geom_rug(data = cag.df %>% 
             filter(Browsed == 1),
           aes(x = CAG, y = Height),sides = "b",color = "#D55E00")+
  geom_ribbon(aes(ymin = lower.CL, ymax = upper.CL, fill = Browsed), alpha = 0.3)+
  geom_line(aes(y = response, color = Browsed), linewidth=0.7)+
  scale_x_continuous(breaks=seq(0,100,20),labels=seq(0,100,20),
                     limits = c(0,105),expand = c(0,0))+
  scale_y_continuous(breaks=seq(0,200,50),labels=seq(0,200,50),
                     limits = c(0,210),expand = c(0, 0))+
  scale_fill_manual(values = c("#0072B2", "#D55E00"),labels = c("Unbrowsed","Browsed"),name = "")+
  scale_color_manual(values = c("#0072B2", "#D55E00"),labels = c("Unbrowsed","Browsed"),name = "")+
  xlab("Leader length (cm)")+
  ylab("Young aspen height (cm)")+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black"),
        axis.text = element_text(size=11, color="black"),
        legend.justification = c("right","bottom"),legend.position = c(1,0.05),
        legend.title = element_blank())
```