---
title: "Explaining analyses from the Montana lake trout food web study"
author: "Charles Wainright"
date: "December 9, 2020"
output: html_document
---

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

This document explains analyses included in the Montana lake trout food web study (hereafter called "the article") to enhance reproducibility of these techniques. Questions regarding these analyses or the R markdown should be directed to Charles Wainright at charles_wainright@fws.gov.

```{r, echo=TRUE, results=FALSE, warning=FALSE, message=FALSE}
library(ggplot2)
library(dplyr)
library(ellipse)
library(gridExtra)
library(rjags)
library(hdrcde)
library(reshape2)
library(tidyr)
library(cowplot)
library(lme4)
library(effects)
library(lmerTest)
library(SIBER)
library(DataCombine)
library(tidyverse)
library(Hmisc)
library(nicheROVER)
library(MixSIAR)
library(cowplot)
library(knitr)
```



###Dataset overview

```{r,echo=FALSE}
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
```

Our dataset has 47 columns and 1459 rows. Each row is one isotope sample. There are 10 lakes and 437 fish samples in this dataset.

```{r}
length(isotope_dataset) #count of columns (variables)
length(isotope_dataset$d13C) #count of rows (samples)
length(unique(isotope_dataset$Lake_name)) #count of lakes
length(subset(isotope_dataset,Sample_type=="Fish")$d13C) #count of fish samples
```
The 4 columns included in the overview below (Table 1) show the most important features of the dataset, including lake name, the scientific name of the organism associated with each sample, and the absolute $\delta$^13^C and $\delta$^15^N for each sample.
```{r, echo = TRUE}
data <- subset(isotope_dataset, select=c(Lake_name,Scientific_name,d13C,d15N))
colnames(data)[1] <- "Lake name"
colnames(data)[2] <- "Scientific name"
colnames(data)[3] <- "$\\delta$^13^C ($\\unicode{x2030}$)"
colnames(data)[4] <- "$\\delta$^15^N ($\\unicode{x2030}$)"
library(knitr)
kable(head(data), align="llcc", caption = "Table 1. Dataset overview.")
```

###Trophic baseline selection

Our first task was to account for spatial and temporal isotopic variability among our study lakes (Post 2002, Guzzo et. al., 2011) by selecting the best $\delta$^13^C and $\delta$^15^N isotopic baseline. To accomplish this task, we measured $\delta$^13^C and $\delta$^15^N variability among primary consumers in our study system and then chose the least variable consumer (Syväranta and Jones 2008) for later isotopic baseline corrections. Since standard error scales inversely with sample size, we started with our dataset's most abundant invertebrates and periphyton as trophic baseline candidates. These candidates included mayflies (Ephemeroptera: Baetidae, Leptophlebiidae, Heptageniidae, and Ephemerellidae), lithophytic periphyton from 0.5 m depth, caddisflies (Trichoptera: Limnephilidae and Leptoceridae), water mites (Hydrachnidia), bulk zooplankton from 10 m depth, and snails (Gastropoda: Physidae, Lynmaeidae, Planorbidae). 

The following code chunks summarize how we selected a trophic baseline for our study lakes. First, we calculated within-lake $\delta$^13^C and $\delta$^15^N variability for each candidate. We provide periphyton as an example of how within-lake $\delta$^13^C and $\delta$^15^N variability was determined for each trophic baseline candidate (Table 2). 

```{r, echo=TRUE}
####################################################################
#Periphyton summary
####################################################################
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_09102020.csv", header=TRUE)
data <- isotope_dataset
data <- data[c("Lake_name","real.name","d13C","d15N")] #$real.name is for convenience because it uses scientific name for invertebrates but uses common name for fish
data <- subset(data,data$real.name%in%c("Periphyton"))
#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Lake_name", "real.name"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C", "d15N", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         n = length(d13C),
                         mean_d13C = mean(d13C),
                         sd_d13C = sd(d13C),
                         se_d13C = sqrt(var(d13C)/length(d13C)),
                         mean_d15N = mean(d15N),
                         sd_d15N = sd(d15N),
                         se_d15N = sqrt(var(d15N)/length(d15N))
  )
 if (printSummary)  print(summary)
  summary
  return(summary)
}
#create summary dataframe using function and grouping specified above
periphyton_isotope_summary<-summariseIsotopeData(df = subset(data),
                                                   grouping = c("Lake_name", "real.name"))
#round numbers
periphyton_isotope_summary$mean_d13C <- round(periphyton_isotope_summary$mean_d13C,2)
periphyton_isotope_summary$se_d13C <- round(periphyton_isotope_summary$se_d13C,2)
periphyton_isotope_summary$mean_d15N <- round(periphyton_isotope_summary$mean_d15N,2)
periphyton_isotope_summary$se_d15N <- round(periphyton_isotope_summary$se_d15N,2)
periphyton_isotope_summary<-subset(periphyton_isotope_summary,select=-c(sd_d13C,sd_d15N))
#re-name columns
colnames(periphyton_isotope_summary)[1] <- "Lake name"
colnames(periphyton_isotope_summary)[2] <- "Candidate"
colnames(periphyton_isotope_summary)[4] <- "Mean $\\delta$^13^C ($\\unicode{x2030}$)"
colnames(periphyton_isotope_summary)[5] <- "SE $\\delta$^13^C ($\\unicode{x2030}$)"
colnames(periphyton_isotope_summary)[6] <- "Mean $\\delta$^15^N ($\\unicode{x2030}$)"
colnames(periphyton_isotope_summary)[7] <- "SE $\\delta$^15^N ($\\unicode{x2030}$)"
#make table
library(knitr)
kable(periphyton_isotope_summary, align="lcccccc", caption = "Table 2. Within-lake isotopic variability of periphyton.")
```

We repeated the within-lake variability process for each trophic baseline candidate. Then, we took the mean of each candidate's within-lake $\delta$^13^C and $\delta$^15^N standard error to determine its among-lake isotopic variability. We repeated this process to determine among-lake $\delta$^13^C and $\delta$^15^N variability for each candidate (Table 3). Among-lake isotopic variability was the basis for our trophic baseline selection because among-lake isotopic variability tells us which trophic baseline candidate had the most consistent $\delta$^13^C and $\delta$^15^N signatures and, therefore, which candidate's sample mean ($\bar{x}$) is most likely to be nearest to its population mean ($\mu$).

