---
title: Variability in krill biomass combined with fixed catch limits affected Pygoscelis
  penguin breeding success in the Western Antarctic Peninsula
author: "Lucas Krüger, Francisco Santa Cruz, Mauricio Mardones, César Cárdenas"
date: "`r Sys.Date()`"
output: html_document
---
#Load necessary packages

```{r, echo=FALSE}
library(plyr)
library(dplyr) # for summarizing data
library(reshape2)
library(ggplot2)  # plotting 
library(sjPlot)
library(patchwork)
library(lubridate) # to deal with time data
library(udunits2) # to deal with geographical data
library(raster)
library(MCMCglmm) #Bayesian generalized linear mixed models
```


#theme for the plots:

```{r}

th<- theme(axis.text=element_text(size=10, face="bold",colour="grey30"),
           axis.title=element_text(size=12,face="bold"),
           legend.text = element_text(size=12),
           panel.grid.major = element_blank(),
           panel.grid.minor = element_blank(),
           panel.spacing.x = unit(1, "lines"))
```

#change working directory

```{r}

setwd("D:/BACKUP_Notebook_06_2022/KidsArentAlright/")

```

#load krill acoustic data from "WG-ASAM Results from the WG-ASAM Intersessional e-Group on Krill Biomass Estimates from Acoustic Surveys; 2021;"

```{r}


krill<-read.csv("D:/BACKUP_Notebook_06_2022/KidsArentAlright/Data/krill_biomass_acoustic_survey.csv") # data from WG-ASAM e-group (2021)


krill2<-subset(krill,
               Strata=="Bransfield"|Strata=="West") ### subset the strata 
                                                    ### where fishing has been consistent

```


#load and process fishing data 

```{r}
kf0<-read.csv("Data/C1_Sect.csv") ### fishing data from CCAMLR secretariat  

kf0$CatchT<-kf0$KrillCaughtKG/1000  ### catch in tones

kf0$ts<- as.POSIXct(strptime(kf0$Day,format="%m/%d/%Y" ,tz="GMT")) 

kf0$month<-month(kf0$ts)

kf0$FSeason<-ifelse(kf0$month>11,
            kf0$CalendarYear,kf0$CalendarYear-1) ### fishing season starts in december
                                                 ### january up to november next year
                                                 ### is assigned to the season of the previous year
                                                 ### krill season however, places december to the
                                                 ### following year, 
                                                 ### therefore fishing season 2018=krill season 2019
                                                 
kf<-kf0

strata<-shapefile("Data/Shapefiles/AMLR_STRATA_.shp")  ### from CCAMLR as well

crs=CRS( "+proj=longlat +datum=WGS84 +no_defs")

kcoords<-coordinates(data.frame(x=kf$FishingLon,y=kf$FishingLat))

kf.spdf<-SpatialPointsDataFrame(kcoords,kf,proj4string = crs)

kf.strata<-over(kf.spdf,strata)

kf.spdf$Strata<-kf.strata$Strata


kf.spdf<-subset(kf.spdf,FishingLon<(-50))  # eliminate the South Orkney Data 

```

#### bring krill and fishery together

