Archeratal_Ircinia_seagrass_nutrient_supplement.Rmd

---
title: "Supplement 2"
output:
  bookdown::pdf_book
---

```{r setup, include=FALSE,warning=FALSE,message=FALSE}
knitr::opts_chunk$set(echo = TRUE)
if(!require(tidyverse))install.packages("tidyverse");library(tidyverse)
if(!require(car))install.packages("car");library(car)
if(!require(lmerTest))install.packages("lmerTest");library(lmerTest)
if(!require(kableExtra))install.packages("kableExtra");library(kableExtra)

theme_set(theme_bw())
source("scripts/03_reimport.R")#imports all the data sets
#load data
mn<-readxl::read_xlsx("Original_data/missing_nutrients.xlsx")
sg_nuts<-bind_rows(sg_nuts,mn)%>%
  select(-mg)%>%
  group_by(treatment,plot,dist,sampling)%>%
  summarize(PN=sum(PN,na.rm = TRUE),PC=sum(PC,na.rm = TRUE),PP=sum(PP,na.rm = TRUE))

#make seagrass nutriet data frame with nvalues filled in
sgn<-sg_nuts %>%
  filter(PC >0 & PN >0 & PP >0)%>%
  mutate(cn=PC/PN,
         cp=PC/PP,
         np=PN/PP)%>%
  pivot_longer(PN:np,names_to = "nut",values_to = "nvalue")%>%
  filter(!is.na(nvalue))

# there are some outliers in PC and PN- look at what these are
filter(sgn,nut=="PC" & nvalue>55)
filter(sgn,nut=="PN" & nvalue>2.5)
# carbon and nitrogen are calculated from the same sample and plot 6 (sampling 2 dist 0), 
# plot 7 sampling 2 dist 0, and plot 8 sampling 2 dist 1
# are outliers for both- makes me think something besides seagrass got in there.
# I'm going to remove these points before going further. 
sgn_no<-bind_rows(sgn %>%
  filter(plot==7 & sampling==2 & dist==0),
  sgn %>%
    filter(plot==8 & sampling==2 & dist==1))

sgn<-sgn %>%
  anti_join(sgn_no)



# create a dataset with offset values and season
sv<-filter(sgn, sampling == 1) %>% 
  rename(snv = nvalue)%>%select(treatment,plot,dist,snv,nut)

dsgn<-left_join(sgn,sv)%>%
  mutate(delta=nvalue-snv,
         season=case_when(
           sampling==2~"Summer",
           sampling==3~"Winter",
           sampling==4~"Summer",
           sampling==5~"Winter"),
         yr=case_when(
           sampling==2~1,
           sampling==3~1,
           sampling==4~2,
           sampling==5~2),
         mnths=case_when(
           sampling==1~0,
           sampling==2~1,
           sampling==3~5,
           sampling==4~12,
           sampling==5~17),
         treatment=case_when(
           treatment=="fake"~"Structure",
           treatment=="blank"~"Control",
           treatment=="real"~"Sponge"))%>%filter(sampling!=1)
```
# Nutrient Analysis  

Unfortunately, some samples for nutrient analysis were lost after collection, resulting in an unbalanced sampling design and even smaller sample sizes(Table \@ref(tab:ntable)). 
```{r ntable,echo=FALSE,warning=FALSE,message=FALSE}
# look at sample sizes to see what makes the most sense
ntable<-dsgn %>%
    filter(nut %in% c("PP", "PC", "PN")) %>%
    group_by(treatment, yr, season, dist, nut) %>%
    summarize(ns = n()) %>%
    pivot_wider(
      names_from = c(treatment, nut),
      values_from = ns) %>%
    ungroup() %>%
    select(
      `Distance (m)` = dist,
      Control_PC,
      Control_PN,
      Control_PP,
      Structure_PC,
      Structure_PN,
      Structure_PP,
      Sponge_PC,
      Sponge_PN,
      Sponge_PP) #%>%
   # # distinct() %>%
   #  htmlTable(
   #    rgroup = c("Summer", "Winter", "Summer", "Winter"),
   #    n.rgroup = c(4, 4, 4, 4),
   #    tspanner = c("Year 1", "Year 2"),
   #    n.tspanner = c(8, 8),
   #    cgroup = c("", "Control", "Structure", "Sponge"),
   #    n.cgroup = c(1, 3, 3, 3),
   #    header=c("Distance(m)","%C","%N","%P","%C","%N","%P","%C"# ,"% N","% P"))

kbl(ntable, caption = "Sample sizes per treatment, nutrient, distance and sampling. The samples used for the analysis presented in the manuscript text are highlighted. Samples were retained for all plots at all distances from the pre-experiment sample period, which was conducted in summer and is not shown here.", booktabs = T,
    col.names = c("Distance(m)","%C","%N","%P","%C","%N","%P","%C","% N","% P")) %>%
  kable_styling(latex_options = "scale_down")%>%
add_header_above(c(" "=1, "Control" = 3, "Structure" = 3,
                   "Sponge" = 3)) %>%
pack_rows("Year 1", 1, 8) %>%
pack_rows("Year 2", 9, 16) %>%
pack_rows("Summer", 1, 4) %>%
pack_rows("Winter", 5, 8)%>%
pack_rows("Summer", 9, 12) %>%
pack_rows("Winter", 13, 16) %>%
  row_spec(9,background="lightgray")
```

