---
title: "Pascoe et al. Decadal change in seabird-driven isotopes on islands with differing invasion histories"
author: "Penny Pascoe"
date: "16/08/2024"
output: html_document
editor_options: 
  chunk_output_type: console
---

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


```{r}
library(tidyverse)
library(dplyr)
library(MuMIn)
library(multcomp)
library(sjPlot)
library(sjlabelled)
library(sjmisc)
library(multcomp)
library(corrplot)
library(RColorBrewer)
library(performance)
library(ggrepel)
library(lme4)
```

#Read in the data
Exact replicate data was matched up manualy in excel as it was simpler than formating all the location IDs from the two sampling events


#### Manuscript Results  ####

# Supplementary material Table 1
Mean d15n values for each ecosystem component on each island across the old and new sampling events
```{r}
df0<-read.csv(file="Marlborough_isotopes.csv", stringsAsFactors = TRUE)
df<-df0%>%filter(sample.type!="guano")

tableS1<-df%>%group_by(treatment, island,colony, spec.short)%>%summarise(
  no.samples=n(),
  m.new=mean(new, na.rm=TRUE), 
  sd.new=sd(new, na.rm=TRUE),
  m.old=mean(old, na.rm=TRUE), 
  sd.old=sd(old, na.rm=TRUE))

```


# Figure 3A - Mean change in δ15N values between 2006/07 and 2022 sampling events (± one standard deviation) for replicate soil, plant (Coprosma repens, C. robusta, Myrsine australis) and spider samples combined for each study island. 
```{r}
df.diff<-df%>%mutate(diff=new-old)

colorBlindGrey8   <- c("#999999","#0072B2", "#E69F00", "#56B4E9",  
                       "#F0E442", "#CC79A7","#009E73",  "#D55E00")


colorBlindGrey2   <- c("#CC79A7",  "#F0E442","#999999", "#56B4E9","#E69F00",  
                      "#0072B2", "#009E73",  "#D55E00")
colorblindcheck::palette_check(colorBlindGrey2, plot = TRUE)

df.diff$island<-factor(df.diff$island, levels=c("Kurupongi","Stephens","Mana","Maud","Nukuwaiata","Whakaterepapanui","Moutiti","Tawhitinui" ))


df.island<-df.diff%>%group_by(island, treatment)%>%summarise(
  n=n(),
  mean.diff=mean(diff, na.rm=TRUE), 
  sd.diff=sd(diff))%>%
  mutate(se.diff=sd.diff/sqrt(n))

  

  (Figure_2A_sd<-ggplot(df.island,aes(x=island, y=mean.diff, fill=island))+
  geom_hline(yintercept=0, color = "black", size=2)+
  scale_shape_manual(values=c(21, 24, 22,23))+    
  scale_colour_manual(values = c("Kurupongi" = "#F0E442",
                                 "Stephens" = "#CC79A7",
                                 "Mana" = "#999999",
                                 "Maud" = "#56B4E9",
                                 "Nukuwaiata"="#E69F00",
                                 "Whakaterepapanui"="#0072B2",
                                 "Moutiti" = "#009E73",
                                 "Tawhitinui" = "#D55E00" ))+ 
  scale_fill_manual(values = c("Kurupongi" = "#F0E442",
                                 "Stephens" = "#CC79A7",
                                 "Mana" = "#999999",
                                 "Maud" = "#56B4E9",
                                 "Nukuwaiata"="#E69F00",
                                 "Whakaterepapanui"="#0072B2",
                                 "Moutiti" = "#009E73",
                                 "Tawhitinui" = "#D55E00" ))+
  theme_bw(base_size = 20)+
  geom_errorbar(aes(ymin=mean.diff-sd.diff, ymax=mean.diff+sd.diff, col=island), lwd=1, width=.7)+
      geom_point(aes( shape=treatment), size=10)+
        ylab("Mean change in d15N") +
  xlab ("Island")+ theme(legend.position = "none")+
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.title.y=element_blank()))

```

