Supp_file_bufferingExpt_v1.Rmd

---
title: "Data for: Do experimental pH increases alter the structure and function of a lowland tropical stream?"
author: "Nicholas S. Marzolf"
date: "Updated January 26, 2022"
output: 
  pdf_document:
---

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

# Introduction
RMarkdown file to accompany "Do experimental pH increases alter the structure and function of a lowland tropical stream?". This file contains the data files to recreate Tables 2-4 and Figures 2-6. The section for each figure loads in the appropriate data file to conduct the analysis and/or recreate the figure.

To recreate this HTML, load the data files into a folder named 'Data' at the same path as the .rmd file, open the .rmd file, and click 'Knit'.

# Load packages
```{r load-packages, echo = TRUE, message=FALSE, warning = FALSE}
library(readxl)    # read excel files
library(dplyr)     # data manipulation
library(tidyr)     # data manipulation
library(ggplot2)   # visualization
library(lubridate) # manage timestamps
library(cowplot)   # multipane figures
library(agricolae) # Tukey HSD test
library(ggsci)     # ggplot colors
library(moments)   # normality
library(emmeans)   # Estimated marginal means
library(car)       # Anova() function
library(effsize)   # Hedge's g effect size for low sample sizes
library(ARTool)    # Align rank transform test
library(vegan)     # Non-metric dimensional scaling
library(ggrepel)   # visualization
library(scales)    # plotting date-times
library(pander)    # session info

source('C:/Users/nmarz/Desktop/Research/R code/theme_nick.R')
theme_set(theme_nick())

```

# Table 2
Nutrients collected weekly in the 4 reaches during the 5 week experiment. 
```{r, nut-summary, warning=FALSE, echo=TRUE}
nutrients <- read.csv('Data/nutrients.csv')     # nutrient data

nutrients %>%
  mutate(din = NOx_N + NH4_N,
         din_srp = (din/14.007)/(PO4_P/30.9737)
  ) %>%
  group_by(Reaches) %>%
  summarise(mean_NO3 = mean(NOx_N*1000),
            se_NO3 = sd(NOx_N*1000)/sqrt(length(NOx_N)),
            corNO3 = cor(y = NOx_N*1000, x = rep(1:6, 3), method = "kendall"),
            mean_NH4 = mean(NH4_N*1000),
            se_NH4 = sd(NH4_N*1000)/sqrt(length(NH4_N)),
            corNH4 = cor(y = NH4_N*1000, x = rep(1:6, 3), method = "kendall"),
            mean_PO4 = mean(PO4_P*1000),
            se_PO4 = sd(PO4_P*1000)/sqrt(length(PO4_P)),
            corPO4 = cor(y = PO4_P*1000, x = rep(1:6, 3), method = "kendall"),
            mean_NP = mean(din_srp),
            se_NP = sd(din_srp)/sqrt(length(din_srp)),
            corNP = cor(y = din_srp, x = rep(1:6, 3), method = "kendall"))

```

# Table 3
Chironomid experiment data
```{r chiro-survial, message=FALSE, warning=FALSE}
chiroGrowth <- read.csv('Data/chiroGrowth.csv') # Chironomid growth

chiroGrowth$Reach = factor(chiroGrowth$Reach)

chiroGrowth %>% 
  group_by(Reach) %>%
  summarise(biomass_mean = mean(biomassFin_mg, na.rm = TRUE),
            biomass_se = sd(biomassFin_mg, na.rm = TRUE)/sqrt(length(biomassFin_mg)),
            per_biomass_mean = mean(delBiomass_per, na.rm = TRUE),
            per_biomass_se = sd(delBiomass_per, na.rm = TRUE)/sqrt(length(biomassFin_mg)),
            per_surv_mean = mean((numSurv/6)*100, na.rm = TRUE),
            per_surv_se = sd((numSurv/6)*100, na.rm = TRUE)/sqrt(length(biomassFin_mg)))
```


