#################################################################################
############### ANALYSES FOR BIOMASS, BREEDING SUCCESS AND DIET  ################
#################################################################################

###### Packages ####
library(mgcv)
library(mgcViz)
library(MuMIn)
library(visreg)
library(psych)
library(ggplot2)
library(jtools)
library(multcomp)
library(car)
library(plyr)
library(corrplot)
library(gridExtra)



###################################################################################
################ 1. GAM analyses for biomass and breeding success  ################

##### Files and data #####
setwd("...")
data <- read.csv("Pelletier2021_biomass_vs_BrS.csv", sep=";")
as.data.frame(data)
str(data)
names(data)
colnames(data)[1] <- "Yr"

subset19792008 <- data[1:8,]
subset20092019 <- data[9:19,]

describe(subset19792008)
describe(subset20092019)


##### Plots #####
plot(data$BrS ~ data$MackTSB)
plot(data$BrS ~ data$HerrFall4TTOT2020)
plot(data$BrS ~ data$HerrSpring4TTOT2020)


##### GAM analyses #####
gam_mod_herr <- gam(BrS ~ s(HerrSpring4TTOT2020) + s(HerrFall4TTOT2020), data = data, method = "REML")
gam_mod_herrSmack <- gam(BrS ~ s(HerrSpring4TTOT2020) + s(MackTSB), data = data, method = "REML")
gam_mod_herrFmack <- gam(BrS ~ s(HerrFall4TTOT2020) + s(MackTSB), data = data, method = "REML")
gam_mod_mack <- gam(BrS ~ s(MackTSB), data = data, method = "REML")
gam_mod_herrF <- gam(BrS ~ s(HerrFall4TTOT2020), data = data, method = "REML")
gam_mod_herrS <- gam(BrS ~ s(HerrSpring4TTOT2020), data = data, method = "REML")


##### Model selection #####
model.sel(gam_mod_herrSmack, gam_mod_herrFmack, gam_mod_herr, gam_mod_herrS, gam_mod_herrF, gam_mod_mack,
          rank = AICc, extra = c(BIC, AICc))

anova(gam_mod_herrFmack, gam_mod_mack, test = "Chisq")

gam.check(gam_mod_mack)


##### Model summary #####
summary(gam_mod_mack)
mean(summary(gam_mod_mack)$s.table[,4])

summary(gam_mod_herrF)
summary(gam_mod_herrS)
summary(gam_mod_herrFmack)
summary(gam_mod_herrSmack)
summary(gam_mod_herr)


##### Plots #####
gam_mod_mack_plot <- ggplot(data, aes(x=MackTSB/1000, y=BrS/100)) +
                      stat_smooth(method = "gam", color = "black", fill="grey60") +
                      geom_point() +
                      ylim(0, 1) +
                      xlab("Spawning-stock biomass of \n Atlantic mackerel (x 1000 t)") + ylab(expression("Breeding success "~(chick.yr^-1))) +
                      theme_apa()
gam_mod_mack_plot

ggsave("Fig_GAM_BreedingSuccess_vs_MackSSB.png", plot = gam_mod_mack_plot,
       path = "...",
       width = 3.5, height = 3.5, units = "in", dpi = 300)


par(mfrow = c(2, 2))
visreg(gam_mod_mack, "MackTSB", gg=F, type="conditional", overlay=T, partial=T, rug=F,
       ylab="Breeding success of northern gannets (%)", xlab="Mackerel abundance (t)", alpha = 0.05,
       line=list(col="black"),
       fill=list(col="#66ccff"),
       points=list(cex=1, pch=16, col="black"))

visreg(gam_mod_herrF, "HerrFall4TTOT2020", gg=F, type="conditional", overlay=T, partial=T, rug=F,
       ylab="Breeding success of northern gannets (%)", xlab="Fall herring abundance (t)", alpha = 0.05,
       line=list(col="black"),
       fill=list(col="#66CCFF"),
       points=list(cex=1, pch=16, col="black"))

visreg(gam_mod_herrS, "HerrSpring4TTOT2020", gg=F, type="conditional", overlay=T, partial=T, rug=F,
       ylab="Breeding success of northern gannets (%)", xlab="Spring herring abundance (t)", alpha = 0.05,
       line=list(col="black"),
       fill=list(col="#66CCFF"),
       points=list(cex=1, pch=16, col="black"))

dev.off()


####################################################################
################       2. Biomass comparison        ################

##### Files and data #####
setwd("...")
biomass <- read.csv("Pelletier2021_biomass.csv", sep=";")
str(biomass)

biomass$Species <- as.factor(biomass$Species)
biomass$Period <- as.factor(biomass$Period)


##### FOR biomass VS SPECIES #####

# Compute the analysis of variance
res.aov1 <- aov(biomass ~ Species, data = biomass)

# Summary of the analysis
summary(res.aov1)

# Pair-wise comparison
summary(glht(res.aov1, linfct = mcp(Species = "Tukey")))

# Homogeneity of variances
plot(res.aov1, 1)
leveneTest(biomass ~ Species, data = biomass) #inequality

