#Load necessary packages

## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## udunits system database read from C:/Users/Lucas Kruger/AppData/Local/R/win-library/4.2/udunits2/share/udunits/udunits2.xml
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## Loading required package: Matrix
## Loading required package: coda
## Loading required package: ape
## 
## Attaching package: 'ape'
## The following objects are masked from 'package:raster':
## 
##     rotate, zoom
## The following object is masked from 'package:dplyr':
## 
##     where

#theme for the plots:

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

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;”

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

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

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
##       Year         Strata              CatchS       
##  Min.   :1980   Length:146         Min.   :    0.7  
##  1st Qu.:1989   Class :character   1st Qu.:  142.4  
##  Median :1999   Mode  :character   Median : 1769.0  
##  Mean   :1999                      Mean   : 6604.6  
##  3rd Qu.:2009                      3rd Qu.: 7939.7  
##  Max.   :2017                      Max.   :48628.1
head(krill) ### amlrbiomass is in ton
##      Strata Year N wtgm2 wtvardens    cv  areakm2  minarea  maxarea strataarea
## 1 Joinville 2019 1 83.01    723.28 32.40 18151.00 18151.00 18151.00      18151
## 2 Joinville 2011 1 76.70   2027.06 58.70 18156.00 18156.00 18156.00      18156
## 3 Joinville 2010 1 17.10    102.83 59.30 18112.24 18112.24 18112.24      18112
## 4 Joinville 2009 1 30.60    149.82 40.00 18153.26 18153.26 18153.26      18153
## 5 Joinville 2008 1 68.80    191.24 20.10 18140.77 18140.77 18140.77      18141
## 6 Joinville 2005 2 11.72     57.55 53.99 18247.16 18172.22 18322.10      18322
##   stratabiomass cvbiomass areaamlr amlrbiomass cvamlrbiomm
## 1       1506746     32.40    18151     1506746       32.40
## 2       1392565     58.70    18151     1392182       58.70
## 3        309719     59.30    18151      310382       59.30
## 4        555490     40.00    18151      555421       40.00
## 5       1248085     20.10    18151     1248789       20.10
## 6        214762     64.72    18151      212756       64.72
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"]) 
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   39166  392888  615157  702795  892260 1975455      22
summary(kmfcr2$amlrbiomass[kmfcr2$Strata=="West"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   53934  752566 1482018 2072445 2837100 5466556      24
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"])
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##      33    2433    7806    8163   13133   20757      20
summary(kmfcr2$CatchS[kmfcr2$Strata=="Bransfield"])
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
##    19.24   241.11   385.02 10862.34 12654.15 48628.09       23
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
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'

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

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)
## 
##  Iterations = 1001:10997
##  Thinning interval  = 3
##  Sample size  = 3333 
## 
##  DIC: 5.664695e+13 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: amlrbiomass ~ Year 
## 
##             post.mean  l-95% CI  u-95% CI eff.samp  pMCMC    
## (Intercept) 100613194 100613071 100613302     3333 <3e-04 ***
## Year           -49248    -49248    -49248     3333 <3e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
length(kmfcr2$amlrbiomass[kmfcr2$Strata=="West"])
## [1] 78
((summary(model1)$solutions[2,4]))/10999 # F values
## [1] 0.3030275
### 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"])
## [1] 76
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)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: 4.138123e+12 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: amlrbiomass ~ Year 
## 
##             post.mean l-95% CI u-95% CI eff.samp  pMCMC    
## (Intercept)  21165506 21165373 21165640     9999 <1e-04 ***
## Year           -10214   -10214   -10214     9999 <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
((summary(model2)$solutions[2,4])/10999)# F values
## [1] 0.9090826
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)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: -182.3794 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: CatchK ~ Year 
## 
##              post.mean   l-95% CI   u-95% CI eff.samp pMCMC
## (Intercept)  -90.15603 -254.18554   86.62903     1448 0.299
## Year           0.04743   -0.04083    0.12917     1448 0.270
(summary(model3)$solutions[2,4])/10999 # F values
## [1] 0.1316058
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)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: -162.0615 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: CatchK ~ Year 
## 
##             post.mean l-95% CI u-95% CI eff.samp  pMCMC    
## (Intercept)  457.0733 255.8857 660.5390     1169 <1e-04 ***
## Year          -0.2252  -0.3262  -0.1248     1170 <1e-04 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
10999/(summary(model4)$solutions[2,4]) # F values
## [1] 9.400214
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)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: -342.2426 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units      2.68   0.6857    4.873    618.2
## 
##  Location effects: CatchK ~ amlrbiomass 
## 
##              post.mean   l-95% CI   u-95% CI eff.samp  pMCMC    
## (Intercept)  4.489e+00  3.654e+00  5.317e+00     1766 <1e-04 ***
## amlrbiomass  4.107e-07 -2.948e-08  8.623e-07     2650 0.0786 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(summary(model5)$solutions[2,4])/10999 # F values
## [1] 0.2409323
plot(model5$Sol)

###——–penguins————-