# Figure 2
Rainfall and pH time-series in the four reaches
```{r Fig2, echo=TRUE, warning=FALSE, results='hide'}

# load the sensor data
CarDown <- read.csv('Data/CarDown.csv') # YSI EXO1 data- Carapa downstream
CarUp <- read.csv('Data/CarUp.csv')     # YSI EXO1 data- Carapa upstream
CarUp$Timestamp <- as.POSIXct(CarUp$Timestamp)
ArbDown <- read.csv('Data/ArbDown.csv') # YSI EXO1 data- ArbSeep downstream
ArbUp <- read.csv('Data/ArbUp.csv')     # YSI EXO1 data- ArbSeep upstream

# hyetographs for the 4 reaches
# the upstream rainfall is the same for both sites
CarUp_rain <- ggplot(CarUp %>% 
                       group_by(Timestamp = as.POSIXct(date(Timestamp))) %>%
                       summarise(Rainfall = sum(Rainfall, na.rm = TRUE)), 
                     aes(x = Timestamp, 
                         y = Rainfall))+
  geom_bar(stat = 'identity', color = 'black')+
  ylim(100, 0)+
  scale_x_datetime(date_breaks = "10 day",
                   date_labels = "%m/%d/%y")+
  labs(y = 'Rainfall (mm)')+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size = 10, face = 'plain'))
CarUp_rain

ArbDown_rain <- ggplot(ArbDown %>% 
                         group_by(Timestamp = date(Timestamp)) %>%
                         summarise(Rainfall = sum(Rainfall, na.rm = TRUE)), 
                       aes(x = as.POSIXct(Timestamp), 
                           y = Rainfall))+
  geom_bar(stat = 'identity', 
           color = 'black')+
  ylim(100, 0)+
  ylab("Rainfall (mm)")+
  scale_x_datetime(date_breaks = "5 day",
                   date_labels = "%m/%d/%y",
                   limits = as.POSIXct(c("2018-06-18", "2018-07-18")))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_text(size = 10))
ArbDown_rain


# pH time-series in the 4 reaches
pH_downstream <- ggplot()+
  geom_point(data = ArbDown, alpha = 0.7,
             aes(x = as.POSIXct(Timestamp), 
                 y = as.numeric(pH), 
                 color = 'ArbSeep'))+
  geom_point(data = CarDown, alpha = 0.7,
             aes(x = as.POSIXct(Timestamp), 
                 y  = as.numeric(pH), 
                 color = 'Carapa'))+
  ylim(4, 8)+
  ylab('pH')+
  xlab('Date')+
  scale_x_datetime(date_breaks = "7 day",
                   date_labels = "%m/%d/%y",
                   limits = as.POSIXct(c("2018-06-18", "2018-07-18")))+
  scale_color_manual(name = 'Downstream Reach',
                     values = c('#868686FF', '#CD534CFF'))+
  theme(legend.position = c(0.22, 0.2),
        legend.background = element_rect(fill = "white", color = "black"),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        axis.text.x = element_text(size = 10))

pH_upstream <- ggplot()+
  geom_point(data = ArbUp, alpha = 0.7,
             aes(x = as.POSIXct(Timestamp), 
                 y = pH, 
                 color = 'ArbSeep'))+
  geom_point(data = CarUp, alpha = 0.7,
             aes(x = as.POSIXct(Timestamp), 
                 y  = pH, 
                 color = 'Carapa'))+
  ylim(4, 8)+
  xlab('Date')+
  scale_color_manual(name = 'Upstream Reach',
                     values = c('#0073C2FF', '#EFC000FF'))+
  scale_x_datetime(date_breaks = "10 day",
                   date_labels = "%m/%d/%y")+
  theme(legend.position = c(0.8, 0.8),
        legend.background = element_rect(fill = "white", color = "black"),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        axis.text.x = element_text(size = 10))

# paired histograms
Car_pairedHist <- ggplot()+
  geom_histogram(data = CarUp, bins = 200, alpha = 0.75,
                 aes(x = pH, fill = 'Upstream'))+
  geom_histogram(data = CarDown, bins = 200, alpha = 0.75,
                 aes(x = as.numeric(pH), fill = 'Downstream'))+
  xlim(4, 8)+
  ylim(0, 600)+
  ylab("Count")+
  scale_fill_manual(labels = c("Upstream", "Downstream"),
                    breaks = c("Upstream", "Downstream"),
                    values = c('#EFC000FF', '#CD534CFF'),
                    name = 'Carapa')+
  theme(legend.justification = c(0.75,0.75),
        legend.position = c(0.9,0.9),
        legend.text = element_text(face = 'plain', size = 10),
        legend.title = element_text(size = 10),
        legend.background = element_rect(color = 'black'),
        axis.text = element_text(size = 8))
Car_pairedHist

Arb_pairedHist <- ggplot()+
  geom_histogram(data = ArbUp, bins = 200, alpha = 0.75,
                 aes(x = pH, fill = 'Upstream'))+
  geom_histogram(data = ArbDown, bins = 200, alpha = 0.75,
                 aes(x = as.numeric(pH), fill = 'Downstream'))+
  xlim(4, 8)+
  ylim(0, 600)+
  ylab("Count")+
  xlab("pH")+
  scale_fill_manual(labels = c("Upstream", "Downstream"),
                    breaks = c("Upstream", "Downstream"),
                    values = c('#0073C2FF', '#868686FF'),
                    name = 'ArbSeep')+
  theme(legend.justification = c(0.75,0.75),
        legend.position = c(0.9,0.9),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10),
        legend.background = element_rect(color = 'black'),
        axis.text = element_text(size = 10))
Arb_pairedHist

# combine into Fig 2
fig2 <- plot_grid(CarUp_rain, ArbDown_rain, 
                  pH_upstream, pH_downstream, 
                  Arb_pairedHist, Car_pairedHist,
                  nrow = 3, ncol = 2,
                  rel_heights = c(0.5,1,1),
                  labels = 'auto', label_fontface = 'plain',
                  label_size = 16,
                  hjust = c(-6.5, -6.5,-7,-6.5,-6.5,-13.5),
                  vjust = c(7,7,1.5,1.5,1.5,1.6),
                  align = 'v')
fig2
```