```{r, echo=FALSE}
rm(list=ls())
graphics.off()
set.seed(1)
####################################################################
#Periphyton summary
####################################################################
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_09102020.csv", header=TRUE)
data <- isotope_dataset
data <- data[c("Lake_name","real.name","d13C","d13C_corr","d15N","d15N_corr")]
data <- subset(data,data$real.name%in%c("Periphyton"))
#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Lake_name", "real.name"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C", "d15N", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         n = length(d13C),
                         mean_d13C = mean(d13C),
                         sd_d13C = sd(d13C),
                         se_d13C = sqrt(var(d13C)/length(d13C)),
                         mean_d15N = mean(d15N),
                         sd_d15N = sd(d15N),
                         se_d15N = sqrt(var(d15N)/length(d15N))
  )
 if (printSummary)  print(summary)
  summary
  return(summary)
}
#create summary dataframe using function and grouping specified above
periphyton_isotope_summary<-summariseIsotopeData(df = subset(data),
                                                   grouping = c("Lake_name", "real.name"))

####################################################################
#Zooplankton summary
####################################################################
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_09102020.csv", header=TRUE)
data <- isotope_dataset
data <- data[c("Lake_name","real.name","d13C","d13C_corr","d15N","d15N_corr")]
data <- subset(data,data$real.name%in%c("Zooplankton"))
#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Lake_name", "real.name"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C", "d15N", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         n = length(d13C),
                         mean_d13C = mean(d13C),
                         sd_d13C = sd(d13C),
                         se_d13C = sqrt(var(d13C)/length(d13C)),
                         mean_d15N = mean(d15N),
                         sd_d15N = sd(d15N),
                         se_d15N = sqrt(var(d15N)/length(d15N))
  )
  
  if (printSummary)  print(summary)
  
  summary
  return(summary)
}
#create summary dataframe using function and grouping specified above
zooplankton_isotope_summary<-summariseIsotopeData(df = subset(data),
                                                   grouping = c("Lake_name", "real.name"))
####################################################################
#Hydrachnidia summary
####################################################################
#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_09102020.csv", header=TRUE)
data <- isotope_dataset
data <- data[c("Lake_name","real.name","d13C","d13C_corr","d15N","d15N_corr")]
data <- subset(data,data$real.name%in%c("Hydrachnidia"))
#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Lake_name", "real.name"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C", "d15N", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         n = length(d13C),
                         mean_d13C = mean(d13C),
                         sd_d13C = sd(d13C),
                         se_d13C = sqrt(var(d13C)/length(d13C)),
                         mean_d15N = mean(d15N),
                         sd_d15N = sd(d15N),
                         se_d15N = sqrt(var(d15N)/length(d15N))
  )
  
  if (printSummary)  print(summary)
  summary
  return(summary)
}
#create summary dataframe using function and grouping specified above
hydrachnidia_isotope_summary<-summariseIsotopeData(df = subset(data),
                                         grouping = c("Lake_name", "real.name"))

####################################################################
#Caddis summary
####################################################################
caddis_lookup <- data.frame(Lake_name=c("Big Salmon","Bowman","Kintla","Trout","Logging","Quartz","Hungry Horse","McDonald","Swan","Lindbergh"),
                            n_base=rep(NA,10),
                            d13C_baseline=rep(NA,10),
                            d13C_base_se=rep(NA,10),
                            d15N_baseline=rep(NA,10),
                            d15N_base_se=rep(NA,10))
#1/10: Big salmon
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Big Salmon")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[1] <- length(data$real.name)
caddis_lookup$d13C_baseline[1] <- mean(data$d13C)
caddis_lookup$d13C_base_se[1] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[1] <- mean(data$d15N)
caddis_lookup$d15N_base_se[1] <- sqrt(var(data$d15N)/length(data$d15N))
#2/10: Bowman
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Bowman")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[2] <- length(data$real.name)
caddis_lookup$d13C_baseline[2] <- mean(data$d13C)
caddis_lookup$d13C_base_se[2] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[2] <- mean(data$d15N)
caddis_lookup$d15N_base_se[2] <- sqrt(var(data$d15N)/length(data$d15N))
#3/10: Kintla
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Kintla")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[3] <- length(data$real.name)
caddis_lookup$d13C_baseline[3] <- mean(data$d13C)
caddis_lookup$d13C_base_se[3] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[3] <- mean(data$d15N)
caddis_lookup$d15N_base_se[3] <- sqrt(var(data$d15N)/length(data$d15N))
#4/10: Trout
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Trout")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[4] <- length(data$real.name)
caddis_lookup$d13C_baseline[4] <- mean(data$d13C)
caddis_lookup$d13C_base_se[4] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[4] <- mean(data$d15N)
caddis_lookup$d15N_base_se[4] <- sqrt(var(data$d15N)/length(data$d15N))
#5/10: Logging
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Logging")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[5] <- length(data$real.name)
caddis_lookup$d13C_baseline[5] <- mean(data$d13C)
caddis_lookup$d13C_base_se[5] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[5] <- mean(data$d15N)
caddis_lookup$d15N_base_se[5] <- sqrt(var(data$d15N)/length(data$d15N))
#6/10: Quartz
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Quartz")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[6] <- length(data$real.name)
caddis_lookup$d13C_baseline[6] <- mean(data$d13C)
caddis_lookup$d13C_base_se[6] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[6] <- mean(data$d15N)
caddis_lookup$d15N_base_se[6] <- sqrt(var(data$d15N)/length(data$d15N))
#7/10: Hungry Horse
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Hungry Horse")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[7] <- length(data$real.name)
caddis_lookup$d13C_baseline[7] <- mean(data$d13C)
caddis_lookup$d13C_base_se[7] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[7] <- mean(data$d15N)
caddis_lookup$d15N_base_se[7] <- sqrt(var(data$d15N)/length(data$d15N))
#8/10: McDonald
data <- isotope_dataset
data <- subset(data,data$Lake_name=="McDonald")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[8] <- length(data$real.name)
caddis_lookup$d13C_baseline[8] <- mean(data$d13C)
caddis_lookup$d13C_base_se[8] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[8] <- mean(data$d15N)
caddis_lookup$d15N_base_se[8] <- sqrt(var(data$d15N)/length(data$d15N))
#9/10: Swan
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Swan")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[9] <- length(data$real.name)
caddis_lookup$d13C_baseline[9] <- mean(data$d13C)
caddis_lookup$d13C_base_se[9] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[9] <- mean(data$d15N)
caddis_lookup$d15N_base_se[9] <- sqrt(var(data$d15N)/length(data$d15N))
#10/10: Lindbergh
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Lindbergh")
data <- subset(data,data$real.name%in%c("Limnephilidae","Leptoceridae"))
caddis_lookup$n_base[10] <- length(data$real.name)
caddis_lookup$d13C_baseline[10] <- mean(data$d13C)
caddis_lookup$d13C_base_se[10] <- sqrt(var(data$d13C)/length(data$d13C))
caddis_lookup$d15N_baseline[10] <- mean(data$d15N)
caddis_lookup$d15N_base_se[10] <- sqrt(var(data$d15N)/length(data$d15N))

####################################################################
#Mayfly summary
####################################################################
mayfly_lookup <- data.frame(Lake_name=c("Big Salmon","Bowman","Kintla","Trout","Logging","Quartz","Hungry Horse","McDonald","Swan","Lindbergh"),
                     n_base=rep(NA,10),
                     d13C_baseline=rep(NA,10),
                     d13C_base_se=rep(NA,10),
                     d15N_baseline=rep(NA,10),
                     d15N_base_se=rep(NA,10))

isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_09102020.csv", header=TRUE)
#1/10: Big salmon
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Big Salmon")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[1] <- length(data$real.name)
mayfly_lookup$d13C_baseline[1] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[1] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[1] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[1] <- sqrt(var(data$d15N)/length(data$d15N))
#2/10: Bowman
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Bowman")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[2] <- length(data$real.name)
mayfly_lookup$d13C_baseline[2] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[2] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[2] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[2] <- sqrt(var(data$d15N)/length(data$d15N))
#3/10: Kintla
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Kintla")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[3] <- length(data$real.name)
mayfly_lookup$d13C_baseline[3] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[3] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[3] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[3] <- sqrt(var(data$d15N)/length(data$d15N))
#4/10: Trout
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Trout")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[4] <- length(data$real.name)
mayfly_lookup$d13C_baseline[4] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[4] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[4] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[4] <- sqrt(var(data$d15N)/length(data$d15N))
#5/10: Logging
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Logging")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[5] <- length(data$real.name)
mayfly_lookup$d13C_baseline[5] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[5] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[5] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[5] <- sqrt(var(data$d15N)/length(data$d15N))
#6/10: Quartz
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Quartz")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[6] <- length(data$real.name)
mayfly_lookup$d13C_baseline[6] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[6] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[6] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[6] <- sqrt(var(data$d15N)/length(data$d15N))
#7/10: Hungry Horse
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Hungry Horse")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[7] <- length(data$real.name)
mayfly_lookup$d13C_baseline[7] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[7] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[7] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[7] <- sqrt(var(data$d15N)/length(data$d15N))
#8/10: McDonald
data <- isotope_dataset
data <- subset(data,data$Lake_name=="McDonald")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[8] <- length(data$real.name)
mayfly_lookup$d13C_baseline[8] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[8] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[8] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[8] <- sqrt(var(data$d15N)/length(data$d15N))
#9/10: Swan
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Swan")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[9] <- length(data$real.name)
mayfly_lookup$d13C_baseline[9] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[9] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[9] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[9] <- sqrt(var(data$d15N)/length(data$d15N))
#10/10: Lindbergh
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Lindbergh")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
mayfly_lookup$n_base[10] <- length(data$real.name)
mayfly_lookup$d13C_baseline[10] <- mean(data$d13C)
mayfly_lookup$d13C_base_se[10] <- sqrt(var(data$d13C)/length(data$d13C))
mayfly_lookup$d15N_baseline[10] <- mean(data$d15N)
mayfly_lookup$d15N_base_se[10] <- sqrt(var(data$d15N)/length(data$d15N))

####################################################################
#Snail summary
####################################################################
snail_lookup <- data.frame(Lake_name=c("Big Salmon","Bowman","Kintla","Trout","Logging","Quartz","Hungry Horse","McDonald","Swan","Lindbergh"),
                            n_base=rep(NA,10),
                            d13C_baseline=rep(NA,10),
                            d13C_base_se=rep(NA,10),
                            d15N_baseline=rep(NA,10),
                            d15N_base_se=rep(NA,10))
#1/10: Big salmon
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Big Salmon")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[1] <- length(data$real.name)
snail_lookup$d13C_baseline[1] <- mean(data$d13C)
snail_lookup$d13C_base_se[1] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[1] <- mean(data$d15N)
snail_lookup$d15N_base_se[1] <- sqrt(var(data$d15N)/length(data$d15N))
#2/10: Bowman
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Bowman")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[2] <- length(data$real.name)
snail_lookup$d13C_baseline[2] <- mean(data$d13C)
snail_lookup$d13C_base_se[2] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[2] <- mean(data$d15N)
snail_lookup$d15N_base_se[2] <- sqrt(var(data$d15N)/length(data$d15N))
#3/10: Kintla
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Kintla")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[3] <- length(data$real.name)
snail_lookup$d13C_baseline[3] <- mean(data$d13C)
snail_lookup$d13C_base_se[3] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[3] <- mean(data$d15N)
snail_lookup$d15N_base_se[3] <- sqrt(var(data$d15N)/length(data$d15N))
#4/10: Trout
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Trout")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[4] <- length(data$real.name)
snail_lookup$d13C_baseline[4] <- mean(data$d13C)
snail_lookup$d13C_base_se[4] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[4] <- mean(data$d15N)
snail_lookup$d15N_base_se[4] <- sqrt(var(data$d15N)/length(data$d15N))
#5/10: Logging
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Logging")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[5] <- length(data$real.name)
snail_lookup$d13C_baseline[5] <- mean(data$d13C)
snail_lookup$d13C_base_se[5] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[5] <- mean(data$d15N)
snail_lookup$d15N_base_se[5] <- sqrt(var(data$d15N)/length(data$d15N))
#6/10: Quartz
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Quartz")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[6] <- length(data$real.name)
snail_lookup$d13C_baseline[6] <- mean(data$d13C)
snail_lookup$d13C_base_se[6] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[6] <- mean(data$d15N)
snail_lookup$d15N_base_se[6] <- sqrt(var(data$d15N)/length(data$d15N))
#7/10: Hungry Horse
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Hungry Horse")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[7] <- length(data$real.name)
snail_lookup$d13C_baseline[7] <- mean(data$d13C)
snail_lookup$d13C_base_se[7] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[7] <- mean(data$d15N)
snail_lookup$d15N_base_se[7] <- sqrt(var(data$d15N)/length(data$d15N))
#8/10: McDonald
data <- isotope_dataset
data <- subset(data,data$Lake_name=="McDonald")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[8] <- length(data$real.name)
snail_lookup$d13C_baseline[8] <- mean(data$d13C)
snail_lookup$d13C_base_se[8] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[8] <- mean(data$d15N)
snail_lookup$d15N_base_se[8] <- sqrt(var(data$d15N)/length(data$d15N))
#9/10: Swan
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Swan")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[9] <- length(data$real.name)
snail_lookup$d13C_baseline[9] <- mean(data$d13C)
snail_lookup$d13C_base_se[9] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[9] <- mean(data$d15N)
snail_lookup$d15N_base_se[9] <- sqrt(var(data$d15N)/length(data$d15N))
#10/10: Lindbergh
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Lindbergh")
data <- subset(data,data$real.name%in%c("Physidae","Planorbidae","Lymnaeidae"))
snail_lookup$n_base[10] <- length(data$real.name)
snail_lookup$d13C_baseline[10] <- mean(data$d13C)
snail_lookup$d13C_base_se[10] <- sqrt(var(data$d13C)/length(data$d13C))
snail_lookup$d15N_baseline[10] <- mean(data$d15N)
snail_lookup$d15N_base_se[10] <- sqrt(var(data$d15N)/length(data$d15N))

snail_lookup <- subset(snail_lookup,Lake_name!="Swan")#No snails were collected in Swan Lake

```