#load data

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

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
##                  SP ncounts
## 1    adelie penguin      72
## 2 chinstrap penguin      81
## 3    gentoo penguin     163
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
##                  SP Npops
## 1    adelie penguin    21
## 2 chinstrap penguin    18
## 3    gentoo penguin    41
### -------------

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

summary(cpn$nests)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1     295     660    1687    2119   77515
ggplot(cpn,aes(nests))+geom_histogram()+xlim(0,10000)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (`stat_bin()`).
## Warning: Removed 2 rows containing missing values (`geom_bar()`).

ggplot(cpn,aes(nests))+geom_histogram()+xlim(0,100) # check small populations
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 277 rows containing non-finite values (`stat_bin()`).
## Removed 2 rows containing missing values (`geom_bar()`).

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)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.9961  1.2000  1.1525  1.3632  1.9419
shapiro.test(cpn$chickspernest) # not normal 
## 
##  Shapiro-Wilk normality test
## 
## data:  cpn$chickspernest
## W = 0.96701, p-value = 1.472e-06
### 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)))
)
## 
##                        MCMC iteration = 0
## 
##                        MCMC iteration = 1000
## 
##                        MCMC iteration = 2000
## 
##                        MCMC iteration = 3000
## 
##                        MCMC iteration = 4000
## 
##                        MCMC iteration = 5000
## 
##                        MCMC iteration = 6000
## 
##                        MCMC iteration = 7000
## 
##                        MCMC iteration = 8000
## 
##                        MCMC iteration = 9000
## 
##                        MCMC iteration = 10000
summary(mc1)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: 615.0645 
## 
##  G-structure:  ~site_id:common_name
## 
##                     post.mean  l-95% CI u-95% CI eff.samp
## site_id:common_name  0.003286 7.379e-11  0.01339     1477
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units         1        1        1        0
## 
##  Location effects: chickspernest ~ season_starting 
## 
##                 post.mean  l-95% CI  u-95% CI eff.samp pMCMC
## (Intercept)       3.38279 -24.84284  31.07375     9999 0.818
## season_starting  -0.00111  -0.01486   0.01297     9999 0.884
(summary(mc1)$solutions[2,4])/10999# F values
## [1] 0.9090826
#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)
## 
##  Iterations = 1001:10999
##  Thinning interval  = 1
##  Sample size  = 9999 
## 
##  DIC: -682.2166 
## 
##  R-structure:  ~units
## 
##       post.mean  l-95% CI  u-95% CI eff.samp
## units 8.112e-06 5.634e-06 1.088e-05     9380
## 
##  Location effects: int ~ Lat 
## 
##              post.mean   l-95% CI   u-95% CI eff.samp  pMCMC  
## (Intercept) -3.773e-02 -7.383e-02 -2.628e-03     9999 0.0414 *
## Lat         -5.865e-04 -1.151e-03 -4.354e-05     9999 0.0408 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(summary(model.lat)$solutions[2,4])/10999 # F values
## [1] 0.9090826
model.sp<-MCMCglmm(int~species, data=rlat,random=NULL, verbose=FALSE,
                    nitt=999, burnin=300, thin=1,family = "gaussian")
summary(model.sp)
## 
##  Iterations = 301:999
##  Thinning interval  = 1
##  Sample size  = 699 
## 
##  DIC: -676.632 
## 
##  R-structure:  ~units
## 
##       post.mean  l-95% CI  u-95% CI eff.samp
## units 8.723e-06 6.145e-06 1.182e-05      699
## 
##  Location effects: int ~ species 
## 
##              post.mean   l-95% CI   u-95% CI eff.samp pMCMC
## (Intercept) -0.0003832 -0.0018715  0.0008477    699.0 0.567
## specieschin  0.0009858 -0.0010290  0.0029413    699.0 0.312
## speciesgent  0.0003777 -0.0011136  0.0020970    628.3 0.612
(summary(model.sp)$solutions[2,4])/10999 # F for chinstrap vs adelie
## [1] 0.06355123
(summary(model.sp)$solutions[3,4])/10999 # F for gentoo vs adelie
## [1] 0.05712695
#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")
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Penguins and fishery

transfer strata information into penguin data

pcoords<-coordinates(data.frame(cpnset[4:5])) 

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