# Normality assumption
plot(res.aov1, 2) #abnormal

# Extract the residuals
aov_residuals1 <- residuals(object = res.aov )

# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals1 ) # abnormal

# Non-parametric alternative to one-way ANOVA test
np.aov1 <- kruskal.test(biomass ~ Species, data = biomass)
print(np.aov1)
summary(glht(data = np.aov1, linfct = mcp(Species = "Tukey")))
pairwise.wilcox.test(biomass$biomass, biomass$Species,
                     p.adjust.method = "BH")

ddply(biomass, .(Species), summarize,  biomass=mean(biomass))


##### FOR biomass VS PERIOD ##### 

# Compute the analysis of variance
res.aov2 <- aov(biomass ~ Period, data = biomass)

# Summary of the analysis
summary(res.aov2)

# Pair-wise comparison
summary(glht(res.aov2, linfct = mcp(Period = "Tukey")))

# Homogeneity of variances
plot(res.aov2, 1)
leveneTest(biomass ~ Period, data = biomass) #inequality

# Normality assumption
plot(res.aov2, 2) #abnormal

# Extract the residuals
aov_residuals2 <- residuals(object = res.aov2)

# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals2 ) # abnormal

# Non-parametric alternative to one-way ANOVA test
np.aov2 <- kruskal.test(biomass ~ Period, data = biomass)
print(np.aov2)
summary(glht(data = np.aov2, linfct = mcp(Period = "Tukey")))
pairwise.wilcox.test(biomass$biomass, biomass$Period,
                     p.adjust.method = "BH")

# Non-parametric alternative to one-way ANOVA test for MACK biomass
np.aovMACK <- kruskal.test(biomass ~ Period, data = biomass[1:19,])
print(np.aovMACK)

# Non-parametric alternative to one-way ANOVA test for HERRF biomass
np.aovHERRF <- kruskal.test(biomass ~ Period, data = biomass[20:38,])
print(np.aovHERRF)

# Non-parametric alternative to one-way ANOVA test for HERRS biomass
np.aovHERRS <- kruskal.test(biomass ~ Period, data = biomass[39:57,])
print(np.aovHERRS)

# Non-parametric alternative to one-way ANOVA test for PERIOD1 (1979-2008)
np.aovPeriod1 <- kruskal.test(biomass ~ Species, data = biomass[biomass$Period==1,])
print(np.aovPeriod1)
pairwise.wilcox.test(biomass[biomass$Period==1,]$biomass, biomass[biomass$Period==1,]$Species,
                     p.adjust.method = "BH")

# Non-parametric alternative to one-way ANOVA test for PERIOD2 (2009-2019)
np.aovPeriod2 <- kruskal.test(biomass ~ Species, data = biomass[biomass$Period==2,])
print(np.aovPeriod2)
pairwise.wilcox.test(biomass[biomass$Period==2,]$biomass, biomass[biomass$Period==2,]$Species,
                     p.adjust.method = "BH")


###################################################################################
################ 3. Differences of breeding success between period ################

##### Files and data #####
BrS <- read.csv("Pelletier2021_BrS.csv", sep=";")
str(BrS)

BrS$Period <- as.factor(BrS$Period)

ddply(BrS, .(Period), summarize,  BrS=mean(BrS))

# Non-parametric alternative to one-way ANOVA test for BrS (between period)
np.aovPeriodBrS <- kruskal.test(BrS ~ Period, data = BrS)
print(np.aovPeriodBrS)


###################################################################################
################ 4. Relationships between diet and breeding success ###############

##### Files and data #####
dietyr <- read.csv("Pelletieretal2021_BrSvsDiet.csv", sep=";")

as.data.frame(dietyr)
str(dietyr)
names(dietyr)

describe(dietyr)


##### Coefficient of variation (CV)
sapply(dietyr, function(x) sd(x) / mean(x) * 100)


#### Mackerel
set.seed(123456)
btcorMACKdiet <- boot(dietyr,
               statistic = function(data, i)
               {cor(data[i, "Diet_pmack"], data[i, "BrS"], method='pearson')},
               R = 10000)