Next, we produced a table of $\delta$^13^C and $\delta$^15^N mean standard error (se) for each trophic baseline candidate in each lake, to determine which candidate had the lowest among-lake isotopic variability. This summary is provided below (Table 3) and in the article's SI Table 5.

```{r, echo=TRUE}
####################################################################
#make comparison table, SI Table 5
####################################################################

compare <- data.frame(Candidate=c("Periphyton","Mayfly","Zooplankton","Caddisfly","Hydrachnidia","Snail"),
                      n=rep(NA,6),
                      mean_d13C_base_se=rep(NA,6),
                      mean_d15N_base_se=rep(NA,6))
#n
compare$n[1] <- sum(periphyton_isotope_summary$n)
compare$n[2] <- sum(mayfly_lookup$n_base)
compare$n[3] <- sum(zooplankton_isotope_summary$n)
compare$n[4] <- sum(caddis_lookup$n_base)
compare$n[5] <- sum(hydrachnidia_isotope_summary$n)
compare$n[6] <- sum(snail_lookup$n)
#d13C_se
compare$mean_d13C_base_se[1] <- mean(periphyton_isotope_summary$se_d13C)
compare$mean_d13C_base_se[2] <- mean(mayfly_lookup$d13C_base_se)
compare$mean_d13C_base_se[3] <- mean(zooplankton_isotope_summary$se_d13C)
compare$mean_d13C_base_se[4] <- mean(caddis_lookup$d13C_base_se)
compare$mean_d13C_base_se[5] <- mean(hydrachnidia_isotope_summary$se_d13C)
compare$mean_d13C_base_se[6] <- mean(snail_lookup$d13C_base_se)
compare$mean_d13C_base_se <- round(compare$mean_d13C_base_se,2)
#d15N_se
compare$mean_d15N_base_se[1] <- mean(periphyton_isotope_summary$se_d15N)
compare$mean_d15N_base_se[2] <- mean(mayfly_lookup$d15N_base_se)
compare$mean_d15N_base_se[3] <- mean(zooplankton_isotope_summary$se_d15N)
compare$mean_d15N_base_se[4] <- mean(caddis_lookup$d15N_base_se)
compare$mean_d15N_base_se[5] <- mean(hydrachnidia_isotope_summary$se_d15N)
compare$mean_d15N_base_se[6] <- mean(snail_lookup$d15N_base_se)
compare$mean_d15N_base_se <- round(compare$mean_d15N_base_se,2)
#re-name columns
colnames(compare)[3] <- "Mean $\\delta$^15^N SE ($\\unicode{x2030}$)"
colnames(compare)[4] <- "Mean $\\delta$^13^C SE ($\\unicode{x2030}$)"
#make table
library(knitr)
kable(compare, align="lccc", caption = "Table 3. Among-lake isotopic variability for all isotope baseline candidates.")
```