identicalCRS(p.spdf,strata)
## [1] TRUE
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))
## Bransfield       West       NA's 
##         72         29        215
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
##     Strata Year amlrbiomass    CatchS catchK
## 163   <NA> 2011          NA    40.000     NA
## 164   <NA> 2013          NA 26111.099     NA
## 165   <NA> 2014          NA  9460.181     NA
## 166   <NA> 2015          NA 20227.648     NA
## 167   <NA> 2016          NA 14755.757     NA
## 168   <NA> 2017          NA  5583.101     NA
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
## Warning: Removed 35 rows containing missing values (`geom_point()`).

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))
##   CatchH    KrillC ncounts
## 1    <5% Krill<med      20
## 2    <5% Krill>med      32
## 3    >5% Krill<med      14
summary(as.factor(pengf$CatchH))
## <5% >5% 
##  52  14
summary(as.factor(pengf$KrillC))
## Krill<med Krill>med 
##        34        32
shapiro.test(pengf$chickspernest) # not normal 
## 
##  Shapiro-Wilk normality test
## 
## data:  pengf$chickspernest
## W = 0.9476, p-value = 0.007453
bartlett.test(pengf$chickspernest,pengf$CatchH) # homogeneity
## 
##  Bartlett test of homogeneity of variances
## 
## data:  pengf$chickspernest and pengf$CatchH
## Bartlett's K-squared = 0.40406, df = 1, p-value = 0.525
bartlett.test(pengf$chickspernest,pengf$KrillC) # homogeneity
## 
##  Bartlett test of homogeneity of variances
## 
## data:  pengf$chickspernest and pengf$KrillC
## Bartlett's K-squared = 2.3833, df = 1, p-value = 0.1226
##-------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)
## 
##                        MCMC iteration = 0
## 
##                        MCMC iteration = 1000
## 
##                        MCMC iteration = 2000
## 
##                        MCMC iteration = 3000
## 
##                        MCMC iteration = 4000
## 
##                        MCMC iteration = 5000
## 
##                        MCMC iteration = 6000
## 
##                        MCMC iteration = 7000
## 
##                        MCMC iteration = 8000
## 
##                        MCMC iteration = 9000
## 
##                        MCMC iteration = 10000
summary(mc2)
## 
##  Iterations = 1001:10991
##  Thinning interval  = 10
##  Sample size  = 1000 
## 
##  DIC: 17.74603 
## 
##  G-structure:  ~site
## 
##      post.mean  l-95% CI u-95% CI eff.samp
## site  0.006433 2.382e-17  0.03399    43.42
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units   0.07095  0.04568  0.09871    157.5
## 
##  Location effects: chickspernest ~ kcch 
## 
##                         post.mean l-95% CI u-95% CI eff.samp  pMCMC    
## (Intercept)               1.19940  1.08142  1.32422   1000.0 <0.001 ***
## kcchKrill<med, Catch<5%  -0.11533 -0.25926  0.04039   1000.0  0.136    
## kcchKrill<med, Catch>5%  -0.27217 -0.46760 -0.09621    244.7  0.008 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(summary(mc2)$solutions[2,4])/10999 # F values
## [1] 0.09091736
#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")
## Warning in MCMCglmm(value ~ variable, data = solm2b, random = NULL, verbose =
## FALSE, : some fixed effects are not estimable and have been removed. Use
## singular.ok=TRUE to sample these effects, but use an informative prior!
summary(mcsite)
## 
##  Iterations = 301:999
##  Thinning interval  = 1
##  Sample size  = 699 
## 
##  DIC: -23550.62 
## 
##  R-structure:  ~units
## 
##       post.mean l-95% CI u-95% CI eff.samp
## units  0.004487 0.004363 0.004626      699
## 
##  Location effects: value ~ variable 
## 
##                    post.mean   l-95% CI   u-95% CI eff.samp  pMCMC   
## (Intercept)       -0.0002107 -0.0094967  0.0090742      699 0.9728   
## variablesite.AITC  0.0103891  0.0001903  0.0215383      699 0.0715 . 
## variablesite.ARDL  0.0284930  0.0178457  0.0375618      699 <0.001 **
## variablesite.BART  0.0251757  0.0150031  0.0356492      699 <0.001 **
## variablesite.ENTR  0.0042627 -0.0059767  0.0146530      699 0.4721   
## variablesite.HALF -0.0160872 -0.0266030 -0.0064322      699 <0.001 **
## variablesite.HANN  0.0291904  0.0199302  0.0405445      699 <0.001 **
## variablesite.LLAN -0.0437257 -0.0538236 -0.0330851      699 <0.001 **
## variablesite.PTHO -0.0091844 -0.0188690  0.0013298      699 0.0973 . 
## variablesite.SHIR -0.0175887 -0.0271777 -0.0072101      699 <0.001 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
site.ranef<-data.frame(mcsite$Sol)

sr2<-melt(site.ranef,id.vars=c("X.Intercept."))
tail(sr2)
##       X.Intercept.          variable       value
## 6286  0.0002079258 variablesite.SHIR -0.01806066
## 6287 -0.0039782162 variablesite.SHIR -0.01246874
## 6288 -0.0003647849 variablesite.SHIR -0.01770216
## 6289 -0.0075576452 variablesite.SHIR -0.01058533
## 6290  0.0041808917 variablesite.SHIR -0.01936429
## 6291 -0.0061162296 variablesite.SHIR -0.01311118
sr2m<-ddply(sr2, c("variable"), summarise,
             int=mean(value),intsd=sd(value),
             intse=intsd/sqrt(length(value)-1))

summary(as.factor(sr2m$variable))
## variablesite.AITC variablesite.ARDL variablesite.BART variablesite.ENTR 
##                 1                 1                 1                 1 
## variablesite.HALF variablesite.HANN variablesite.LLAN variablesite.PTHO 
##                 1                 1                 1                 1 
## variablesite.SHIR 
##                 1
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!