## Pascoe et al (2022). Temporal and spatial variability in stable isotope values on seabird islands: what, where and when to sample 

## Libraries
library(tidyverse)
library(lme4)
library(multcomp)
library(corrplot)
library(MuMIn)
library(performance)
library(RColorBrewer)
library(ggprism)

## Read in data

#### Ecosystem component variation

## Do plant species differ
P1<-lmer(dn15~ spec.short +(1|island)+(1|seabirds_present)+(1|round), data=iso_data%>%filter(sample.type=="plant"))
drop1(P1, test="Chisq")
summary(P1)
# Yes - Which ones?
summary(glht(P1, linfct = mcp(spec.short = "Tukey")))




## Do spider family, age, sex differ?
S1<-lmer(dn15 ~ family+sex+age +(1|island)+(1|seabirds_present)+(1|round), data=iso_data%>%filter(sample.type=="spider"& sex!="unk" & age!="unk"))
drop1(S1, test="Chisq")
#only family has an effect so drop the others
S1.1<-lmer(dn15 ~ family+(1|island)+(1|seabirds_present)+(1|round), data=iso_data%>%filter(sample.type=="spider"))
summary(glht(S1.1, linfct = mcp(family = "Tukey")))
#Lycosidae and Miturgidae are similar



## does guano vary between islands or seasons?
gg<-iso_data%>%filter(sample.type=="guano")
f.g<-lm(dn15~island+round+seabirds_present, data=gg)
anova(f.g)
# No, it doesn't


## Select just the spiders and plants that are similar
rr<-iso_data%>%filter(spec.short %in%c("soil","Pte.esc", "Poa.poi", "Lycosidae", "Miturgidae","guano"))



## Mean and range of d15n values for each sample type
rr%>%group_by(spec.short)%>%summarise(min=min(dn15), max=max(dn15), mean=mean(dn15), sd=sd(dn15))



## correlation between d15n in the ecosystem components collected at a site

corr<-rr%>%group_by(sample.type, island, seabirds_present, round)%>%summarise(dn15=mean(dn15, na.rm=TRUE))%>%pivot_wider(names_from = sample.type, values_from=dn15)

#soil - plant
cor.test(corr$soil, corr$plant, method = "pearson")

#Soil - spider
cor.test(corr$soil, corr$spider, method = "pearson")

#Plant spider
cor.test(corr$plant, corr$spider,method = "pearson")

#guano-plant
cor.test(corr$plant, corr$guano, method = "pearson")

#guano-soil
cor.test(corr$soil, corr$guano,method = "pearson")

#guano-spider
cor.test(corr$spider, corr$guano, method = "pearson")




### Spatial and Temporal variation

#Remove non-pitfall locations and Tasman
r1<-rr %>% filter(island!="Tasman", point.type!="oppertunistic")

# Set the order of sample type
r1$sample.type<-  factor(r1$sample.type, levels=c("soil","plant","spider","guano"))
r1$sample.type<- factor(r1$sample.type)
r1$seabirds_present <- factor(r1$seabirds_present, levels = c("no", "close", "colony"))

# Remove guano as it is consistent across the islands and seasons
r1.ng<- r1%>%filter(sample.type!="guano")%>%droplevels()

# Are any of the variables highly correlated?
corrplot(cor(r1.ng%>%dplyr::select(round, tot_rain_30,max_rain_30)),
         method = "number",
         type = "upper") # No, they look fine


# Fit full model
F1 <- lmer (dn15~
              seabirds_present + round + tot_rain_30 + max_rain_30 + sample.type + 
              (seabirds_present*sample.type) + (sample.type*round)+ 
              (1|island), 
            data = r1.ng, na.action=na.fail)

dredge(F1, rank="AIC") ## Temporal terms are not important predictors

# plot change in d15n with season
ggplot(r1, aes(x=as.character(round), y=dn15))+geom_boxplot()+facet_wrap(~sample.type)+theme_classic() # no change


#Select best model
F.best <- lmer(dn15~
                 seabirds_present+sample.type+
                 (seabirds_present*sample.type)+
                 (1|island), 
               data = r1.ng)
summary(F.best)
r2_nakagawa(F.best)

#Plot the raw and predicted results (Paper Fig 3)
predicted_values<-modelr::data_grid(r1.ng, dn15, seabirds_present, sample.type, island)%>%
  modelr::add_predictions(F.best)

mean_g<-r1%>%filter(sample.type=="guano")%>%summarise(mean_guano=mean(dn15), sd_guano = sd(dn15))


myColors <- brewer.pal(3,"Dark2")
names(myColors) <- levels(r1.ng$sample.type)
colScale <- scale_colour_manual(name = "sample.type",values = myColors)

xlabs <- c("Outside", "Near", "Colony")