For subsequent isotope analyses, we used mayflies (Ephemeroptera: Baetidae, Leptophlebiidae, Heptageniidae, and Ephemerellidae) as our lake-specific trophic baselines because mayflies had the most consistent isotope signature of any isotope baseline candidate (Table 3, lowest standard error; Article's SI Table 5). Mayflies are a common choice for isotope baseline in rivers (e.g., Jardine et. al., 2014; Rogosch and Olden, 2020) and can be used as $\delta$^13^C and $\delta$^15^N baselines in lakes (Vander Zanden and Rasmussen, 1999; de Eyto et. al., 2020). Therefore, mayflies were the best available choice for trophic baseline in this study.

###SI Table 2, lake summary table

With mayflies selected as the trophic baseline, we added lake-specific mayfly $\delta$^13^C and $\delta$^15^N baselines to a table summarizing lake characteristics. The final result was Table 4 (Article's SI Table 2).

Lake summary table:
```{r, echo=TRUE}
#########################################################################################################
#Lake table
#########################################################################################################
#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)

#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Lake_name"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C", "d15N", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         Lake_surface_elev_m = mean(Lake_surface_elev_m),
                         Surface_area_ha = mean(Surface_area_ha),
                         Lake_max_depth_m = mean(Lake_max_depth_m),
                         Secchi_m = mean(Secchi_m),
                         Current_conversion = mean(Current_conversion)
                         
  )
  
  if (printSummary)  print(summary)
  
  summary
  return(summary)
}

#create summary dataframe using function and grouping specified above
isotope_summary_df<-summariseIsotopeData(df = subset(isotope_dataset),
                                         grouping = c("Lake_name"))
head(isotope_summary_df)

##NOT RUN but useful
#use binomial linear regression coefficients to calculate modelled "timestep" for each lake
#isotope_summary_df$prop <- isotope_summary_df$Current_conversion/(1-isotope_summary_df$Current_conversion) #for use in logit link function
#isotope_summary_df$log <- log10(isotope_summary_df$prop) #logit link function
#isotope_summary_df$Time <- (isotope_summary_df$log+1.97331)/.09127 #binomial linear regression coefficients (calculated in Figure 7, below)


```

Isotopic baseline table:
```{r, echo=TRUE}
#########################################################################################################
#Isotopic baselines for lake table
#########################################################################################################
lookup <- data.frame(Lake_name=c("Big  Salmon","Bowman","Kintla","Trout","Logging","Quartz","Hungry Horse","McDonald","Swan","Lindbergh"),
                     n_base=rep(NA,10),
                     d13C_baseline=rep(NA,10),
                     d13C_base_se=rep(NA,10),
                     d15N_baseline=rep(NA,10),
                     d15N_base_se=rep(NA,10)) #empty dataframe to receive calculated values

isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
#1/10: Big salmon
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Big Salmon")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[1] <- length(data$real.name)
lookup$d13C_baseline[1] <- mean(data$d13C)
lookup$d13C_base_se[1] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[1] <- mean(data$d15N)
lookup$d15N_base_se[1] <- sqrt(var(data$d15N)/length(data$d15N))
#2/10: Bowman
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Bowman")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[2] <- length(data$real.name)
lookup$d13C_baseline[2] <- mean(data$d13C)
lookup$d13C_base_se[2] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[2] <- mean(data$d15N)
lookup$d15N_base_se[2] <- sqrt(var(data$d15N)/length(data$d15N))
#3/10: Kintla
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Kintla")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[3] <- length(data$real.name)
lookup$d13C_baseline[3] <- mean(data$d13C)
lookup$d13C_base_se[3] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[3] <- mean(data$d15N)
lookup$d15N_base_se[3] <- sqrt(var(data$d15N)/length(data$d15N))
#4/10: Trout
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Trout")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[4] <- length(data$real.name)
lookup$d13C_baseline[4] <- mean(data$d13C)
lookup$d13C_base_se[4] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[4] <- mean(data$d15N)
lookup$d15N_base_se[4] <- sqrt(var(data$d15N)/length(data$d15N))
#5/10: Logging
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Logging")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[5] <- length(data$real.name)
lookup$d13C_baseline[5] <- mean(data$d13C)
lookup$d13C_base_se[5] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[5] <- mean(data$d15N)
lookup$d15N_base_se[5] <- sqrt(var(data$d15N)/length(data$d15N))
#6/10: Quartz
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Quartz")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[6] <- length(data$real.name)
lookup$d13C_baseline[6] <- mean(data$d13C)
lookup$d13C_base_se[6] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[6] <- mean(data$d15N)
lookup$d15N_base_se[6] <- sqrt(var(data$d15N)/length(data$d15N))
#7/10: Hungry Horse
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Hungry Horse")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[7] <- length(data$real.name)
lookup$d13C_baseline[7] <- mean(data$d13C)
lookup$d13C_base_se[7] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[7] <- mean(data$d15N)
lookup$d15N_base_se[7] <- sqrt(var(data$d15N)/length(data$d15N))
#8/10: McDonald
data <- isotope_dataset
data <- subset(data,data$Lake_name=="McDonald")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[8] <- length(data$real.name)
lookup$d13C_baseline[8] <- mean(data$d13C)
lookup$d13C_base_se[8] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[8] <- mean(data$d15N)
lookup$d15N_base_se[8] <- sqrt(var(data$d15N)/length(data$d15N))
#9/10: Swan
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Swan")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[9] <- length(data$real.name)
lookup$d13C_baseline[9] <- mean(data$d13C)
lookup$d13C_base_se[9] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[9] <- mean(data$d15N)
lookup$d15N_base_se[9] <- sqrt(var(data$d15N)/length(data$d15N))
#10/10: Lindbergh
data <- isotope_dataset
data <- subset(data,data$Lake_name=="Lindbergh")
data <- subset(data,data$real.name%in%c("Baetidae","Leptophlebiidae","Heptageniidae","Ephemerellidae"))
lookup$n_base[10] <- length(data$real.name)
lookup$d13C_baseline[10] <- mean(data$d13C)
lookup$d13C_base_se[10] <- sqrt(var(data$d13C)/length(data$d13C))
lookup$d15N_baseline[10] <- mean(data$d15N)
lookup$d15N_base_se[10] <- sqrt(var(data$d15N)/length(data$d15N))

head(lookup) #lake-specific isotopic baselines
```


```{r, echo=TRUE}
#########################################################################################################
#Combine lake table and isotopic baseline tables. SI Table 2.
#########################################################################################################
SI_table_2 <- merge(isotope_summary_df,lookup,by.x="Lake_name",by.y="Lake_name") #combine
SI_table_2$d13C_baseline <- round(SI_table_2$d13C_baseline,2) #round
SI_table_2$d15N_baseline <- round(SI_table_2$d15N_baseline,2) #round
SI_table_2 <- subset(SI_table_2,select=-c(n_base,d13C_base_se,d15N_base_se))
#re-name columns
colnames(SI_table_2)[1] <- "Lake name"
colnames(SI_table_2)[2] <- "Surface elev (m)"
colnames(SI_table_2)[3] <- "Surface area (ha)"
colnames(SI_table_2)[4] <- "Max depth (m)"
colnames(SI_table_2)[5] <- "Secchi (m)"
colnames(SI_table_2)[6] <- "Current conversion"
colnames(SI_table_2)[7] <- "Baseline <br> $\\delta$^13^C ($\\unicode{x2030}$)"
colnames(SI_table_2)[8] <- "Baseline <br> $\\delta$^15^N ($\\unicode{x2030}$)"
#make table
library(knitr)
kable(SI_table_2, align="lcccccc", caption = "Table 4. Lake summary table with isotopic baselines added (Article's SI Table 2).")
```

###SI Table 1, summary of baseline-corrected fish $\delta$^13^C and $\delta$^15^N

Next, we corrected each sample's empirical $\delta$^13^C and $\delta$^15^N using ${\delta}X_{corrected}$ = ${\delta}X_{empirical}$ - ${\delta}X_{baseline}$ described in Hobson et. al. (2012). Then we produced Table 5 (Article's SI Table 1) as a summary of our fish dataset.

```{r, echo=FALSE}

rm(list=ls())
graphics.off()
set.seed(1)
```
```{r, echo=TRUE, results=TRUE}
#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
isotope_dataset$d13C_corr <- isotope_dataset$d13C - isotope_dataset$d13C_baseline #calculate lake-specific baseline-corrected d13C
isotope_dataset$d15N_corr <- isotope_dataset$d15N - isotope_dataset$d15N_baseline #calculate lake-specific baseline-corrected d15N
```

```{r, echo=FALSE}
rm(list=ls())
graphics.off()
set.seed(1)
#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")
```

```{r, echo=TRUE, results=TRUE}
#########################################################################################################
#Fish summary table (SI Table 1)
#########################################################################################################
#write summary function
#summary function based on: https://github.com/clquezada/tRophicPosition/blob/master/R/summariseIsotopeData.R
summariseIsotopeData <- function (df = NULL, grouping = c("Invasion_category", "Functional_group"),
                                  printSummary = FALSE, ...){
  
  d13C <- d15N <- NULL
  
  if (is.null(checkmate::checkNames(df, c("d13C_corr", "d15N_corr", grouping))))
    stop("Check the grouping variable or the names in your dataframe")
  
  summary <- plyr::ddply(df, grouping, plyr::summarise,
                         n = length(d13C_corr),
                         mean_d13C = mean(d13C_corr),
                         sd_d13C = sd(d13C_corr),
                         se_d13C = sqrt(var(d13C_corr)/length(d13C_corr)),
                         mean_d15N = mean(d15N_corr),
                         sd_d15N = sd(d15N_corr),
                         se_d15N = sqrt(var(d15N_corr)/length(d15N_corr)),
                         meanL_mm = mean(Total_length_mm),
                         seL_mm = sqrt(var(Total_length_mm)/length(Total_length_mm))
                         
  )
  
  if (printSummary)  print(summary)
  
  summary
  return(summary)
}

#create summary dataframe using function and grouping specified above
isotope_summary_df<-summariseIsotopeData(df = subset(data),
                                         grouping = c("Invasion_category", "Functional_group"))
#round
isotope_summary_df$mean_d13C <- round(isotope_summary_df$mean_d13C,2)
isotope_summary_df$sd_d13C <- round(isotope_summary_df$sd_d13C,2)
isotope_summary_df$se_d13C <- round(isotope_summary_df$se_d13C,2)
isotope_summary_df$mean_d15N <- round(isotope_summary_df$mean_d15N,2)
isotope_summary_df$sd_d15N <- round(isotope_summary_df$sd_d15N,2)
isotope_summary_df$se_d15N <- round(isotope_summary_df$se_d15N,2)
isotope_summary_df$meanL_mm <- round(isotope_summary_df$meanL_mm,0)
isotope_summary_df$seL_mm <- round(isotope_summary_df$seL_mm,0)
isotope_summary_df <- subset(isotope_summary_df,select=-c(sd_d13C,sd_d15N))
#re-name columns
colnames(isotope_summary_df)[1] <- "Invasion<br>category"
colnames(isotope_summary_df)[2] <- "Functional<br>group"
colnames(isotope_summary_df)[4] <- "Mean <br> $\\delta$^13^C ($\\unicode{x2030}$)"
colnames(isotope_summary_df)[5] <- "SE <br> $\\delta$^13^C ($\\unicode{x2030}$)"
colnames(isotope_summary_df)[6] <- "Mean <br> $\\delta$^15^N ($\\unicode{x2030}$)"
colnames(isotope_summary_df)[7] <- "SE <br> $\\delta$^15^N ($\\unicode{x2030}$)"
colnames(isotope_summary_df)[8] <- "Mean total<br>length (mm)"
colnames(isotope_summary_df)[9] <- "SE total <br> length (mm)"
#make table
library(knitr)
kable(subset(isotope_summary_df), align="llccccccc", caption = "Table 5. Summary of fish isotope data (Article SI Table 1).")
```

###Figure 1, baseline-corrected fish $\delta$^13^C and $\delta$^15^N
Figure 1 is a scatterplot of fish baseline-corrected $\delta$^13^C and $\delta$^15^N in each invasion category. Ellipses are 95% confidence intervals around mean $\delta$^13^C and $\delta$^15^N for each fish group.
```{r, echo=TRUE, warning=FALSE}
#########################################################################################################
#Figure 1: facet wrap plot of d13C and d15N with 95% confidence interval ellipses 
#########################################################################################################
rm(list=ls())
graphics.off()
set.seed(1)
library(ggplot2)

#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")

#Put variables in correct order
data$Invasion_category <- factor(data$Invasion_category, levels= c("Reference","Mid","Late"))
data$Functional_group <- factor(data$Functional_group, levels= c("Bull trout","Lake trout","Littoral forage fish","Generalist fish","Pelagic forage fish"))
head(subset(data, select=c(d13C_corr,d15N_corr,Functional_group,Invasion_category)))
```

```{r, echo=TRUE, fig.width=12}
#Plot Figure 1

ggplot(data,aes(x=d13C_corr, y=d15N_corr, shape=Functional_group))+
  facet_wrap(~Invasion_category,ncol=5)+
  geom_point(size = 3)+
  stat_ellipse(aes(x=d13C_corr, y=d15N_corr, fill=Functional_group),
               geom="polygon",
               alpha = 0.4, #transparency
               level = 0.95)+ #confidence interval
  theme_bw()+
  labs(x=expression(~delta^13*'C'~'(\211)'),
       y=expression(~delta^15*'N'~'(\211)'),
       fill = "Species",
       shape = "Species",
       group = "Species")+
  theme(axis.text = element_text(size = 20, color = "black"),
        axis.title.x = element_text(size = 20, vjust=-2),
        axis.title.y = element_text(size = 20, vjust=2),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 20),
        strip.text = element_text(size = 20))
```

###Figure 2, SEA.b, ellipse area in SIBER

This figure uses R package "SIBER" to calculate a posterior distribution of estimates of the isotopic dietary breadth population parameter.

```{r, echo = TRUE, fig.width=14, fig.height=6, results=TRUE, warning=FALSE, message=FALSE}

#https://cran.r-project.org/web/packages/SIBER/vignettes/Introduction-to-SIBER.html
#https://cran.r-project.org/web/packages/SIBER/vignettes/Ellipse-Overlap.html
#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")

#prep data to meet siber specs
siber_df <- data[c("d13C_corr", "d15N_corr", "Functional_group", "Invasion_category")] #select siber's columns from full dataset
names(siber_df) <- c("iso1","iso2","group","community") #re-name columns to match siber
head(siber_df)
```
```{r, echo = TRUE, fig.width=14, fig.height=6, results=FALSE, warning=FALSE, message=FALSE}
siber_object <- createSiberObject(siber_df) #create siber object (the input for the mcmc)
group.ML <- groupMetricsML(siber_object) #calculate TA, SEA, SEA.c

# define siber mcmc model parameters
parms <- list()
parms$n.iter <- 2 * 10^4   # number of iterations to run the model
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10     # thin the posterior by this many
parms$n.chains <- 2        # run this many chains

# define siber mcmc model priors
priors <- list()
priors$R <- 1 * diag(2) #uninformative
priors$k <- 2 #uninformative
priors$tau.mu <- 1.0E-3 #uninformative

#run siber mcmc
ellipses.posterior <- siberMVN(siber_object, parms, priors) #run the model


#calculate TA, SEA, SEAc for each group
group.ML <- groupMetricsML(siber_object) #ellipse area, NOT ellipse coordinates
#draw posteriors for siber ellipses for each group.community combination
SEA.B <- siberEllipses(ellipses.posterior) #motherlode
#SEA.B credible intervals:
cr.p <- c(0.95, 0.99) # vector of quantiles
#calculate credible intervals of SEA.B
SEA.B.credibles <- lapply(# call to hdrcde:hdr using lapply()
  as.data.frame(SEA.B), 
  function(x,...){tmp<-hdrcde::hdr(x)$hdr},
  prob = cr.p) #vectors of the credible intervals for each group's sea.b

#siber has a built-in SEA.B plot: https://cran.r-project.org/web/packages/SIBER/vignettes/Introduction-to-SIBER.html
#but we'll re-do siberDensityPlot in ggplot

#draw posteriors for siber ellipses area for each group.community combination
SEA.B_2 <- SEA.B # copy siber's output list
colnames(SEA.B_2) <- colnames(group.ML) #re-name columns in new list
SEA.B_2 <- as.data.frame(SEA.B_2) #make dataframe from list
SEA.B_2 <- melt(SEA.B_2) #re-arrange variables in list
SEA.B_2$combo <- SEA.B_2$variable #re-name $combo (from melt()) into $variable, as required for separate()
SEA.B_2 <- separate(data = SEA.B_2,
                    col = variable,
                    into = c("invasion.category","Species"),
                    sep = "\\.",
                    extra = "merge") #use separate() to produce two useful columns
SEA.B_2$invasion.category <- factor(SEA.B_2$invasion.category,
                                    levels=c('Reference',
                                             'Mid',
                                             'Late')) #re-order variables

SEA.B_2$Species <- factor(SEA.B_2$Species,levels=c("Bull trout","Lake trout","Littoral forage fish","Generalist fish","Pelagic forage fish")) #re-order variables
```
```{r, echo = TRUE, fig.width=14, fig.height=6, results=TRUE, warning=FALSE, message=TRUE}
head(SEA.B_2)
#########################################################################################################
#Figure 2: SEA.b, ellipse area in SIBER 
#########################################################################################################
library(ggplot2)
ggplot(data = SEA.B_2,
       aes(x=invasion.category, y=value, fill=invasion.category)) +
  geom_boxplot()+
  scale_fill_manual(values=c("white","light gray","grey32"))+
  theme_bw()+
  scale_y_continuous(limits=c(0,25))+
  labs(y = expression("Isotopic dietary breadth ( \211"^2~")"),
       x = "\nInvasion")+
  facet_wrap(~Species,
             ncol = 5)+
  theme(axis.text.x = element_text(size=20, angle = 90, color = "black", vjust=0.5),
        axis.text.y = element_text(size=20, color = "black"),
        axis.title.x = element_text(size=20),
        axis.title.y = element_text(size=20, vjust=3),
        legend.position = "none",
        strip.text = element_text(size=20, color = "black"))

```

###Figure 3, mixed effects linear modelling for $\delta$^15^N
The mixed effects models need to be done in two parts: 1) all fish other than lake trout 2) lake trout. There are no lake trout in reference lakes so we must run the models separately and then combine.