As a result our ability to detect differences in one global model was greatly reduced. 
Therefore, we first examined each nutrient for a seasonal effect using separate mixed effects models with season as the fixed effect and plot as the random effect.
There was a significant seasonal effect for all nutrients (Nitrogen, Table \@ref(tab:nseason); Phosphorus, Table \@ref(tab:pseason); Carbon, Table \@ref(tab:cseason)). 

```{r nseason,echo=FALSE,message=FALSE,warning=FALSE}
# check for season effect
Nseason<-lmer(nvalue~season+(1|plot), data = dsgn%>%
                filter(nut=="PN")) 
n2<-signif(Anova(Nseason,Type="III"),2)
kbl(n2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the nitrogen content of seagrass in our study.",booktabs=TRUE)
```


```{r pseason,echo=FALSE,message=FALSE,warning=FALSE}
Pseason<-lmer(nvalue~season+(1|plot), data = dsgn%>%
                filter(nut=="PP")) 
p2<-signif(Anova(Pseason,Type="III"),2)
kbl(p2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the phosphorus content of seagrass in our study.",booktabs=TRUE)
```


```{r cseason,echo=FALSE,message=FALSE,warning=FALSE}
Cseason<-lmer(nvalue~season+(1|plot), data = dsgn%>%
                filter(nut=="PC")) 
c2<-signif(Anova(Cseason,Type="III"),2)
kbl(c2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the carbon content of seagrass in our study.",booktabs=TRUE)

```

We assumed that if there is an effect of sponges on seagrass nutrient content it would first appear closest to the sponge. 
So, we next examined the effect of treatment at each distance individually using a mixed effects model with seagrass nutrient content as the response, treatment and season as fixed effects, plot as a random effect, and the pre-experiment seagrass nutrient content as an offset.
For Nitrogen, there was a significant treatment effect at distance 0 m (Table \@ref(tab:nd0)), but treatment was no longer significant at 0.5, 1, and 2 m distances (Tables \@ref(tab:nd05),\@ref(tab:nd1),\@ref(tab:nd2) respectively)

```{r nd0,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
Nd0<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PN")%>%
            filter(dist==0)) 
nd0<-signif(Anova(Nd0,Type="III"),2)
kbl(nd0,caption="Results of the mixed effects model to determine if there is an effect of treatment on the nitrogen content of seagrass in our study at the center of each plot.",booktabs=TRUE)
```



```{r nd1,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
Nd1<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PN")%>%
            filter(dist==1)) 
nd1<-signif(Anova(Nd1,Type="III"),2)
kbl(nd1,caption="Results of the mixed effects model to determine if there is an effect of treatment on the nitrogen content of seagrass in our study at 1 m from the center of each plot.",booktabs=TRUE)
```

```{r nd2,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
Nd2<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PN")%>%
            filter(dist==2)) 
nd2<-signif(Anova(Nd2,Type="III"),2)
kbl(nd2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the nitrogen content of seagrass in our study at 2 m from the center of each plot..",booktabs=TRUE)
```

For phosphorus, there was no significant treatment effect at any distance (0 m, Table \@ref(tab:pd0); 0.5, Table \@ref(tab:pd05); 1 m, Table \@ref(tab:pd1); 2 m, Table \@ref(tab:pd2))

```{r pd0,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
pd0<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PP")%>%
            filter(dist==0)) 
pd0<-signif(Anova(pd0,Type="III"),2)
kbl(pd0,caption="Results of the mixed effects model to determine if there is an effect of treatment on the phosphorus content of seagrass in our study at the center of each plot.",booktabs=TRUE)
```

```{r pd05,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
pd05<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PP")%>%
            filter(dist==0.5)) 
pd05<-signif(Anova(pd05,Type="III"),2)
kbl(pd05,caption="Results of the mixed effects model to determine if there is an effect of treatment on the phosphorus content of seagrass in our study at 0.5 m from the center of each plot.",booktabs=TRUE)
```

```{r pd1,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
pd1<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PP")%>%
            filter(dist==1)) 
pd1<-signif(Anova(pd1,Type="III"),2)
kbl(pd1,caption="Results of the mixed effects model to determine if there is an effect of treatment on the phosphorus content of seagrass in our study at 1 m from the center of each plot.",booktabs=TRUE)
```