# Figure 3B - Change in replicate soil, plant and spider samples displayed separately and broken up by island restoration treatment (facets), sample type (colours) and seabird colony status at the sampling location (shapes). 
```{r}
col2<-(c("#228B22","#FFFF00", "#00BFFF","#DC143C","#FFA500"))
colorblindcheck::palette_check(col2, plot = TRUE)

#add sample row id to each row for plot
df.diff2<-df.diff %>%filter(island!="Moutiti")%>% arrange(treatment,spec.short, colony) %>%group_by(treatment)%>% mutate(id2 = row_number())


 (Figure2_points<-ggplot(df.diff2,aes(x=id2, y=diff, fill=spec.short, shape=colony))+
  geom_hline(yintercept=0, color = "black", size=2)+
  geom_point(size=6)+
  scale_shape_manual(values=c(21, 24))+
  scale_fill_manual(values=col2)+
       ylab("change in d15N") +
  xlab ("Replicated Samples")+ 
    ylim(-6, 6)+
  theme_bw(base_size = 25)+
  facet_wrap(~treatment, scales = "free_x")+ theme(legend.position = "none"))

```

# Models
### Long-term δ15N change on individual islands

There was no significant change in average δ15N between the two sampling events for any of the islands (LRT, G1=1.29, P=0.255)
```{r}
df_long<-df%>%pivot_longer(cols=c('new', 'old',"old.sd", "new.sd"), values_to = 'd15n', names_to='round')
df_long_novic<-df_long%>%filter(island!="Moutiti"&round%in%c("new","old"))
df_long_novic$round<-factor(df_long_novic$round, levels = c("old", "new"))
fit.noVIc<-lmer(d15n~round+(1|island), data=df_long_novic)
drop1(fit.noVIc, test="Chisq")

plot(fit.noVIc)


df_long%>%filter(round%in%c("new", "old"))%>%dplyr::select(spec.short, treatment, d15n)%>%group_by(spec.short, treatment)%>%mutate(m=mean(d15n), sd=sd(d15n))

```

## Table 2
```{r}
summary(fit.noVIc)
r2_nakagawa(fit.noVIc)
```

Excluding Moutiti which increased by 8.5 ± 4.4 ‰ following the unexpected recolonisation of seabirds
```{r}
(df.island<-df.diff%>%group_by(island, treatment)%>%summarise(
  n=n(),
  mean.diff=mean(diff, na.rm=TRUE), 
  sd.diff=sd(diff))%>%
  mutate(se.diff=sd.diff/sqrt(n)))


```

# For individual replicated samples, there was also no effect of sample type (ANOVA, F4,57=0.44, P=0.78), colony status (ANOVA, F1,57=0.87, P=0.36) or island restoration treatment (ANOVA, F3,57=0.99, P=0.40) on the change in δ15N values between the sampling events (Table 3).

```{r}
 ## Compate the model with and without island
 d4<-df.diff%>%filter(island!="Moutiti")
  fit4.lm.island<-lm(diff~spec.short+colony+treatment+island, d4, na.action=na.fail)
 fit4.lm3.noisland<-lm(diff~spec.short+colony+treatment, d4)
 AIC(fit4.lm.island,fit4.lm3.noisland)


 #AIC says including island doesn't improve the model
 summary(fit4.lm3.noisland)
 anova(fit4.lm3.noisland)
```


# Restoration history variance
). Never-invaded islands showed the smallest variance in δ15N  (¬± 1.2 ‰ SD) and actively restored islands had the largest variance (± 2.7 ‰ SD)
```{r}
# Spread of points
df.diff%>%group_by(treatment)%>%summarise(
  n=n(),
  mean.diff=mean(diff, na.rm=TRUE), 
  sd.diff=sd(diff))%>%
  mutate(se.diff=sd.diff/sqrt(n))
```



## Guano variability
Is the guano signature similar between the islands?
Mana samples were unreliable as not dried for ages - exclude
```{r}
ggplot(df0%>%filter(spec.short=="guano" & island!="Mana"), aes(x=island, y=new, fill=island)) + 
geom_boxplot()+theme_classic(base_size = 20)+
    ylab({delta}^15*N~'\u2030') +
    xlab ("")
```


Excluding Mana samples, guano δ15N values did not statistically differ between the islands (ANOVA: F4,10=2.38, P = 0.12)
```{r}

df0$island<-as.factor(df0$island)
guano.df<-df0%>%filter(spec.short=="guano"& island!="Mana")

f.gg<-lm(new~island, data=guano.df, na.action = na.fail)

anova(f.gg)
summary(f.gg)

##mean guano = 11.4 +1 1.60SD
guano.df%>%group_by(sample.type)%>%summarise(m=mean(new), sd=sd(new))

```