# Figure 3
Leaf litter and coarse woody debris decomposition rates.
```{r Fig3, echo = TRUE, warning = FALSE}
OM <- read.csv('Data/OM.csv')                   # leaf litter and coarse woody debris decomposition rate data

OM$Reach = factor(OM$Reach,
                  levels = c('ArbSeep Upstream',
                             'ArbSeep Downstream',
                             'Carapa Upstream',
                             'Carapa Downstream'))

# leaf litter boxplot; letters are Tukey post-hoc test
plot_LL_k <- OM %>%
  select(Stream, Treatment, Reach, Bag_mesh, k_LL_d) %>%
  group_by(Reach, Bag_mesh) %>%
  summarise(mean_k = mean(k_LL_d, na.rm = TRUE),
            se_k = sd(k_LL_d, na.rm = TRUE)/length(k_LL_d)) %>%
  ggplot(., aes(x = Reach, y = mean_k, 
                color = Bag_mesh))+
  geom_point(size = 3,
             position=position_dodge(width=0.5))+
  geom_errorbar(aes(ymin = mean_k - se_k,
                    ymax = mean_k + se_k),
                width = 0.3, 
                position=position_dodge(width=0.5))+
  ylim(.035, 0.05)+
  scale_color_jco(name = 'Mesh Size')+
  ylab(expression(paste(k[LL], ' (',d^-1,')')))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = 'none')

# woody debris boxplot
plot_CWD_k <- OM %>%
  select(Stream, Treatment, Reach, Bag_mesh, k_CWD_yr) %>%
  group_by(Reach, Bag_mesh) %>%
  summarise(mean_k = mean(k_CWD_yr, na.rm = TRUE),
            se_k = sd(k_CWD_yr, na.rm = TRUE)/length(k_CWD_yr)) %>%
  ggplot(., aes(x = Reach, y = mean_k, 
                color = Bag_mesh))+
  geom_point(size = 3,
             position=position_dodge(width=0.5))+
  geom_errorbar(aes(ymin = mean_k - se_k,
                    ymax = mean_k + se_k),
                width = 0.3, 
                position=position_dodge(width=0.5))+
  #ylim(.035, 0.05)+
  scale_color_jco(name = 'Mesh Size')+
  ylab(expression(paste(k[CWD], ' (',y^-1,')')))+
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = c(0.8, 0.75),
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "black"))

# combine plots
plot_grid(plot_LL_k, plot_CWD_k, 
          ncol = 2, 
          labels = 'auto',
          align = 'hv')
```