First, run the model for all fish besides lake trout.

```{r, echo=TRUE, results=TRUE}
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")

#d15N model first part
data$Invasion_category <- factor(data$Invasion_category, levels= c("Reference","Mid","Late")) #re-order variables for plotting
data <- subset(data,real.name!="Lake trout") #must exclude lake trout because lake trout are not present in all invasion categories
head(subset(data, select=c(Lake_name,Invasion_category,Functional_group,d15N_corr)))
```
```{r, echo=TRUE, results=FALSE}
d15N_lmer <- lmer(d15N_corr ~ Invasion_category*Functional_group+ #predict d15N_corr for each unique combination of invasion category and functional group
                   (1|Lake_name), #lake is random effect
                 data=data)
#summary(d15N_lmer)
invasion <- predictorEffect("Invasion_category", d15N_lmer)#compute lmer results for plotting #https://www.rdocumentation.org/packages/effects/versions/4.2-0/topics/predictorEffects
invasion <- as.data.frame(invasion) #make results a dataframe for ggplot
```

Then run the second part of the model to generate lake trout results.

```{r, echo=TRUE, results=TRUE}
#d15N model second part
data <- subset(isotope_dataset,Functional_group%in%c("Lake trout","Bull trout")) #now we include lake trout
data <- subset(data,Invasion_category != "Reference") #there are no lake trout in reference lakes
data$Invasion_category <- factor(data$Invasion_category, levels = c("Mid","Late")) #re-order variables
head(subset(data, select=c(Lake_name,Invasion_category,Functional_group,d15N_corr)))
```
```{r, echo=TRUE, results=FALSE}
d15N_lmer <- lmer(d15N_corr ~ Invasion_category*Functional_group+ #predict d15N_corr for each unique combination of invasion category and functional group
                   (1|Lake_name), #lake is random effect
                 data=data)

invasion2 <- predictorEffect("Invasion_category", d15N_lmer)#compute lmer results for plotting #https://www.rdocumentation.org/packages/effects/versions/4.2-0/topics/predictorEffects
invasion2 <- as.data.frame(invasion2) #make results a dataframe for ggplot
invasion2 <- subset(invasion2,Functional_group!="Bull trout")#exclude bull trout from invasion2 so we aren't duplicating
```

Finally, combine all results and plot.

```{r, echo=TRUE, fig.width=14, fig.height=6}
#combine parts 1 and 2
test <- rbind(invasion,invasion2)#combine
test$Functional_group <- factor(test$Functional_group, levels= c("Bull trout","Lake trout","Littoral forage fish","Generalist fish","Pelagic forage fish"))#re-order variables
test$Invasion_category <- factor(test$Invasion_category, levels= c("Reference","Mid","Late"))#re-order variables
head(test)
#Plot figure 3
ggplot(test, aes (x= Invasion_category, y=fit, group=Functional_group))+
  geom_line()+
  geom_point(size = 4) +
  geom_errorbar(width=.1, aes(ymin=fit-se, ymax=fit+se)) +
  labs(y= expression(paste(delta^15, "N"~'(\211)')), x="\nInvasion")+
  facet_wrap(~Functional_group, ncol=5)+
  theme_bw()+
  theme(axis.text.x = element_text(size = 20, angle=90, color = "black", vjust=0.5),
        axis.text.y = element_text(size = 20, color = "black"),
        axis.title.y = element_text(size = 20, vjust = 4),
        axis.title.x = element_text(size = 20),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        strip.text = element_text(size = 20))

```

To generate p-values, cycle through each functional group, like this:

```{r, echo=TRUE, results=TRUE}
#p-values for bull trout d15N
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")
data$Invasion_category <- factor(data$Invasion_category, levels= c("Reference","Mid","Late"))#re-order variables for plotting
data$Functional_group <- factor(data$Functional_group, levels= c("Bull trout","Generalist fish","Pelagic forage fish","Lake trout","Littoral forage fish")) #Produce p-value for the first factor in this list. Change the order of levels to generate each p-value.
data <- subset(data,real.name!="Lake trout")
d15N_lmer <- lmer(d15N_corr ~ Invasion_category*Functional_group+ #predict d15N_corr for each unique combination of invasion category and functional group
                   (1|Lake_name), #lake is random effect
                 data=data)
#summary(d15N_lmer) #bull trout d15N p-values comparing against Reference are from the model summary. Specifically, the Pr(>|t|) column and rows Invasion_categoryMid and Invasion_categoryLate.
```

````{r, echo=TRUE}
summary(d15N_lmer)[["coefficients"]][50] #p-value for Mid-invasion bull trout d15N compared to reference bull trout d15N
summary(d15N_lmer)[["coefficients"]][51] #p-value for Late-invasion bull trout d15N compared to reference bull trout d15N
```

###Figure 4, mixed effects linear modelling for $\delta$^13^C

Like the $\delta$^15^Nd15N model, the $\delta$^13^C mixed effects models need to be done in two parts: 1) all fish other than lake trout 2) lake trout. There are no lake trout in reference lakes so we must run the models separately and then combine.

First, run the model for all fish besides lake trout.

```{r, echo=TRUE, results=TRUE}
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")