(Figure3<-predicted_values %>% 
    ggplot(aes(seabirds_present, pred,  fill=sample.type))+
    annotate(geom="rect",xmin=-Inf, xmax=Inf, ymin=(mean_g$mean_guano - mean_g$sd_guano), ymax = (mean_g$mean_guano + mean_g$sd_guano), fill="gray26",alpha=0.2)+
    geom_hline(yintercept=mean_g$mean_guano, color = "gray26", size=1)+
    geom_jitter(data = r1.ng, aes(seabirds_present, dn15, col = sample.type, alpha=0.8),width = 0.2)+
    geom_boxplot(aes(middle = mean(dn15)))+ #Box plots are mean predicted values
    theme_classic(base_size = 20) +
    ylab({delta}^15*N~'\u2030') +
    xlab ("Sampling Location")+
    guides(alpha="none", col = "none")+guides(fill=guide_legend(title="Ecosystem Component"))+ scale_x_discrete(labels= xlabs)+
    scale_color_manual(values = c("plant" = "yellowgreen",
                                  "soil"="cyan3",
                                  "spider" ="mediumpurple3"))+
    scale_fill_manual(values=c("cyan3","yellowgreen",  "mediumpurple2"))+  scale_y_continuous(guide = "prism_minor", 
                                                                                              limits = c(-15, 30),
                                                                                              expand = c(0, 0),
                                                                                              minor_breaks = seq(-14, 30, 2))+ theme(legend.title = element_blank(),
                                                                                                                                     prism.ticks.length.y = unit(2.5, "pt")))

## Spatial variation with transect data
# Including Tasman Island, exclude Maatsuyker Island and use a combination of pitfall and transect points

r_dist<-rr%>%
  filter(!is.na(dist.col) & point.type!="opportunistic"  & sample.type!="guano" )%>%
  droplevels()

cor.test(r_dist$dist.coast, r_dist$DEM, method = "pearson") # elevation and distance to coast are correlated

#Fit full model
F.dist <- lmer(dn15~
                 log(dist.cat+1)*dist.coast+
                 (1|island)+(1|sample.type), data=r_dist, na.action=na.fail)

drop1(F.dist, test="Chisq")

dredge(F.dist) #Just distance to colony is important

F.dist.best<-lmer(dn15~log(dist.cat+1)+(1|island)+(1|sample.type), data=r_dist, na.action=na.fail)
summary(F.dist.best)

## plot the model output (Fiure 4)

# Make a grid with observed and predicted values for d15n
grid <- r_dist%>%   modelr::data_grid( dist.cat, sample.type, island)
grid
grid_pred <- grid %>%
  modelr::add_predictions(F.dist.best)
grid_pred


(Figure4<-ggplot(data = r_dist) +
    geom_point(aes(x = dist.cat, y = dn15, col=sample.type), pch=19) +
    # Overlay grid predictions on ggplot
    geom_smooth(aes(x = (dist.cat+1), y = pred, col=sample.type), data = grid_pred, size = 1, method="glm", formula = 'y~log(x)', se=TRUE, fullrange=FALSE)+
    theme_classic(base_size = 20)+
    ylab(~{delta}^15*N~'\u2030') +
    xlab ("Distance to Colony (m)")+ 
    labs(col = "Ecosystem Component")+
    scale_color_manual(values = c("plant" = "yellowgreen",
                                  "soil"="cyan3",
                                  "spider" ="mediumpurple3"))+
    scale_fill_manual(values=c("plant" = "yellowgreen",
                               "soil"="cyan3",
                               "spider" ="mediumpurple3"))+  scale_y_continuous(guide = "prism_minor", 
                                                                                limits = c(-15, 30),
                                                                                expand = c(0, 0),
                                                                                minor_breaks = seq(-14, 30, 2))+ theme(legend.title = element_blank(),
                                                                                                                       prism.ticks.length.y = unit(2.5, "pt")))



##Seabird colony density
r_col<-rr %>% filter(seabirds_present=="colony" & !is.na(burrows), point.type=="pitfall")
C1.1<-lmer(dn15~ burrows_m2+(1|island)+(burrows_m2|sample.type), data=r_col, na.action=na.fail)
r2_nakagawa(C1.1)
summary(C1.1)

## Plot the model output - Figure 5
# Make a grid with observed and predicted values for d15n
grid_col <- r_col%>%   modelr::data_grid( burrows_m2, sample.type, island)
grid_col
grid_col_pred <- grid_col %>%
  modelr::add_predictions(C1.1)
grid_col_pred

#raw data with the modeled line based on the predicted values from the model
(Figure5<-ggplot(data = r_col) +
    geom_point(aes(x = burrows_m2, y = dn15, col=sample.type), pch=19) +
    # Overlay grid_col predictions on ggplot
    geom_smooth(aes(x = burrows_m2, y = pred, col=sample.type), data = grid_col_pred, size = 1, method="glm", formula = 'y~x', se=TRUE, fullrange=FALSE)+
    theme_classic(base_size = 20)+
    ylab(~{delta}^15*N~'\u2030') +
    xlab ("Number of burrow entrances (1m"^2~')')+ 
    labs(col = "Ecosystem Component")+
    scale_color_manual(values = c("plant" = "yellowgreen",
                                  "soil"="cyan3",
                                  "spider" ="mediumpurple3",
                                  "guano" ="gray26"))+
    scale_fill_manual(values=c("plant" = "yellowgreen",
                               "soil"="cyan3",
                               "spider" ="mediumpurple3",
                               "guano" ="gray26"))+  scale_y_continuous(guide = "prism_minor", 
                                                                        limits = c(0, 30),
                                                                        expand = c(0, 0),
                                                                        minor_breaks = seq(2, 30, 2))+ theme(legend.title = element_blank(),
                                                                                                             prism.ticks.length.y = unit(2.5, "pt")))