Variability in krill biomass combined with fixed catch limits affected Pygoscelis penguin breeding success in the Western Antarctic Peninsula
Lucas Krüger, Francisco Santa Cruz, Mauricio Mardones, César Cárdenas
2023-03-17
#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.1head(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.72kr<-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 22summary(kmfcr2$amlrbiomass[kmfcr2$Strata=="West"])## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 53934 752566 1482018 2072445 2837100 5466556 24kmfcr2$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 20summary(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 23kmfcr2$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 ' ' 1length(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] 76model2<-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.9090826plot(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.1316058plot(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 ' ' 110999/(summary(model4)$solutions[2,4]) # F values## [1] 9.400214plot(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.2409323plot(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 163pops<-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 77515ggplot(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.9419shapiro.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 = 10000summary(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.9090826model.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] TRUEp.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 215kmfcr<-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 NApeng.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 14summary(as.factor(pengf$CatchH))## <5% >5%
## 52 14summary(as.factor(pengf$KrillC))## Krill<med Krill>med
## 34 32shapiro.test(pengf$chickspernest) # not normal ##
## Shapiro-Wilk normality test
##
## data: pengf$chickspernest
## W = 0.9476, p-value = 0.007453bartlett.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.525bartlett.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 = 10000summary(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 ' ' 1site.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.01311118sr2m<-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
## 1ggplot(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!