# Figure 4
Macroinvertebrate abundance per g of organic matter AFDM
```{r macros-long-plot, warning=FALSE, echo=TRUE}
macros <- read.csv('Data/macroinverts.csv')     # macroinvertebrate abundance for W1 and W5 data

# assign NAs to 0
macros[is.na(macros)] <- 0

# Family counts
glimpse(macros)

# Make factors
macros$Week = factor(macros$Week)
macros$Organic_Type = factor(macros$Organic_Type)
macros$Mesh_Size = factor(macros$Mesh_Size)

# Subset data
macro_fam <- macros %>%
  mutate(Hydracarina = Hydracarina + Water.Mite) %>%
  select(-Water.Mite) %>%
  select(Sample:Reach, Chironomidae:Hydracarina) %>%
  mutate(Reaches = factor(interaction(Site, Reach, sep = ' '))) %>%
  rename(OM = Organic_Type,
         Mesh = Mesh_Size)

macro_fam$Reaches <- ordered(macro_fam$Reaches, 
                             levels = c("Carapa Upstream", "Carapa Downstream", 
                                        "ArbSeep Upstream", "ArbSeep Downstream"))


# convert data to mean abundance by family for LL litter bags
mean_fam_abun_LL = macro_fam %>% 
  dplyr::rowwise(Sample) %>% 
  dplyr::mutate(sum = rowSums(across(Chironomidae:Hydracarina))) %>%
  dplyr::group_by(Week, OM, Mesh, Reaches) %>%
  dplyr::summarise(mean = mean(sum),
                   se = sd(sum)/sqrt(length(sum)),
                   len = length(sum)) %>%
  dplyr::filter(OM == 'Leaf')

# these data are the leaf litter AFDM remaining at Week 5 from the OM data file. Here, I load them in as a new column
mean_fam_abun_LL$LL_AFDM = c(4.4, 3.97, 4.28, 4.33, 4.07, 4.07, 4.4, 3.93,
                             0.59, 1.14, 0.91, 1.22, 0.63, 1.13, 0.91, 1.12)

# calculate the mean family abundance per g AFMD remaining
mean_fam_abun_LL = mean_fam_abun_LL %>%
  mutate(mean_abund_AFDM = mean/LL_AFDM,
         se_abund_AFDM = se/LL_AFDM)

# plot mean family abundance per g AFDM remaining
abundance_LL = ggplot(data = mean_fam_abun_LL, 
                      aes(x = Week, y = mean_abund_AFDM))+
  geom_point(size = 3)+
  geom_errorbar(aes(ymin = mean_abund_AFDM - se_abund_AFDM, 
                    ymax = mean_abund_AFDM + se_abund_AFDM),
                width = 0.2)+
  facet_grid(Mesh ~ Reaches)+
  scale_y_continuous(limit = c(0, 200),
                     breaks = c(0, 50, 100, 150, 200),
                     labels = c(0, 50, 100, 150, 200))+
  #scale_y_log10(limit = c(0.1, 300),
  #              breaks = c(0.1, 1, 10, 100, 300),
  #              labels = c(0.1, 1, 10, 100, 300))+
  ylab(expression(paste("Mean abundance (# g LL ", AFDM^-1,")")))+
  theme(strip.background = element_rect(color = "black", fill = "white"),
        strip.text = element_text(color = "black", size = 12))

# repeat the calculation for CWD litter bags
mean_fam_abun_CWD <- macro_fam %>% 
  dplyr::rowwise(Sample) %>% 
  mutate(sum = rowSums(across(Chironomidae:Hydracarina))) %>%
  group_by(Week, OM, Mesh, Reaches) %>%
  summarise(mean = mean(sum),
            se = sd(sum)/sqrt(length(sum))) %>%
  filter(OM == 'Wood')

# these are the CWD AFDM data, as described above
mean_fam_abun_CWD$CWD_AFDM <- c(36.8, 30.8, 38.3, 34.7, 40.2, 30.7, 40.9, 34.8,
                                35.6, 28.6, 25.8, 28.4, 23.7, 29.2, 24.4, 29.6)

# calculate mean family abundance per g AFDM
mean_fam_abun_CWD <- mean_fam_abun_CWD %>%
  mutate(mean_abund_AFDM = mean/CWD_AFDM,
         se_abund_AFDM = se/CWD_AFDM)

# create the plot
abundance_CWD <- ggplot(data = mean_fam_abun_CWD, 
                        aes(x = Week, y = mean_abund_AFDM))+
  geom_point(size = 3)+
  geom_errorbar(aes(ymin = mean_abund_AFDM - se_abund_AFDM, 
                    ymax = mean_abund_AFDM + se_abund_AFDM),
                width = 0.2)+
  facet_grid(Mesh ~ Reaches)+
  scale_y_continuous(limit = c(-1, 8),
                     breaks = c(0, 2, 4, 6, 8),
                     labels = c(0, 2, 4, 6, 8))+
  ylab(expression(paste("Mean abundance (# g CWD  ", AFDM^-1,")")))+
  theme(strip.background = element_rect(color = "black", fill = "white"),
        strip.text = element_text(color = "black", size = 12))

# combine OM abundance plots
Fig4 <- plot_grid(abundance_LL, abundance_CWD, 
                  align = 'hv',
                  nrow = 2, 
                  labels = 'auto', label_fontface = 'plain', label_size = 18,
                  hjust = -2, vjust = 0.75)
Fig4

```