#d13C model first part
data$Invasion_category <- factor(data$Invasion_category, levels= c("Reference","Mid","Late"))#re-order variables for plotting
data <- subset(data,real.name!="Lake trout") #must exclude lake trout because lake trout are not present in all invasion categories
head(subset(data, select=c(Lake_name,Invasion_category,Functional_group,d13C_corr)))
```
```{r, echo=TRUE, results=FALSE}
d13C_lmer <- lmer(d13C_corr ~ Invasion_category*Functional_group+ #predict d13Ccorr for each unique combination of invasion category and functional group
                   (1|Lake_name), #lake is random effect
                 data=data)
#summary(d13C_lmer)
invasion <- predictorEffect("Invasion_category", d13C_lmer)#compute lmer results for plotting #https://www.rdocumentation.org/packages/effects/versions/4.2-0/topics/predictorEffects
invasion <- as.data.frame(invasion) #make results a dataframe for ggplot
```

Then run the second part of the model to generate lake trout results.

```{r, echo=TRUE, results=TRUE}
#d13C model second part
data <- subset(isotope_dataset,Functional_group%in%c("Lake trout","Bull trout")) #now we include lake trout
data <- subset(data,Invasion_category != "Reference") #there are no lake trout in reference lakes
data$Invasion_category <- factor(data$Invasion_category, levels = c("Mid","Late")) #re-order variables for plotting
head(subset(data, select=c(Lake_name,Invasion_category,Functional_group,d13C_corr)))
```
```{r, echo=TRUE, results=FALSE}
d13C_lmer <- lmer(d13C_corr ~ Invasion_category*Functional_group+ #predict d13Ccorr for each unique combination of invasion category and functional group
                   (1|Lake_name), #lake is random effect
                 data=data)
#summary(d13C_lmer)#p-values
invasion2 <- predictorEffect("Invasion_category", d13C_lmer) #compute lmer results for plotting #https://www.rdocumentation.org/packages/effects/versions/4.2-0/topics/predictorEffects
invasion2 <-as.data.frame(invasion2) #make results a dataframe for ggplot
invasion2 <- subset(invasion2,Functional_group!="Bull trout")#exclude bull trout from invasion2 so we aren't duplicating
```

```{r, echo=TRUE, fig.width=14, fig.height=6}
#combine parts 1 and 2
test <- rbind(invasion,invasion2)
test$Functional_group <- factor(test$Functional_group, levels= c("Bull trout","Lake trout","Littoral forage fish","Generalist fish","Pelagic forage fish")) #re-order variables for plotting
test$Invasion_category <- factor(test$Invasion_category, levels= c("Reference","Mid","Late")) #re-order variables for plotting
head(test)
#Plot figure 4
ggplot(test, aes (x= Invasion_category, y=fit, group=Functional_group))+
  geom_line()+
  geom_point(size = 4) +
  geom_errorbar(width=.1, aes(ymin=fit-se, ymax=fit+se)) +
  labs(y= expression(paste(delta^13, "C"~'(\211)')), x="\nInvasion")+
  facet_wrap(~Functional_group, ncol=5)+
  theme_bw()+
  theme(axis.text.x = element_text(size = 20, angle=90, color = "black", vjust=0.5),
        axis.text.y = element_text(size = 20, color = "black"),
        axis.title.y = element_text(size = 20, vjust = 4),
        axis.title.x = element_text(size = 20),
        legend.text = element_text(size = 15),
        legend.title = element_text(size = 15),
        strip.text = element_text(size = 20))
```

To generate p-values, cycle through each functional group, like this:

```{r, echo=TRUE, results=TRUE}
#p-values for bull trout
rm(list=ls())
graphics.off()
set.seed(1)
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")
data$Invasion_category <- factor(data$Invasion_category, levels= c("Reference","Mid","Late"))
data$Functional_group <- factor(data$Functional_group, levels= c("Bull trout","Generalist fish","Pelagic forage fish","Lake trout","Littoral forage fish"))#Produce p-value for the first factor in this list. Change the order of levels to generate each p-value.
data <- subset(data,real.name!="Lake trout")
d13C_lmer <- lmer(d13C_corr ~ Invasion_category*Functional_group+
                   (1|Lake_name),
                 data=data)
#summary(d13C_lmer) #bull trout d13C p-values comparing against Reference are from the model summary. Specifically, the Pr(>|t|) column and rows Invasion_categoryMid and Invasion_categoryLate.
```

```{r echo=TRUE}
summary(d13C_lmer)[["coefficients"]][50] #p-value for Mid-invasion bull trout d13C compared to reference bull trout d13C
summary(d13C_lmer)[["coefficients"]][51] #p-value for Late-invasion bull trout d13C compared to reference bull trout d13C
```

###Figure 5, isotope niche overlap 

This analysis uses R package "nicheROVER" to estimate diet overlap between bull trout and lake trout or vice versa in mid- and late-invasion lakes. Run the model separately for mid-invasion lakes and late-invasion lakes.

First, run the model for late-invasion lakes.

```{r, echo = TRUE, results = TRUE}
rm(list=ls())
graphics.off()
set.seed(1)

#load data
isotope_dataset <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/isotope_dataset_MT_lake_trout_food_web_11072020.csv", header=TRUE)
data <- subset(isotope_dataset,Sample_type=="Fish")

#STEP 1: RUN THIS SECTION FOR Invasion_category == Late

#prep data
nicheROVER_df <- subset(data, Invasion_category %in% "Late" & real.name %in% c("Bull trout","Lake trout")) #select Late blt and lkt
#eliminate all extra columns (nicheROVER requires only 3 columns as inputs)
nicheROVER_df <- nicheROVER_df[c("Functional_group","d13C_corr", "d15N_corr")] #columns must be in this order (species, d13C, d15N)
names(nicheROVER_df) <- c("species","D13C","D15N") #re-name columns to match nicherover
head(nicheROVER_df)
```
```{r, echo = TRUE, results = FALSE}
#aggregate data as required for nicheROVER
aggregate(nicheROVER_df[2:3], nicheROVER_df[1], mean) #point estimates of mean isotope values

#create fish.par object
nsamples <- 1000
system.time({
  fish.par <- tapply(1:nrow(nicheROVER_df), nicheROVER_df$species, function(ii) niw.post(nsamples = nsamples, 
                                                                                         X = nicheROVER_df[ii, 2:3]))
})
#remove NULLs from fish.par list
#fish.par should be a list of length = number of species (or whatever factor you're aggregating by)
#https://stackoverflow.com/questions/33004238/r-removing-null-elements-from-a-list
fish.par <- fish.par[-which(sapply(fish.par, is.null))] #get rid of NULLs from our fish.par list


#estimate niche overlaps
nsamples <- 1000

# Overlap calculation.  use nsamples = nprob = 10000 (1e4) for higher
# accuracy.  the variable over.stat can be supplied directly to the
# overlap.plot function

over.stat <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = c(0.95,.99))

# The mean overlap metrics calculated across iteratations for both niche
# region sizes (alpha = .95 and alpha = .99) can be calculated and displayed
# in an array.
over.mean <- apply(over.stat,
                   c(1:2, 4),
                   mean) * 100
round(over.mean, 2)

over.cred <- apply(over.stat * 100,
                   c(1:2, 4),
                   quantile,
                   prob = c(0.025, 0.975), 
                   na.rm = TRUE)
round(over.cred[, , , 1])  # display alpha = .95 niche region

#generate array of overlap probabilities for species A vs species B using over.mean and and over.cred
over.stat <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = 0.95)

#re-do the package's plotting function to prep data for ggplot
over.hist <- apply(over.stat, 1:2, function(x) {
  if (any(is.na(x))) 
    return(NULL)
  else {
    tmp <- x * 100
  }
  tmp
})
```
```{r, echo = TRUE, results = TRUE}
#How much do bull trout overlap with lake trout?
data1<-as.data.frame(over.hist[[1,2]])#make dataframe from over.hist (the model's output)
colnames(data1)<-"Overlap"

#Add column to data1, add labels
BLT_Late <- data1
BLT_Late$Species <- "Bull trout"
BLT_Late$Invasion <- "Late"

head(BLT_Late)

#How much do lake trout overlap bull trout?
data1<-as.data.frame(over.hist[[2,1]])#make dataframe from over.hist (the model's output)
colnames(data1)<-"Overlap"

#Add column to data1, label rows
LKT_Late <- data1
LKT_Late$Species <- "Lake trout"
LKT_Late$Invasion <- "Late"

head(LKT_Late)
```

Next, run the same model for mid-invasion lakes.
```{r, echo = TRUE, results = FALSE}
#STEP 2: RUN THIS SECTION FOR Invasion_category == Mid

#prep data
nicheROVER_df <- subset(data, Invasion_category %in% "Mid" & real.name %in% c("Bull trout","Lake trout")) #select Mid blt and lkt #select blt and lkt
#eliminate all extra columns (nicheROVER requires only 3 columns as inputs)
nicheROVER_df <- nicheROVER_df[c("Functional_group","d13C_corr", "d15N_corr")] #columns must be in this order (species, d13C, d15N)
names(nicheROVER_df) <- c("species","D13C","D15N") #re-name columns to match nicherover