```{r}

kf.spdf$Year<-kf.spdf$CalendarYear

### months with krill estimation

kfMo<-subset(kf.spdf,month=="12"|month=="1"|month=="2"|month=="3") 

fcr2<-ddply((data.frame(kfMo)), c("Year","Strata"), summarise,
           CatchS=sum(CatchT))  ### total catch for each year and strata

summary(fcr2) ### CatchS is catch in Ton

head(krill) ### amlrbiomass is in ton

kr<-data.frame(krill[1],krill[2],krill[14]) # keep only needed variables

kmfcr<-merge(kr,fcr2,by=c("Strata","Year"),all=T)  # merge krill and fishery data

kmfcr2<-subset(kmfcr,Year<2019 & Year>1994) 

# check 1st quantiles to fill missing values
summary(kmfcr2$amlrbiomass[kmfcr2$Strata=="Bransfield"]) 
summary(kmfcr2$amlrbiomass[kmfcr2$Strata=="West"])

kmfcr2$amlrbiomass[is.na(kmfcr2$amlrbiomass) & kmfcr2$Strata=="Bransfield"]<-392888# 1st qu
kmfcr2$amlrbiomass[is.na(kmfcr2$amlrbiomass) & kmfcr2$Strata=="West"]<-752566# 1st qu

# substitute missing values by the mean
summary(kmfcr2$CatchS[kmfcr2$Strata=="West"])
summary(kmfcr2$CatchS[kmfcr2$Strata=="Bransfield"])

kmfcr2$CatchS[is.na(kmfcr2$CatchS) & kmfcr2$Strata=="West"]<-8163 # mean
kmfcr2$CatchS[is.na(kmfcr2$CatchS) &  kmfcr2$Strata=="Bransfield"]<-8163 # mean

kmfcr2$Strata<-factor(kmfcr2$Strata,levels=c("West","Bransfield"))


kmfcr2$CatchK<-kmfcr2$CatchS/kmfcr2$amlrbiomass # harvest rate 

plot2a<-ggplot(subset(kmfcr2,Strata=="Bransfield"|Strata=="West"),
               aes(Year,(amlrbiomass/1000)))+
  geom_hline(yintercept=155/100,linetype="dotted",colour="red")+
  geom_hline(yintercept=1251.728,linetype="dashed")+
  geom_hline(yintercept=802.284,linetype="dashed")+
  geom_smooth(method="lm",se=F,fullrange=F,formula=y~x)+xlim(1995,2018)+
  geom_point(size=2)+theme_bw()+facet_wrap(~Strata)+
  ggtitle(label="a. krill biomass")+
  scale_colour_manual(values=c("red","blue"))+
  ylab("million tonnes")+xlab("") # figure X


plot2b<-ggplot(subset(kmfcr2,Strata=="Bransfield"|Strata=="West"),
               aes(Year,(CatchK*100)))+
  geom_hline(yintercept=5,linetype="dotted",colour="red")+geom_point()+
  
  geom_smooth(method="lm",se=F,fullrange=F,span=0.9)+xlim(1995,2018)+
  geom_point(size=2)+theme_bw()+facet_wrap(~Strata)+
  ggtitle(label="b. krill caught biomass")+
  scale_colour_manual(values=c("red","blue"))+
  ylab("% of catched biomass") # figure X


plot2c<-ggplot(subset(kmfcr2,Strata=="Bransfield"|Strata=="West"),
       aes((amlrbiomass/1000),(CatchK*100)))+
  geom_hline(yintercept=5,linetype="dashed")+geom_point()+
  geom_vline(xintercept=1251.728,linetype="dashed")+
  geom_vline(xintercept=802.284,linetype="dashed")+
  geom_smooth(method="gam",se=F,fullrange=F,span=0.9)+
  geom_point(size=2)+theme_bw()+facet_wrap(~Strata)+
  ggtitle(label="c. availability vs catch")+
  scale_colour_manual(values=c("red","blue"))+
  ylab("% of catched biomass")+ 
  xlab("Krill biomass (million tonnes)") # figure X

plot2a/plot2b/plot2c # figure 2

```

### Bayesian generalized linear mixed models. PS: results can vary slightly each run