# Figure 5
Species richness
```{r spec-richness, warning=FALSE, echo=TRUE}
macro_fam$numFam = specnumber(select(macro_fam, Chironomidae:Hydracarina))

Fig5 <- ggplot(macro_fam %>% 
                 group_by(Week, OM, Reaches, Mesh) %>%
                 summarise(mean_rich = mean(numFam, na.rm = TRUE),
                           se_rich = sd(numFam, na.rm = TRUE)/length(numFam)), 
               aes(x = Week, 
                   y = mean_rich, 
                   color = OM))+
  geom_point(size = 3,
             position=position_dodge(width=0.5))+
  geom_errorbar(aes(ymin = mean_rich - se_rich,
                    ymax = mean_rich + se_rich),
                width = 0.2,
                position=position_dodge(width=0.5))+
  facet_grid(Mesh~Reaches)+
  ylab("Macroinvertebrate Richness")+
  scale_y_continuous(limits = c(0, 12),
                     breaks = c(0, 6, 12),
                     labels = c(0, 6, 12))+
  scale_color_jco(name = 'Litter Type')+
  theme(strip.background = element_rect(color = "black", fill = "white"),
        strip.text = element_text(color = "black"))
Fig5
```

# Figure 6
Non-metric dimensional scaling of the various treatments in the experiment
```{r NMDS, warning=FALSE, echo=TRUE}
macro_fam_W5 <- macro_fam %>%
  filter(Week == 'W5') %>%
  mutate_at(vars(Chironomidae:Hydracarina), 
            funs(log10(1+.)))

nmds <- metaMDS(macro_fam_W5 %>%
                  select(Chironomidae:Hydracarina),
                distance = 'bray',     # Bray-Curtis dissimilarity index
                k = 3,                 # specify 3 dimensions
                autotransform = FALSE) # we did this manually
nmds$stress # 0.17, not too bad
```