# Space for time substitution analysis

#read in data
```{r}
sft.df<-read.csv(file="space_for_time.csv")

```
# Group by island/time since eradication/sample type
```{r}
ER_long<-sft.df%>%filter(treatment%in%c("ER", "AR"))%>%
    filter(sample.type%in%c("plant", "soil","spider"))%>%filter(island!="Maud")

INV<-sft.df%>%filter(treatment=="IN" )%>%filter(!(island=="Moutiti" & round  %in%c(2022,2006)))%>%
group_by( treatment,sample.type, spec_short )%>% 
  summarise( m.d15n=mean(d15n, na.rm=TRUE),
  sd.d15n=sd(d15n, na.rm=TRUE))%>%
    filter(sample.type%in%c("plant", "soil","spider"))

NI<-sft.df%>%filter(treatment=="NI" )%>%
group_by(treatment,sample.type, spec_short )%>% 
  summarise( m.d15n=mean(d15n, na.rm=TRUE),
  sd.d15n=sd(d15n, na.rm=TRUE))%>%
    filter(sample.type%in%c("plant", "soil","spider"))
```


# Figure 4 plots
```{r}

# Soil
(P_soil_long<-ggplot()+
#SD for INV_exaded and never INV_exaded islands
    geom_rect(aes(xmin=0, xmax=35, ymin=(INV%>%filter(spec_short == "soil" ))$m.d15n-(INV%>%filter(spec_short ==
    "soil" ))$sd.d15n, ymax = (INV%>%filter(spec_short == "soil" ))$m.d15n+(INV%>%filter(spec_short == "soil" ))$sd.d15n),
    fill='#D55E00', alpha=0.4)+
    geom_rect(aes(xmin=0, xmax=35, ymin=(NI%>%filter(spec_short == "soil" ))$m.d15n-(NI%>%filter(spec_short ==
    "soil" ))$sd.d15n,  ymax=(NI%>%filter(spec_short == "soil" ))$m.d15n+(NI%>%filter(spec_short == "soil"
    ))$sd.d15n),fill="#009E73",alpha=0.4)+
# trend line for ER islands
    # geom_smooth(data=ER_long%>%filter(spec_short == "soil" ), aes(x=t_erad, y=d15n), method="lm",se = TRUE,
    # fullrange=FALSE, col="black",fill="light grey", alpha=.9)+
# mean never INV_exaded and INV_exaded island lines
    geom_segment(data=INV%>%filter(spec_short == "soil" ), aes(x=0,xend=35,y=m.d15n,yend=m.d15n), col='#D55E00', lwd=1)+
    geom_segment(data=NI%>%filter(spec_short == "soil" ), aes(x=0,xend=35, y=m.d15n, yend=m.d15n), col="#009E73", lwd=1)+
# eradicated island points
    geom_point(data=ER_long%>%filter(spec_short == "soil" ), aes(x=t_erad, y=d15n, fill=island, size=5), shape=21)+
    scale_fill_manual(values = c("Mana" = "#999999",
                                "Nukuwaiata"="#E69F00",
                                "Whakaterepapanui"="#0072B2"))+ 
# appearance
    ggtitle("soil") + xlab(NULL)+ ylab(NULL) +theme_classic(base_size = 25) + xlim(0,35) + ylim(-10, 24)+
    theme(legend.position = "none"))


# Spider
(P_spider_long<-ggplot()+
#SD for INV_exaded and never INV_exaded islands
    geom_rect(aes(xmin=0, xmax=35, ymin=(INV%>%filter(spec_short == "spider" ))$m.d15n-(INV%>%filter(spec_short ==
    "spider" ))$sd.d15n, ymax = (INV%>%filter(spec_short == "spider" ))$m.d15n+(INV%>%filter(spec_short == "spider" ))$sd.d15n),
    fill='#D55E00', alpha=0.4)+
    geom_rect(aes(xmin=0, xmax=35, ymin=(NI%>%filter(spec_short == "spider" ))$m.d15n-(NI%>%filter(spec_short ==
    "spider" ))$sd.d15n,  ymax=(NI%>%filter(spec_short == "spider" ))$m.d15n+(NI%>%filter(spec_short == "spider"
    ))$sd.d15n),fill="#009E73",alpha=0.4)+
# trend line for ER islands
    # geom_smooth(data=ER_long%>%filter(spec_short == "spider" ), aes(x=t_erad, y=d15n), method="lm",se = TRUE,
    # fullrange=FALSE, col="black",fill="light grey", alpha=.9)+
# mean never INV_exaded and INV_exaded island lines
    geom_segment(data=INV%>%filter(spec_short == "spider" ), aes(x=0,xend=35,y=m.d15n,yend=m.d15n), col='#D55E00', lwd=1)+
    geom_segment(data=NI%>%filter(spec_short == "spider" ), aes(x=0,xend=35, y=m.d15n, yend=m.d15n), col="#009E73", lwd=1)+
# eradicated island points
    geom_point(data=ER_long%>%filter(spec_short == "spider" ), aes(x=t_erad, y=d15n, fill=island, size=5), shape=21)+
    scale_fill_manual(values = c("Mana" = "#999999",
                                "Nukuwaiata"="#E69F00",
                                "Whakaterepapanui"="#0072B2"))+ 
# appearance
    ggtitle("spider") + xlab(NULL)+ ylab(NULL) +theme_classic(base_size = 25) + xlim(0,35) + ylim(-10, 24)+
    theme(legend.position = "none"))



# Plant
(P_Cop.rep_long<-ggplot()+
#SD for INV_exaded and never INV_exaded islands
    geom_rect(aes(xmin=0, xmax=35, ymin=(INV%>%filter(spec_short == "Cop.rep" ))$m.d15n-(INV%>%filter(spec_short ==
    "Cop.rep" ))$sd.d15n, ymax = (INV%>%filter(spec_short == "Cop.rep" ))$m.d15n+(INV%>%filter(spec_short == "Cop.rep" ))$sd.d15n),
    fill='#D55E00', alpha=0.4)+
    geom_rect(aes(xmin=0, xmax=35, ymin=(NI%>%filter(spec_short == "Cop.rep" ))$m.d15n-(NI%>%filter(spec_short ==
    "Cop.rep" ))$sd.d15n,  ymax=(NI%>%filter(spec_short == "Cop.rep" ))$m.d15n+(NI%>%filter(spec_short == "Cop.rep"
    ))$sd.d15n),fill="#009E73",alpha=0.4)+
# trend line for ER islands
    # geom_smooth(data=ER_long%>%filter(spec_short == "Cop.rep" ), aes(x=t_erad, y=d15n), method="lm",se = TRUE,
    # fullrange=FALSE, col="black",fill="light grey", alpha=.9)+
# mean never INV_exaded and INV_exaded island lines
    geom_segment(data=INV%>%filter(spec_short == "Cop.rep" ), aes(x=0,xend=35,y=m.d15n,yend=m.d15n), col='#D55E00', lwd=1)+
    geom_segment(data=NI%>%filter(spec_short == "Cop.rep" ), aes(x=0,xend=35, y=m.d15n, yend=m.d15n), col="#009E73", lwd=1)+
# eradicated island points
    geom_point(data=ER_long%>%filter(spec_short == "Cop.rep" ), aes(x=t_erad, y=d15n, fill=island, size=5), shape=21)+
    scale_fill_manual(values = c("Mana" = "#999999",
                                "Nukuwaiata"="#E69F00",
                                "Whakaterepapanui"="#0072B2"))+ 
# appearance
    ggtitle("Cop.rep") + xlab(NULL)+ ylab(NULL) +theme_classic(base_size = 25) + xlim(0,35) + ylim(-10, 24)+
    theme(legend.position = "none"))


#### model the trends of a linear model through the eradicated samples
#Raw data models

#plant
fit_c.rep_long<-lm(d15n~t_erad, data=ER_long%>%filter(spec_short == "Cop.rep" ))
drop1(fit_c.rep_long, test="Chisq")
summary(fit_c.rep_long) # no significant trend

#soil
fit_soil_long<-lm(d15n~t_erad, data=ER_long%>%filter(spec_short == "soil" ))
drop1(fit_soil_long, test="Chisq")
summary(fit_soil_long)# no significant trend

#spider
fit_spider_long<-lm(d15n~t_erad, data=ER_long%>%filter(spec_short == "spider" ))
drop1(fit_spider_long, test="Chisq")
summary(fit_spider_long)# no significant trend

```