```{r}
prior.m <- list(
  R=list(V=1, n=1, fix=1))

model1<-MCMCglmm(amlrbiomass~Year, data=subset(kmfcr2,Strata=="West"),random=NULL, verbose=FALSE,
                 nitt=10999, burnin=1000, thin=3,family = "gaussian",prior=prior.m)
summary(model1)

length(kmfcr2$amlrbiomass[kmfcr2$Strata=="West"])

((summary(model1)$solutions[2,4]))/10999 # F values
### effective sample size (eff. samp) is the number of the iteration when model converged
### an F measure was calculated dividing  eff. samp by number of iterations.


plot(model1$Sol)


length(kmfcr2$amlrbiomass[kmfcr2$Strata=="Bransfield"])

model2<-MCMCglmm(amlrbiomass~Year, data=subset(kmfcr2,Strata=="Bransfield"),random=NULL, verbose=FALSE,
                 nitt=10999, burnin=1000, thin=1,family = "gaussian",prior=prior.m)

summary(model2)



((summary(model2)$solutions[2,4])/10999)# F values

plot(model2$Sol)

model3<-MCMCglmm(CatchK~Year, data=subset(kmfcr2,Strata=="West"),random=NULL, verbose=FALSE,
                 nitt=10999, burnin=1000, thin=1,family = "exponential",prior=prior.m)

summary(model3)

(summary(model3)$solutions[2,4])/10999 # F values

plot(model3$Sol)

model4<-MCMCglmm(CatchK~Year, data=subset(kmfcr2,Strata=="Bransfield"),random=NULL, verbose=FALSE,
                 nitt=10999, burnin=1000,thin=1,family = "exponential",prior=prior.m)

summary(model4)

10999/(summary(model4)$solutions[2,4]) # F values

plot(model4$Sol)

model5<-MCMCglmm(CatchK~amlrbiomass, data=subset(kmfcr2,Strata=="Bransfield"|Strata=="West"),random=NULL, verbose=FALSE,
                 nitt=10999, burnin=1000, thin=1,family = "exponential")

summary(model5)

(summary(model5)$solutions[2,4])/10999 # F values

plot(model5$Sol)



```
###--------penguins-------------

#load data

### data from MAPPPD (penguinmap.org): Humphries et al. 2017 https://doi.org/10.1017/S0032247417000055