Now let's begin the visualization process. First, we will extract points from the ordination and bind them to the experimental treatment data not run in the NMDS matrix. We will then look at the taxa that are driving the dissimilarity in the matrix and plot the vectors on the NMDS points
```{r nmds-env-vectors, warning=FALSE, echo = TRUE}
df <- data.frame(x = nmds$points[,1],
                 y = nmds$points[,2],
                 z = nmds$points[,3])
df<-cbind(macro_fam_W5[,c(1:6, 40)],
          df)

fit<-(envfit(nmds, macro_fam_W5[,-c(1:6)],
             perm = 9999))
scrs <- as.data.frame(scores(fit, 'vectors'))
scrs$pvals = fit$vector$pvals
scrs_sig<-subset(scrs, pvals <= 0.05)
scrs_sig$env.variables = row.names(scrs_sig)

# plot nmds with significant family vectors
nmds_vect <- ggplot(data = df, aes(x = x, y = y))+                             
  geom_point(size = 3) + 
  geom_segment(data = scrs_sig,
               aes(x = 0, xend = NMDS1, 
                   y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm")), 
               colour = "black") +
  geom_label_repel(data = scrs_sig, 
                   aes(NMDS1, NMDS2, label = env.variables),
                   size = 3)+
  labs(x = "NMDS axis 1", y = "NMDS axis 2")+
  annotate("text", label = "Stress = 0.17", x = -0.85, y = -1.4, fontface = 3)+
  lims(x = c(-1, 1), y = c(-1.5, 1))
nmds_vect
```

Now we will calculate 95% confidence interval ellipses for the four study reaches, two mesh bag sizes, and 2 organic matter types, and plot the ellipses on the NMDS to determine the similarity of the assemblage as explained by those treatments. We will need the veganCovEllipse.R function to help us. Briefly, we will create a dummy data frame that is run through a for loop that calculates the 95% CI ellipse for the treatment factor of interest. For the sake of space, I will show the process for the first 2 axes (x and y) but can repeated to visualize the NMDS using other axes (substitute y for x to change from axis 2 to 1, and substitute z for y to change from axis 2 to axis 3):
```{r, nmds-ellipses, warning=FALSE, echo=TRUE}

# source('veganCovEllipse.R')

ell_reaches_12 <- data.frame() # create a dummy data frame for ellipses based on reach for axes 1 and 2

for(g in unique(df$Reaches)) {
  ell_reaches_12 <- rbind(ell_reaches_12,
                          cbind(
                            as.data.frame(
                              with(df[df$Reaches == g,],
                                   veganCovEllipse(cov.wt(cbind(x, y),
                                                          wt = rep(1/length(x),
                                                                   length(x)))$cov,
                                                   center = c(mean(x),
                                                              mean(y)
                                                   )
                                   )
                              )
                            ), 
                            Reaches = g)
  )
}

# repeat for OM and Mesh treatments
ell_OM_12 <- data.frame()

for(g in unique(df$OM)) {
  ell_OM_12 <- rbind(ell_OM_12,
                     cbind(
                       as.data.frame(
                         with(df[df$OM == g,],
                              veganCovEllipse(cov.wt(cbind(x, y),
                                                     wt = rep(1/length(x),
                                                              length(x)))$cov,
                                              center = c(mean(x),
                                                         mean(y)
                                              )
                              )
                         )
                       ),
                       OM = g)
  )
}

ell_mesh_12 <- data.frame()
for(g in unique(df$Mesh)) {
  ell_mesh_12 <- rbind(ell_mesh_12,
                       cbind(
                         as.data.frame(
                           with(df[df$Mesh == g,],
                                veganCovEllipse(cov.wt(cbind(x, y),
                                                       wt = rep(1/length(x),
                                                                length(x)))$cov,
                                                center = c(mean(x),
                                                           mean(y)
                                                )
                                )
                           )
                         ),
                         Mesh = g)
  )
}

```