btcorMACKdiet
boot.pval(btcorMACKdiet, type = "norm", theta_null = 0)
median(btcorMACKdiet$t)
boot.ci(btcorMACKdiet, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI.
plot(density(btcorMACKdiet$t))
abline(v = mean(btcorMACKdiet$t0), lty = "dashed", col = "grey60")


#### Herring
btcorHERRdiet <- boot(dietyr,
                      statistic = function(data, i)
                      {cor(data[i, "Diet_pherr"], data[i, "BrS"], method='pearson')},
                      R = 10000)
btcorHERRdiet
boot.pval(btcorHERRdiet, type = "norm", theta_null = 0)
median(btcorHERRdiet$t)
boot.ci(btcorHERRdiet, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI.
plot(density(btcorHERRdiet$t))
abline(v = mean(btcorHERRdiet$t0), lty = "dashed", col = "grey60")


#### Shannon index
btcorDIVdiet <- boot(dietyr,
                      statistic = function(data, i)
                      {cor(data[i, "Diet_diversity"], data[i, "BrS"], method='pearson')},
                      R = 10000)
boot.pval(btcorDIVdiet, type = "norm", theta_null = 0)
median(btcorDIVdiet$t)
boot.ci(btcorDIVdiet, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI.
plot(density(btcorDIVdiet$t))
abline(v = mean(btcorDIVdiet$t0), lty = "dashed", col = "grey60")


# Sand lance
btcorSANDdiet <- boot(dietyr,
                     statistic = function(data, i)
                     {cor(data[i, "Diet_psand"], data[i, "BrS"], method='pearson')},
                     R = 10000)
boot.pval(btcorSANDdiet, type = "norm", theta_null = 0)
boot.pval(btcorSANDdiet, type = "basic", theta_null = 0)
boot.pval(btcorSANDdiet, type = "perc", theta_null = 0)
median(btcorSANDdiet$t)
boot.ci(btcorDIVdiet, type = c("norm", "basic", "perc", "bca")) #bootstrapped CI.
plot(density(btcorSANDdiet$t))
abline(v = mean(btcorSANDdiet$t0), lty = "dashed", col = "grey60")


#### Plots: linear trend + confidence interval
BrSuccvsDietDiv <- ggplot(dietyr, aes(x=Diet_diversity, y=BrS)) +
                   geom_smooth(method=lm , color="black", fill="grey", se=TRUE) +
                    geom_point() +
                    ylim(0, 1) +
                    xlab("Diet Diversity (Shannon index)") +
                    theme_apa() +
                    theme(axis.title.y=element_blank(), axis.text.y = element_blank())
BrSuccvsDietDiv

BrSuccvsDietSand <- ggplot(dietyr, aes(x=Diet_psand, y=BrS)) +
  geom_smooth(method=lm , color="black", fill="grey", se=TRUE) +
  geom_point() +
  ylim(0, 1) +
  xlab("Proportion of sandlance in diet") + ylab(expression("Breeding success "~(chick.yr^-1))) +
  theme_apa()
BrSuccvsDietSand

BrSuccvsDietMack <- ggplot(dietyr, aes(x=Diet_pmack, y=BrS)) +
  geom_smooth(method=lm , color="black", fill="grey", se=TRUE) +
  geom_point() +
  ylim(0, 1) +
  xlab("Proportion of Atlantic mackerel in diet") + ylab(expression("Breeding success "~(chick.yr^-1))) +
  theme_apa()
BrSuccvsDietMack

BrSuccvsDietHerr <- ggplot(dietyr, aes(x=Diet_pherr, y=BrS)) +
  geom_smooth(method=lm , color="black", fill="grey", se=TRUE) +
  geom_point() +
  ylim(0, 1) +
  xlab("Proportion of Atlantic herring in diet") +
  theme_apa() +
  theme(axis.title.y=element_blank(), axis.text.y = element_blank())
BrSuccvsDietHerr

BrSuccvsDietCape <- ggplot(dietyr, aes(x=Diet_pcap, y=BrS)) +
  geom_smooth(method=lm , color="black", fill="grey", se=TRUE) +
  geom_point() +
  xlab("Proportion of capelin in diet") + ylab(expression("Breeding success "~(chick.yr^-1))) +
  theme_apa()
BrSuccvsDietCape

Diet_plots <- grid.arrange(BrSuccvsDietMack, BrSuccvsDietHerr,
            BrSuccvsDietSand, BrSuccvsDietDiv, ncol=2, nrow = 2, widths = c(1.2,1))

ggsave("Fig_BreedingSuccess_vs_diet.png", plot = Diet_plots,
       path = "...",
       width = 7, height = 7, units = "in", dpi = 300)



#############################################################################
################    5. Difference of diet between species    ################

DIET <- read.csv("Pelletier2021_diet.csv", sep=";")
str(DIET)

DIET$Species <- as.factor(DIET$Species)
DIET$Year <- as.factor(DIET$Year)

ddply(DIET, .(Species), summarize,  PropDiet=mean(PropDiet))

# Compute the analysis of variance
res.aov5 <- aov(PropDiet ~ Species, data = DIET)

# Summary of the analysis
summary(res.aov5)

# Pair-wise comparison
summary(glht(res.aov5, linfct = mcp(Species = "Tukey")))

# Homogeneity of variances
plot(res.aov5, 1)
leveneTest(PropDiet ~ Species, data = DIET) #equality

# Normality assumption
plot(res.aov5, 2) #abnormal

# Extract the residuals
aov_residuals5 <- residuals(object = res.aov5)

# Run Shapiro-Wilk test
shapiro.test(x = aov_residuals5) # abnormal

# Non-parametric alternative to one-way ANOVA test for BrS (between period)
np.aovDIET <- kruskal.test(PropDiet ~ Species, data = DIET)
print(np.aovDIET)

summary(glht(data = np.aovDIET, linfct = mcp(Species = "Tukey")))
pairwise.wilcox.test(DIET$PropDiet, DIET$Species,
                     p.adjust.method = "BH",)