```{r}
peng<-read.csv("D:/BACKUP_Notebook_06_2022/KidsArentAlright/Data/CountQuery_V_4_0.csv") # query for area 48.1 

nests<-subset(peng,count_type=="nests",
              select=c(site_id,longitude_epsg_4326,latitude_epsg_4326,
                       common_name,day,month,season_starting,penguin_count,
                       count_type,vantage,accuracy))


chicks<-subset(peng,count_type=="chicks",
               select=c(site_id,longitude_epsg_4326,latitude_epsg_4326,
                        common_name,day,month,season_starting,penguin_count,
                        count_type,vantage,accuracy))
nestsM<-ddply(nests, c("common_name","site_id","season_starting"), summarise,
              nests=max(penguin_count),
              Lat=mean(latitude_epsg_4326),
              Long=mean(longitude_epsg_4326))
chicksM<-ddply(chicks, c("common_name","site_id","season_starting"), summarise,
               chicks=min(penguin_count),
               Lat=mean(latitude_epsg_4326),
               Long=mean(longitude_epsg_4326))

cpn<-merge(nestsM,chicksM)
cpn$chickspernest<-cpn$chicks/cpn$nests  ### measure of breeding success (chicks raised per nest)

cpn2<-data.frame(site=cpn$site_id,SP=cpn$common_name,
                 Year=cpn$season_starting,
                 Long=as.numeric(cpn$Long),
                 Lat=as.numeric(cpn$Lat),
                 chickspernest=cpn$chickspernest) # cpn stands for Chicks Per Nest

#### this is for summarizing number of populations and number of counts

cpnset<-cpn2 

counts<-ddply(cpnset, c("SP"), summarise,
              ncounts=length(chickspernest))

counts

pops<-ddply(cpn2, c("SP","site"), summarise,
            pops=length(chickspernest),
            Lat=mean(Lat),Long=mean(Long))


pops2<-ddply(pops, c("SP"), summarise,
             Npops=length(pops))

pops2


### -------------

cpn$chickspernest<-cpn$chicks/cpn$nests  ### measure of breeding success (chicks raised per nest)

summary(cpn$nests)

ggplot(cpn,aes(nests))+geom_histogram()+xlim(0,10000)
ggplot(cpn,aes(nests))+geom_histogram()+xlim(0,100) # check small populations
ggplot(cpn,aes(nests,chickspernest))+geom_point() # a single large population will be excluded too


### eliminate pops with success = 2, as they were two cases of one or two nests
### eliminate the single population with 80 000 nests

cpn<-subset(cpn,chickspernest<2 & nests<40000)
summary(cpn$chickspernest)


shapiro.test(cpn$chickspernest) # not normal 

### test if chicks per nest varied in time using sites and species as random factors

mc1<-MCMCglmm(chickspernest~season_starting, random=~site_id:common_name, 
              
              family="gaussian",
              data=cpn, nodes="ALL",  scale=TRUE, nitt=10999, thin=1, burnin=1000, pr=T,
              pl=FALSE, verbose=TRUE, DIC=TRUE, singular.ok=FALSE, saveX=TRUE,
              saveZ=TRUE, saveXL=TRUE, slice=FALSE, ginverse=NULL, trunc=FALSE,
              prior= list(
                R=list(V=1, n=1, fix=1, nu=0.002),
                G=list(G1=list(V=diag(1),
                               nu=0.002,
                               alpha.mu=rep(0, 1),
                               alpha.V=diag(1)*3)))
)

summary(mc1)

(summary(mc1)$solutions[2,4])/10999# F values

#plot(mc1$Sol) # diagnostic plots

sol<-data.frame(mc1$Sol) # random effects
solm<-melt(sol,id.vars=c("X.Intercept.","season_starting"))
solm$site_id<-substring(solm$variable,first=21,last=24)
solm$species<-substring(solm$variable,first=26,last=29)
sites<-ddply(cpn, c("site_id"), summarise,
             pops=length(chickspernest),
             Lat=mean(Lat),Long=mean(Long))

ranef<-ddply(solm, c("site_id","species"), summarise,
             int=mean(value),intsd=sd(value),
             intse=intsd/sqrt(length(value)-1))

rlat<-merge(ranef,sites,by="site_id")

model.lat<-MCMCglmm(int~Lat, data=rlat,random=NULL, verbose=FALSE,
                    nitt=10999, burnin=1000, thin=1,family = "gaussian")

summary(model.lat)


(summary(model.lat)$solutions[2,4])/10999 # F values

model.sp<-MCMCglmm(int~species, data=rlat,random=NULL, verbose=FALSE,
                    nitt=999, burnin=300, thin=1,family = "gaussian")
summary(model.sp)

(summary(model.sp)$solutions[2,4])/10999 # F for chinstrap vs adelie
(summary(model.sp)$solutions[3,4])/10999 # F for gentoo vs adelie

#plot(model.lat$Sol) # diagnostic plots


ggplot(cpn,aes(season_starting,chickspernest,colour=common_name,shape=common_name,linetype=common_name))+
  geom_smooth(method="lm",se=F)+
  geom_point(size=1,alpha=0.25)+
  xlab("Year")+ylab("Chick raised per nest")+
  theme_bw()+
  scale_colour_manual(values=c("blue","green3","red"))+
  scale_linetype_manual(values=c("dashed","dotted","dotdash"))+
  ggtitle(label="a. Temporal trends")+

ggplot()+
  geom_hline(yintercept=0,linetype="dashed")+
  geom_smooth(method="lm",data=rlat,aes(x=-1*Lat,y=int),colour="black")+
  geom_smooth(method="lm",se=F,data=rlat,aes(x=-1*Lat,y=int,colour=species,linetype=species))+
  geom_errorbar(data=rlat,aes(x=-1*Lat,y=int,ymin=int-intse,ymax=int+intse,colour=species,linetype=species),alpha=0.25)+
  geom_point(data=rlat,aes(-1*Lat,int,colour=species,shape=species),size=2,alpha=0.25)+theme_bw()+
  xlab("Latitude (°S)")+
  ylab("Random intercept")+
  scale_colour_manual(values=c("blue","green3","red"))+
  scale_linetype_manual(values=c("dashed","dotted","dotdash"))+
  ggtitle(label="b. Geographical trends")

```

### Penguins and fishery
### transfer strata information into penguin data