Now we will plot the NMDS using the 95% CI ellipses we just calculated for reaches, organic matter, and mesh size. Then, we will combine the plots into Figure 6
```{r nmds-plots, warning = FALSE, echo = TRUE}
# by reach
nmds_reaches <- ggplot(data = df, 
                       aes(x = x, y = y))+                             
  geom_point(size = 2, aes(color = Reaches)) + 
  geom_path(data = ell_reaches_12, 
            aes(x = x, y = y, color = Reaches),
            size = 1)+
  labs(x = "NMDS axis 1", 
       y = "NMDS axis 2")+
  guides(color = guide_legend(ncol = 2))+
  theme(legend.position = c(0.6, 0.1),
        legend.background = element_blank(),
        legend.direction = 'horizontal')+
  scale_color_jco(name = element_blank())+
  lims(x = c(-1, 1), y = c(-1.5, 1))
nmds_reaches

# by organic matter type
nmds_OM <- ggplot(data = df, 
                  aes(x = x, y = y))+                             
  geom_point(size = 2, aes(color = OM)) + 
  geom_path(data = ell_OM_12, 
            aes(x = x, y = y, color = OM),
            size = 1)+
  labs(x = "NMDS axis 1", 
       y = "NMDS axis 2")+
  guides(color = guide_legend(ncol = 2))+
  theme(legend.position = c(0.6, 0.1),
        legend.background = element_blank(),
        legend.direction = 'horizontal')+
  scale_color_jco(name = element_blank())+
  lims(x = c(-1, 1), y = c(-1.5, 1))
nmds_OM

# by mesh size
nmds_mesh <- ggplot(data = df, 
                    aes(x = x, y = y))+                             
  geom_point(size = 2, aes(color = Mesh)) + 
  geom_path(data = ell_mesh_12, 
            aes(x = x, y = y, color = Mesh),
            size = 1)+
  labs(x = "NMDS axis 1", 
       y = "NMDS axis 2")+
  guides(color = guide_legend(ncol = 2))+
  theme(legend.position = c(0.6, 0.1),
        legend.background = element_blank(),
        legend.direction = 'horizontal')+
  scale_color_jco(name = element_blank())+
  lims(x = c(-1, 1), y = c(-1.5, 1))
nmds_mesh

# combine plots
Fig6 <- plot_grid(nmds_vect, nmds_reaches,
                  nmds_OM, nmds_mesh,
                  align = 'hv',
                  labels = 'auto', label_fontface = 'plain',
                  hjust = -7.5)
Fig6
```