```{r pd2,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
pd2<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PP")%>%
            filter(dist==2)) 
pd2<-signif(Anova(pd2,Type="III"),2)
kbl(pd2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the phosphorus content of seagrass at 2 m in our study.",booktabs=TRUE)
```

For carbon, the pattern was similar to nitrogen. Treatment was significant at 0 m (Table \@ref(tab:cd0)) but not at any other distance (0.5 m, Table \@ref(tab:cd05); 1 m, Table \@ref(tab:cd1); and 2 m, Table \@ref(tab:cd2))

```{r cd0,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
cd0<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PC")%>%
            filter(dist==0)) 
cd0<-signif(Anova(cd0,Type="III"),2)
kbl(cd0,caption="Results of the mixed effects model to determine if there is an effect of treatment on the carbon content of seagrass in our study at the center of each plot.",booktabs=TRUE)
```

```{r cd05,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
cd05<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PC")%>%
            filter(dist==0.5)) 
cd05<-signif(Anova(cd05,Type="III"),2)
kbl(cd05,caption="Results of the mixed effects model to determine if there is an effect of treatment on the carbon content of seagrass in our study at 0.5 m from the center of each plot.",booktabs=TRUE)
```

```{r cd1,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
cd1<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PC")%>%
            filter(dist==1)) 
cd1<-signif(Anova(cd1,Type="III"),2)
kbl(cd1,caption="Results of the mixed effects model to determine if there is an effect of treatment on the carbon content of seagrass in our study at 1 m from the center of each plot.",booktabs=TRUE)
```

```{r cd2,echo=FALSE,message=FALSE,warning=FALSE}
# check for year effect
cd2<-lmer(nvalue~treatment+season+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PC")%>%
            filter(dist==2)) 
cd2<-signif(Anova(cd2,Type="III"),2)
kbl(cd2,caption="Results of the mixed effects model to determine if there is a seasonal pattern to the carbon content of seagrass in our study at 2 m from the center of each plot.",booktabs=TRUE)
```

Consequently, we chose to narrow our primary analysis to only the distance closest to the center of the plot. At this distance we only have a full dataset (i.e., no missing samples) for the winter of year 1 and the summer of year 2. As we already know there is a seasonal effect, and our preliminary data was from the summer, for the primary analysis presented in the main text, we chose to examine the change in nutrient levels in seagrass closest to the center of our plots one year into the experiment (i.e., the summer of year 2). 

For comparison purposes we show below what the results are if we include the winter sampling in the model so that there is an interaction between treatment and months into the experiment (Tables \@ref(tab:altnutn),\@ref(tab:altnutp),\@ref(tab:altnutc). The interpretation for the results does not change for nitrogen and phosphorus. However, due to our small sample size, by including the extra terms in this model we have lost the power to detect differences in carbon content. 


```{r altnutn,echo=FALSE,message=FALSE,warning=FALSE}
dn<-lmer(nvalue~treatment*mnths+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PN")%>%
           filter(mnths %in% c(5,12))%>%
            filter(dist==0)%>%
           mutate(mnths=factor(mnths),
                  treatment=factor(treatment),
                  treatment = relevel(treatment, ref = "Sponge"))) 

dns<-signif(summary(dn)$coefficient,2)

kbl(dns,caption="Results of the mixed effects model to determine if there is an effect of treatment on seagrass nitrogen content at 5 and 12 months into the experiment.",booktabs=TRUE)
```

```{r altnutp,echo=FALSE,message=FALSE,warning=FALSE}
dp<-lmer(nvalue~treatment*mnths+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PP")%>%
           filter(mnths %in% c(5,12))%>%
            filter(dist==0)%>%
           mutate(mnths=factor(mnths),
                  treatment=factor(treatment),
                  treatment = relevel(treatment, ref = "Sponge"))) 

dps<-signif(summary(dp)$coefficient,2)

kbl(dps,caption="Results of the mixed effects model to determine if there is an effect of treatment on seagrass phosphorus content at 5 and 12 months into the experiment.",booktabs=TRUE)
```

```{r altnutc,echo=FALSE,message=FALSE,warning=FALSE}
dc<-lmer(nvalue~treatment*mnths+(1|plot),
          offset=snv,
          data = dsgn%>%
            filter(nut=="PC")%>%
           filter(mnths %in% c(5,12))%>%
            filter(dist==0)%>%
           mutate(mnths=factor(mnths),
                  treatment=factor(treatment),
                  treatment = relevel(treatment, ref = "Sponge"))) 

dcs<-signif(summary(dc)$coefficient,2)

kbl(dcs,caption="Results of the mixed effects model to determine if there is an effect of treatment on seagrass carbon content at 5 and 12 months into the experiment.",booktabs=TRUE)
```