#aggregate data as required for nicheROVER
aggregate(nicheROVER_df[2:3], nicheROVER_df[1], mean) #point estimates of mean isotope values

#create fish.par object
nsamples <- 1000
system.time({
  fish.par <- tapply(1:nrow(nicheROVER_df), nicheROVER_df$species, function(ii) niw.post(nsamples = nsamples, 
                                                                                         X = nicheROVER_df[ii, 2:3]))
})
#remove NULLs from fish.par list
#fish.par should be a list of length = number of species (or whatever factor you're aggregating by)
#https://stackoverflow.com/questions/33004238/r-removing-null-elements-from-a-list
fish.par <- fish.par[-which(sapply(fish.par, is.null))] #get rid of NULLs from our fish.par list


#estimate niche overlaps
nsamples <- 1000

# Overlap calculation.  use nsamples = nprob = 10000 (1e4) for higher
# accuracy.  the variable over.stat can be supplied directly to the
# overlap.plot function

over.stat <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = c(0.95,.99))

# The mean overlap metrics calculated across iteratations for both niche
# region sizes (alpha = .95 and alpha = .99) can be calculated and displayed
# in an array.
over.mean <- apply(over.stat,
                   c(1:2, 4),
                   mean) * 100
round(over.mean, 2)

over.cred <- apply(over.stat * 100,
                   c(1:2, 4),
                   quantile,
                   prob = c(0.025, 0.975), 
                   na.rm = TRUE)
round(over.cred[, , , 1])  # display alpha = .95 niche region

#generate array of overlap probabilities for species A vs species B using over.mean and and over.cred
over.stat <- overlap(fish.par, nreps = nsamples, nprob = 1000, alpha = 0.95)

#re-do the package's plotting function to prep data for ggplot
over.hist <- apply(over.stat, 1:2, function(x) {
  if (any(is.na(x))) 
    return(NULL)
  else {
    tmp <- x * 100
  }
  tmp
})
```
```{r, echo = TRUE, results = TRUE}

#How much do bull trout overlap with lake trout?
data1<-as.data.frame(over.hist[[1,2]])#make dataframe from over.hist (the model's output)
colnames(data1)<-"Overlap"

#Add column to data1, label rows
BLT_Mid <- data1
BLT_Mid$Species <- "Bull trout"
BLT_Mid$Invasion <- "Mid"
head(BLT_Mid)
#How much do lake trout overlap bull trout?
data1<-as.data.frame(over.hist[[2,1]])#make dataframe from over.hist (the model's output)
colnames(data1)<-"Overlap"

#Add column to data1, label rows
LKT_Mid <- data1
LKT_Mid$Species <- "Lake trout"
LKT_Mid$Invasion <- "Mid"
head(LKT_Mid)
```

Then, combine results for mid-invasion and late-invasion lake results:

```{r, echo = TRUE}
BLT_both <- rbind(BLT_Mid,BLT_Late) #combine BLT rows
LKT_both <- rbind(LKT_Mid,LKT_Late) #combine LKT rows
#or
Mid_both <- rbind(BLT_Mid,LKT_Mid) #combine mid rows
Late_both <- rbind(BLT_Late,LKT_Late) #combine late rows
```
Then, plot comparisons:
```{r, echo = TRUE, fig.width=12, fig.height=6}
#BLT_plot
BLT_both$Invasion <- factor(BLT_both$Invasion, levels=c("Mid","Late")) #re-order variables
head(BLT_both)
ggplot(data=BLT_both, aes(x=Overlap,y=..scaled.., fill=Invasion, color = Invasion))+ 
  geom_density(alpha=0.5)+
  scale_fill_manual(values=c("gray","black"))+
  scale_color_manual(values=c("gray","black"))+
  scale_x_continuous(limits=c(0,100),breaks=c(0,50,100))+ #for percent instead of proportional overlap
  scale_y_continuous(limits=c(0,1),breaks=c(0,0.5,1))+
  labs(title="Bull trout overlap with lake trout",
       x="\nDiet Overlap (%)",
       y="Scaled Posterior Density\n",
       fill="Invasion")+
  theme_bw()+
  geom_vline(xintercept = quantile(subset(BLT_both,BLT_both$Invasion=="Late")$Overlap, prob = 0.5), linetype="solid", size=2)+ #Late
  geom_vline(xintercept = quantile(subset(BLT_both,BLT_both$Invasion=="Mid")$Overlap, prob = 0.5), linetype="dotted", size=2)+ #Mid
  theme(axis.text = element_text(size = 30, color = "black"),
    axis.title = element_text(size = 30, color = "black"),
    title = element_text(size=30),
    legend.text = element_text(size=30))

#LKT_plot
LKT_both$Invasion <- factor(LKT_both$Invasion, levels=c("Mid","Late")) #re-order variables
head(LKT_both)
ggplot(data=LKT_both, aes(x=Overlap,y=..scaled.., fill=Invasion, color = Invasion))+ 
  geom_density(alpha=0.5)+
  scale_fill_manual(values=c("gray","black"))+
  scale_color_manual(values=c("gray","black"))+
  scale_x_continuous(limits=c(0,100),breaks=c(0,50,100))+ #for percent instead of proportional overlap
  scale_y_continuous(limits=c(0,1),breaks=c(0,0.5,1))+
  labs(title="Lake trout overlap with bull trout",
       x="\nDiet Overlap (%)",
       y="Scaled Posterior Density\n",
       fill="Invasion")+
  theme_bw()+
  geom_vline(xintercept = quantile(subset(LKT_both,LKT_both$Invasion=="Late")$Overlap, prob = 0.5), linetype="solid", size=2)+ #Late
  geom_vline(xintercept = quantile(subset(LKT_both,LKT_both$Invasion=="Mid")$Overlap, prob = 0.5), linetype="dotted", size=2)+ #Mid
  theme(axis.text = element_text(size = 30, color="black"),
    axis.title = element_text(size = 30, color = "black"),
    title = element_text(size=30),
    legend.text = element_text(size=30))