# Table 4
Effect size values, measured as Hedges' g.
```{r effect-sizes,warning=FALSE, echo=TRUE}

# Macroinvertebrate abundance
fam_long <- 
  macro_fam %>%
  pivot_longer(cols = Chironomidae:Hydracarina,
               names_to = 'family', values_to = 'count') %>%
  dplyr::select(Sample:Reach, Reaches, family, count)

cohen.d(data = fam_long %>% 
          dplyr::filter(Site == "Carapa"),
        count ~ Reach, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

cohen.d(data = fam_long %>% 
          dplyr::filter(Site == "ArbSeep"),
        count ~ Reach, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

# Macroinvertebrate richness
cohen.d(data = macro_fam %>% 
          dplyr::filter(Site == "Carapa"),
        numFam ~ Reach, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

cohen.d(data = macro_fam %>% 
          dplyr::filter(Site == "ArbSeep"),
        numFam ~ Reach, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

# leaf litter decomposition rates
cohen.d(data = dplyr::filter(OM, Stream == "Carapa"),
        k_LL_d ~ Treatment, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

cohen.d(data = dplyr::filter(OM, Stream == "Arb Seep"),
        k_LL_d ~ Treatment, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

# CWD decomposition rates
cohen.d(data = dplyr::filter(OM, Stream == "Carapa"),
        k_CWD_yr ~ Treatment, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

cohen.d(data = dplyr::filter(OM, Stream == "Arb Seep"),
        k_CWD_yr ~ Treatment, na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

# Chironomid %growth
cohen.d(data = chiroGrowth %>% 
          dplyr::filter(Reach %in% c("Carapa Upstream", "Carapa Downstream")),
        delBiomass_per ~ Reach, 
        na.rm = TRUE,
        hedges.correction = TRUE,
        conf.level = 0.95)

cohen.d(data = chiroGrowth %>% 
          dplyr::filter(Reach %in% c("ArbSeep Upstream", "ArbSeep Downstream")),
        delBiomass_per ~ Reach,
        hedges.correction = TRUE,
        conf.level = 0.95)
``` 

# Supplementary Tables and Figures
## Table S1
```{r supptable1, }
LL_nest <- lm(k_LL_d ~ (Stream/Treatment)*Bag_mesh,
              data = OM)
summary(LL_nest)
anova(LL_nest)
Anova(LL_nest, type = 'III')

print(HSD.test(LL_nest, trt = 'Bag_mesh'))
print(HSD.test(LL_nest, trt = 'Stream'))
print(HSD.test(LL_nest, trt = 'Treatment'))

CWD_nest <- lm(k_CWD_yr ~ (Stream/Treatment)*Bag_mesh,
               data = OM)
summary(CWD_nest)
Anova(CWD_nest, type = 'III')

```

## Table S2
```{r supptable2, warming = FALSE, echo = TRUE}
biomass_nested <- lm(data = chiroGrowth,
                     log10(delBiomass_per) ~ Stream/Trt_Reach)
summary(biomass_nested)
Anova(biomass_nested, type = 'III')

biomass_chg_nested <- lm(data = chiroGrowth,
                         log10(biomassFin_mg) ~ (Stream/Trt_Reach))
summary(biomass_chg_nested)
Anova(biomass_chg_nested, type = 'III')
```

## Fig S1
```{r suppfig1, warning=FALSE, echo=TRUE}
family_long <- 
  macro_fam %>%
  #mutate(Hydracarina = Hydracarina + Water.Mite) %>%
  #select(-Water.Mite) %>% 
  pivot_longer(cols = Chironomidae:Hydracarina,
               names_to = 'family', values_to = 'count') %>%
  dplyr::select(Sample:Mesh, Reaches, family, count)


Fig_S1 <- ggplot(data = family_long %>% 
                   dplyr::group_by(family, Week) %>% 
                   dplyr::summarise(sum = sum(count)),
                 aes(x = reorder(family, sum), y = sum))+
  geom_bar(stat = 'identity', color = 'black')+
  geom_text(aes(label = sum), hjust = 2, color = 'white')+
  coord_flip()+
  facet_grid(.~Week)+
  scale_y_log10(limits = c(1, 5000),
                breaks = c(1, 10, 100, 1000, 5000),
                labels = c(1, 10, 100, 1000, 5000))+
  xlab('Taxon')+
  ylab("Abundance")
Fig_S1
```

## Fig S2
```{r suppfig2, warning=FALSE, echo=TRUE}
k_pH = read_excel("Data/decomp and pH.xlsx", 
                  sheet = "Sheet2")
ggplot(k_pH, aes(x = pH, y = k))+
  geom_point(aes(color = Species),
             size = 3)+
  scale_color_jco()+
  ylim(0, 0.15)+
  xlim(4, 7)+
  ylab(expression(paste("k (", d^-1,")")))+
  theme(legend.text = element_text(face = 'italic'))
```

# Session Info
```{r session-info}
pander(sessionInfo())
```