```{r}
pcoords<-coordinates(data.frame(cpnset[4:5])) 

p.spdf<-SpatialPointsDataFrame(pcoords,cpnset,proj4string = crs)

identicalCRS(p.spdf,strata)

p.strata<-over(p.spdf,strata)

p.spdf$Strata<-p.strata$Strata

p.spdf$StrataP<-p.spdf$Strata


####merge penguins and fishery
p.df<-data.frame(p.spdf)

peng<-data.frame(p.df[1:8]) ## maintain in the data only variables to analyze

summary(as.factor(peng$StrataP))

kmfcr<-merge(kr,fcr2,by=c("Strata","Year"),all=T) # re-merge krill and fishery data

kmfcr$catchK<-(kmfcr$CatchS/kmfcr$amlrbiomass)*100 # caught krill biomass (harvest rate)

tail(kmfcr) ### summarized by strata and year, catch from Dec to Mar

peng.f<-merge(peng,kmfcr,all.x=T)  ## merge and maintain all penguin data

pengf<-subset(peng.f,StrataP==Strata)#if penguin and fishing strata coincides, then it is right

pengf$CatchK<-pengf$CatchS/pengf$amlrbiomass

ggplot(pengf,aes(CatchK,chickspernest))+ #figure supp2

  geom_point()+
  theme_bw()+th+
  xlab("% of caught biomass")+ylab("chicks raised per nest")+
  geom_vline(xintercept=0.05,linetype="dashed") #5% catch

pengf$CatchH<-ifelse(pengf$CatchK>0.05,">5%","<5%")

pengf$KrillC<-ifelse(pengf$amlrbiomass<767000,"Krill<med","Krill>med")

pengf<-na.omit(pengf)
ddply(pengf, c("CatchH","KrillC"), summarise,
      ncounts=length(chickspernest))

summary(as.factor(pengf$CatchH))
summary(as.factor(pengf$KrillC))


shapiro.test(pengf$chickspernest) # not normal 

bartlett.test(pengf$chickspernest,pengf$CatchH) # homogeneity

bartlett.test(pengf$chickspernest,pengf$KrillC) # homogeneity

##-------permutation tests-------

pengf$kcch<-paste(pengf$KrillC,pengf$CatchH,sep=", Catch")

pengf$kcch<-factor(pengf$kcch,levels=c("Krill>med, Catch<5%","Krill<med, Catch<5%","Krill<med, Catch>5%"))

mc2<-MCMCglmm(chickspernest~kcch, random=~site, 
              
              family="gaussian",
              data=pengf, nodes="ALL",  scale=TRUE, nitt=10999, thin=10, burnin=1000, pr=T,
              pl=FALSE, verbose=TRUE, DIC=TRUE, singular.ok=FALSE, saveX=TRUE,
              saveZ=TRUE, saveXL=TRUE, slice=FALSE, ginverse=NULL, trunc=FALSE)




summary(mc2)

(summary(mc2)$solutions[2,4])/10999 # F values


#plot(mc2$Sol) # diagnostic plots
sol2<-data.frame(mc2$Sol) # random effects

solm2<-melt(sol2,id.vars=c("X.Intercept."))

solm2b<-solm2[2001:11000,] # subset only the sites

b = data.frame(X.Intercept.=gl(1,170),variable=c("ref"),value=as.numeric(gl(1, 170))*0) # add a zero reference for comparison
solm2b<-rbind(solm2b,b)

mcsite<-MCMCglmm(value~variable, data=solm2b,random=NULL, verbose=FALSE,
         nitt=999, burnin=300, thin=1,family = "gaussian")

summary(mcsite)
site.ranef<-data.frame(mcsite$Sol)

sr2<-melt(site.ranef,id.vars=c("X.Intercept."))
tail(sr2)

sr2m<-ddply(sr2, c("variable"), summarise,
             int=mean(value),intsd=sd(value),
             intse=intsd/sqrt(length(value)-1))

summary(as.factor(sr2m$variable))

ggplot(pengf,aes(kcch,chickspernest))+geom_boxplot()+
  geom_hline(yintercept=1,linetype="dashed",colour="grey50")+
  ylab("Chicks raised per nest")+
  theme_bw()+
  xlab("Categories")+
  ggtitle(label="a. Catched krill biomass")+
  coord_flip()+
  
ggplot(sr2m,aes(reorder(variable,-int),int))+
  geom_hline(yintercept=0,linetype="dashed",colour="grey50")+
  geom_errorbar(aes(ymin=int-intsd,ymax=int+intsd))+
  geom_point(size=2)+
  coord_flip()+theme_bw()+
  xlab("Breeding sites")+ylab("Random intercept")

 

### that's all folks!

```