```

#SI Table 4
Finally, produce a summary table for the niche overlap outputs:

```{r}
#Produce SI Table 4
ls_SI_4 <- vector(mode = "list")
#BLT Mid
ls_SI_4$BLT_in_LKT$Mid$CrI_2.5 <- quantile(subset(Mid_both,Mid_both$Species=="Bull trout")$Overlap, prob = 0.025)
ls_SI_4$BLT_in_LKT$Mid$CrI_50 <- quantile(subset(Mid_both,Mid_both$Species=="Bull trout")$Overlap, prob = 0.5)
ls_SI_4$BLT_in_LKT$Mid$Mean <- mean(subset(Mid_both,Mid_both$Species=="Bull trout")$Overlap)
ls_SI_4$BLT_in_LKT$Mid$CrI_97.5 <- quantile(subset(Mid_both,Mid_both$Species=="Bull trout")$Overlap, prob = 0.975)
#BLT Late
ls_SI_4$BLT_in_LKT$Late$CrI_2.5 <- quantile(subset(Late_both,Mid_both$Species=="Bull trout")$Overlap, prob = 0.025)
ls_SI_4$BLT_in_LKT$Late$CrI_50 <- quantile(subset(Late_both,Mid_both$Species=="Bull trout")$Overlap, prob = 0.5)
ls_SI_4$BLT_in_LKT$Late$Mean <- mean(subset(Late_both,Late_both$Species=="Bull trout")$Overlap)
ls_SI_4$BLT_in_LKT$Late$CrI_97.5 <- quantile(subset(Late_both,Late_both$Species=="Bull trout")$Overlap, prob = 0.975)
#LKT Mid
ls_SI_4$LKT_in_BLT$Mid$CrI_2.5 <- quantile(subset(Mid_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.025)
ls_SI_4$LKT_in_BLT$Mid$CrI_50 <- quantile(subset(Mid_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.5)
ls_SI_4$LKT_in_BLT$Mid$Mean <- mean(subset(Mid_both,Mid_both$Species=="Lake trout")$Overlap)
ls_SI_4$LKT_in_BLT$Mid$CrI_97.5 <- quantile(subset(Mid_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.975)
#LKT Late
ls_SI_4$LKT_in_BLT$Late$CrI_2.5 <- quantile(subset(Late_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.025)
ls_SI_4$LKT_in_BLT$Late$CrI_50 <- quantile(subset(Late_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.5)
ls_SI_4$LKT_in_BLT$Late$Mean <- mean(subset(Late_both,Late_both$Species=="Lake trout")$Overlap)
ls_SI_4$LKT_in_BLT$Late$CrI_97.5 <- quantile(subset(Late_both,Mid_both$Species=="Lake trout")$Overlap, prob = 0.975)

stats <- c("Cr_2.5_mid","Cr_50_mid","Mean_mid","Cr_97.5_mid","Cr_2.5_late","Cr_50_late","Mean_late","Cr_97.5_late")
overlaps <- c("BLT_in_LKT","LKT_in_BLT")

SI_4<-data.frame(matrix(unlist(ls_SI_4), #make a dataframe out of the list ls_SI_4
                        nrow=length(overlaps), #one row for each overlap
                        byrow=T),
                 stringsAsFactors=FALSE)
row.names(SI_4) <- overlaps
colnames(SI_4) <- stats
SI_4 <- round(SI_4,1)
library(knitr)
kable(SI_4, align="cccccccc", caption = "Table 6. Summary of isotope niche overlap between lake trout and bull trout in mid- and late-invasion lakes (Article SI Table 4).")
```

#Figure 7, conversion
Conversion uses long-term lake trout and bull trout abundance data to predict the timing of lake trout invasion and persistence of bull trout. Given generic binomial linear form (below), we need R to solve for coefficients ${\beta_0}$ and $\beta$, given lake trout and bull trout abundance data.

$$
Predicted~conversion = \hat{y} = \frac{1}{1+e^{-(\beta_0+\beta\cdot~x)}}
$$
To solve for ${\beta_0}$ and $\beta$, we use binomial linear regression.

```{r, echo=TRUE, warning=FALSE}
rm(list=ls())
graphics.off()
set.seed(1)

data <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/conversion_09302020_2.csv", header=TRUE)
data$conv <- data$LKT/(data$LKT+data$BLT) #calculate empirical conversion for each survey record
data$timestep <- data$year-1955 #conv is ~0 in 1955
tail(data)
conv_lm <- glm(conv ~ timestep, data = data, family = binomial) #lm using timestep
#conv_lm2 <- glm(conv ~ year, data = data, family = binomial) #lm using year
summary(conv_lm) #Estimates of (Intercept) and timestep are the binomial linear regression coefficients. 
```
Specifically,${\beta_0}$ is
```{r}
summary(conv_lm)[["coefficients"]][1]
```
and ${\beta}$ is 
```{r}
summary(conv_lm)[["coefficients"]][2]
```
Producing this:
$$
Predicted~conversion = \hat{y} = \frac{1}{1+e^{-(-3.251+0.091\cdot~x)}}
$$
Then, check the model's fit:
```{r, warning=FALSE,echo=TRUE}
library(rsq)
rsq(conv_lm) #r-squared
```

Finally, plot empirical and predicted conversion versus time:
```{r,echo=TRUE, results=TRUE, warning = FALSE, fig.width=12, fig.height=6}

#Figure 7 (both timeline and year)
ggplot(data = data, aes(x = year,y = conv))+
  geom_jitter(width = 0.2, height = 0, size = 5)+
  stat_smooth(method = "glm",
              method.args = list(family = "binomial"),
              se = TRUE,
              size = 3,
              color = "black",
              fullrange = TRUE)+ #extend prediction beyond empirical data
  theme_bw()+
  labs(y = "Conversion\n")+
  scale_y_continuous(breaks = c(0,0.5,1))+
  scale_x_continuous(name = "\nSurvey year",
                     breaks = c(1955,1970,1985,2000,2015,2030,2045),
                     limits = c(1954,2046),
                     sec.axis = sec_axis(trans = ~.-1955, #conv is ~0 in 1955
                                       name = "Invasion timeline (years)\n",
                                       breaks = c(0,15,30,45,60,75,90)))+
  theme(axis.text.x = element_text(size = 20, color = "black"),
        axis.text.y = element_text(size = 20, color = "black"),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        plot.title = element_text(size = 29))
```

#Figure 8, NMDS ordination for littoral macroinvertebrates

We used non-metric dimensional scaling ordination to measure differences in macroinvertebrate communities among lake trout invasion categories.

```{r, results=FALSE, warning = FALSE, echo=FALSE, message = FALSE}
rm(list=ls())
graphics.off()
set.seed(1)
library(vegan)
library(remotes)
library(ggplot2)
library(dplyr)
library(metR)
library(directlabels)
```

```{r, results=TRUE, warning = FALSE, echo=TRUE, message = FALSE}
#load data
data <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/transect_data.csv", header=TRUE)
head(data[,1:5])
lake <- read.csv("C:/Users/wainr/Documents/nature_08132020/final_datasets/lake_data.csv", header=TRUE)
head(lake)

#community statistics
lake$richness<-rowSums(data[,2:48] > 0) #richness
lake$div<-diversity(data[,2:48]) #diversity
lake$even<-lake$div/log(lake$richness) #evenness
lake$abundance<-rowSums(data[,2:48]) #abundance
species<-data[,2:48]

#calculate nmds scores for communities
nms<-metaMDS(species, k=3)
```
```{r, results=TRUE, warning = FALSE, echo=TRUE, message = FALSE}
nms[["stress"]] #ordination stress (goodness of fit) #http://strata.uga.edu/8370/handouts/nmsTutorial.pdf
#permanova for p-value Pr(>F)
adonis(species~Invasion_category, data = lake)
p_value <- adonis(species~Invasion_category, data = lake)
p_value[["aov.tab"]][["Pr(>F)"]]
#stressplot
stressplot(nms)
```

```{r, results=FALSE, warning = FALSE, echo=TRUE, message = FALSE, fig.show='hide'}
#prep nmds scores for plotting
#ordisurf example: https://chrischizinski.github.io/rstats/ordisurf/
#ordisurf example: https://oliviarata.wordpress.com/2014/07/17/ordinations-in-ggplot2-v2-ordisurf/
ordi <- ordisurf(nms ~ lake$Timestep_years) #create the ordisurf object
ordi.grid <- ordi$grid #extracts the ordisurf object
str(ordi.grid) #it's a list though - cannot be plotted as is
ordi.lake <- expand.grid(x = ordi.grid$x, y = ordi.grid$y) #get x and ys
ordi.lake$z <- as.vector(ordi.grid$z) #unravel the matrix for the z scores
ordi.lake.na <- data.frame(na.omit(ordi.lake)) #gets rid of the nas
ordi.lake.na #looks ready for plotting!
lake.NMDS.data <- lake #Copy lake data, so we can add nms scores. there are other ways of doing this. But this is the way I do it for ease of plotting
lake.NMDS.data$NMDS1 <- nms$points[ ,1] #put the NMDS scores for the plots into a lake dataframe. you could put them into the existing one if you preferred.
lake.NMDS.data$NMDS2 <- nms$points[ ,2]

lake.NMDS.data$Invasion_category <- factor(lake.NMDS.data$Invasion_category,
                                   levels=c('Reference',
                                            'Mid',
                                            'Late')) #re-order variables
```

```{r, echo=TRUE, warning=FALSE, fig.width=12, fig.height=6}
ggplot(lake.NMDS.data, aes(x = NMDS1, y = NMDS2))+
  geom_point(aes(shape = Invasion_category, size=Invasion_category), alpha = 0.6)+ #plots the NMDS points, with shape by invasion
  theme_minimal() + #for aesthetics
  stat_ellipse(aes(x=NMDS1, y=NMDS2, group=Invasion_category, fill = Invasion_category),
               type = "norm",
               geom = "polygon",
               alpha = 0.3,
               level = 0.95)+
  labs(shape = "Invasion", size="Invasion", group="Invasion", fill="Invasion")+#another way to set the labels, in this case, for the color legend
  geom_hline(yintercept=0, linetype="dashed")+
  geom_vline(xintercept=0, linetype="dashed")+
  theme(legend.key = element_blank(),  #removes the box around each legend item
        legend.position = "bottom", #legend at the bottom
        legend.direction = "horizontal",
        legend.box = "horizontal",
        legend.box.just = "centre",
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size=20),
        axis.title.x = element_text(size=20),
        axis.title.y = element_text(size=20),
        plot.title = element_text(size=30),
        legend.text = element_text(size=20),
        legend.title = element_text(size=25),
        panel.background = element_blank(), 
        panel.grid.major = element_blank(),  #remove major-grid labels
        panel.grid.minor = element_blank(),  #remove minor-grid labels
        plot.background = element_blank())
```

# Bibliography
Bunn, S. E., Leigh, C. & Jardine, T. D. Diet-tissue fractionation of $\delta$^15^N  by consumers from streams and rivers. Limnology and Oceanography 58, 765-773 (2013).\
\
de Eyto, E. et al. Characterisation of salmonid food webs in the rivers and lakes of an Irish peatland ecosystem. Proceedings of the Royal Irish Academy 120, 1-17 (2020).\
\
Guzzo, M. M., Haffner, G. D., Sorge, S., Rush, S. A. & Fisk, A. T. Spatial and temporal variabilities of $\delta$^13^C and $\delta$^15^N within lower trophic levels of a large lake: implications for estimating trophic relationships of consumers. Hydrobiologia 675, 41-53 (2011).\
\
Hobson, K. A., Ofukany, A., Soto, D. X. & Wassenaar, L. I. An isotopic baseline ($\delta$^13^C, $\delta$^15^N) for fishes of Lake Winnipeg: Implications for investigating impacts of eutrophication and invasive species. Journal of Great Lakes Research 38, 58-65 (2012).\
\
Post, D. M. Using stable isotopes to estimate trophic position: models, methods, and assumptions. Ecology 83, 703-718 (2002).\
\
Rogosch, J. S. & Olden, J. D. Invaders induce coordinated isotopic niche shifts in native fish species. Can. J. Fish. Aquat. Sci. 1-11 (2020).\
\
Syväranta, J. & Jones, R. I. Changes in feeding niche widths of perch and roach following biomanipulation, revealed by stable isotope analysis. Freshwater Biology 53, 425-434 (2008).\
\
Vander Zanden, M. J. & Rasmussen, J. B. Primary Consumer $\delta$^13^C and $\delta$^15^N and the Trophic Position of Aquatic Consumers. Ecology 80, 1395-1404 (1999).








