Appendix B
###########################################################################
################### 2 Habitat 4 Species Scenario ###############################
###########################################################################
rm(list=ls())
require("mgcv")
require("nFactors")
#Define Habitat
# 1 DeepReef
# 2 ShallowReef
# Select between control = 0 and test =1 scenario
scenario = 1
# set up
nspec <- 4
reps <- 500
years <- 20
year <- c(1:years)
# tweedie model parameters
p=1.3
disp=10
# determine the probability for a habitat switch during trip t
pSwitch <- 0.25
# Set up data generation
# Species order
# KOB, GLBK, HAKE, PANG
sp.names = c("KOB","GLBK","HAKE","PANG")
#random r
r <- runif(nspec)*0.1*(-1)^round((runif(nspec)*4),0)
#random biomass trends
B1 <- c(200*exp(rnorm(3,0,0.5)),50*exp(rnorm(1,0,0.5)))
### get biomass trajectories for 20 years
B <- mat.or.vec(years,nspec)
for (i in 1:nspec)
{
B[,i] <- B1[i] * exp((year-1) * r[i])
}
biomass_dat <- data.frame(year,B)
ylim = c(0,max(biomass_dat*1.05))
windows(width=7, height=5)
matplot(year, B,type="b",ylim=ylim, ylab="Biomass",pch=1:4,main="Biomass dynamics")
# Catchabilties for species i and habitat h
# KOB, GLBK, HAKE, PANG
q1 <- c(1.0,0.9,0.05,0.1)
q2 <- c(0.02,0.1,1.0,0.8)
# Effort error
error <- 0.2
#dataset structure
sim_year <- rep(c(1:20),each = reps)
KOB <- rep(B[,1], each = reps)
GLBK <- rep(B[,2], each = reps)
HAKE <- rep(B[,3], each = reps)
PANG <- rep(B[,4], each = reps)
n <- length(sim_year)
habitats <- c(1:2)
# Produce pi(h,y)
if(scenario==0)
{
vecP1 <- c(rep(0.3,20))*exp(rnorm(years,0,error))
vecP2 <- 1 - (vecP1)
}
if(scenario==1)
{
vecP1 <- c(0.70+0*(c(1:10)-1),0.3+0*(c(1:10)-1))*exp(rnorm(years,0,error))
vecP2 <- 1 - (vecP1)
}
#effort trends pi(h,y)
effort<- (rbind(vecP1,vecP2))
rand.effort <- effort[sample(1:2,replace=F),]
windows(height=4,width=8)
par(mfrow=c(1,1),mex=0.75, cex=1, cex.lab=1.1,mar=c(5,5,4,1))
plot(year,vecP1,type="b",pch=0, xlab="Year",ylab=expression(pi[h][","][y]),main="Effort dynamics",ylim=c(0,1))
lines(year,vecP2,type="b",pch=1,lty=2)
#probabilities for habitat to trip assignments
prob1 <- rep((rand.effort[1,]), each=reps)
# assigng up to 2 habitats to each fishing trip
habitatA <- ifelse(runif(n)<prob1,1,2)
habitatB <- ifelse(runif(n)<prob1,1,2)
###########################################################################
######################## Generate dataset ##################################### ###########################################################################
#Determine the fraction of effort (f(h,t))
propA <- 1/(1+exp(-(runif(n)-0.5)/0.5))
propB <- 1-propA
sim_dat0 <- data.frame(sim_year, habitatA, habitatB,propA,propB,KOB,GLBK,HAKE,PANG)
######## generate vectors with habitats 1 - 2
habA <- as.vector(sim_dat0$habitatA)
habB <- as.vector(sim_dat0$habitatB)
#######Generate CPUE records as a function of f(q,B,f)
Biomass <- sim_dat0[,6:(5+nspec)]
qA <- (if(habA[1] == 1) (rTweedie((q1*t(Biomass[1,])),p=p,disp)) else if(habA[1] == 2)
(rTweedie((q2*t(Biomass[1,])),p=p,disp)))
qB <- (if(runif(1)>pSwitch) (qA) else if(habB[1] == 1) (rTweedie((q1*t(Biomass[1,])),p=p,disp)) else
(rTweedie((q2*t(Biomass[1,])),p=p,disp)))
for(i in 2:n)
{
rqA <- (if(habA[i] == 1) (rTweedie((q1*t(Biomass[i,])),p=p,disp)) else
if(habA[i] == 2)(rTweedie((q2*t(Biomass[i,])),p=1.2,disp)))
qA <- rbind(qA,rqA)
qB <- rbind(qB,if(runif(1)>pSwitch) (rqA) else if(habB[i] == 1) (rTweedie((q1*t(Biomass[i,])),p=p,disp)) else if(habB[i] == 2)(rTweedie((q2*t(Biomass[i,])),p=p,disp)))
}
Bname <- data.frame(ifelse(Biomass>0,1,1))
cpueA <- propA*qA*Bname
cpueB <- propB*qB*Bname
cpue <- (cpueA+cpueB)
### create simulationed dataset
dat <- data.frame(year=sim_year,habitatA,habitatB,cpue)
ncols <- 3+nspec
dat1 <- subset(dat,apply(dat[,4:ncols], 1,sum)>0)
###########################################################################
########################### PCA ###########################################
###########################################################################
#extract catch composition
spec_dat <- dat1[,4:ncols]
spec_dat$sum <- apply(spec_dat, 1,sum)
mmult_dat <- subset(spec_dat,sum>0)
## prepare species composition for PCA
pca_dat <- (sqrt(sqrt(mmult_dat[,1:nspec]/mmult_dat$sum)))
# run PCA
pca<-prcomp(pca_dat, scale=T)
# Predict PCA Loadings
pr_pca <- predict(pca,pca_dat)
pcs <- data.frame(pr_pca[,1:nspec])
dataset <- cbind(dat1,pcs)
# get Eigenvalue
eig = summary(pca)$sdev^2
# run Optimal Coordinates test
OCtest = nScree(eig)
windows()
plotnScree(OCtest)
# retain number of PCs in combination with Eigenvalue > 1
nPC = min(OCtest$Components$nkaiser,OCtest$Components$noc)
###########################################################################
#################### Run standardization models for each species ###################
###########################################################################
# make CPUE histogram
br=50
windows(width=10, height=7)
par(mfrow=c(2,2),mex=0.65, cex=0.6, cex.lab=1,mar=c(5,5,1,1))
for (i in 1:nspec)
{
xlim <- c(0,0.5*max(cpue[,i]))
br <- round((max(cpue[,i])/log(max(cpue[,i])/20)),0)
hist(cpue[,i],breaks=br,col="gray",xlab="CPUE",main="",xlim=xlim)
legend('topright',paste(names(cpue[i])),bty="n")
}
#SPECIES LIST: KOB, GLBK, HAKE, PANG
# number of note points for splines
kn = 6
# run standardization models for each species
for(i in 1:nspec)
{
B <- biomass_dat[,i+1]
dataset$y <- dataset[,i+3]
#null model
m0 <- gam((y)~+1,data=dataset,family=Tweedie(p=p))
# fitting year-effect (Nominal)
m_base <- gam(y~as.factor(year),data=dataset,family=Tweedie(p=p,link="log"))
# select DPC variant
if(nPC==0)
m_dpc <- m_base
if(nPC==1)
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn),data = dataset,family=Tweedie(p=p,link="log"))
if(nPC==2) ##!!!!
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn)+s(PC2,bs="cr",k=kn),data = dataset,family=Tweedie(p=p,link="log"))
# calc. normalized trend for known biomass
trueB <- B/(mean(B))
## create prediction dataset
sel_pcs<- apply(dataset[,(4+nspec):(7+nspec)],2,FUN=mean)
tn <-1
tactic <- tn
for( j in 2:20)
{
sel_pcs <- rbind(apply(dataset[,(4+nspec):(7+nspec)],2,FUN=mean))
}
pred_dat <- data.frame(year=year,sel_pcs)
# predict nominal CPUE
bc <- (predict(m_base,pred_dat,type="response"))
bc_mu <- as.matrix(bc)
norm.bc_mu <- bc_mu/mean(bc_mu)
# predict DPC CPUE
pc <- (predict(m_dpc,pred_dat,type="response"))
pc_mu <- as.matrix(pc)
norm.pc_mu <- pc_mu/mean(pc_mu)
# fit log-linear trends
lm.bc <- lm(log(norm.bc_mu)~year)
lm.pc <- lm(log(norm.pc_mu)~year)
# get rhat
slope.bc <- exp(predict(lm.bc,data.frame(year)))
slope.pc <- exp(predict(lm.pc,data.frame(year)))
# Plot
if(i==1) windows(width=7, height=7.5)
if(i==1) par(mfrow=c(4,2),mex=0.65, cex=0.6, cex.lab=1,mar=c(5,5,1,1))
ylim = c(0,3)
plot(year,trueB,ylim=ylim,col=c(2),lty=c(1),lwd=c(2),type=("l"), ylab="rel. Abundance",xlab="Year",main="")
points(year,norm.bc_mu,cex=1.1)
lines(year,slope.bc,lty=5,lwd=2)
text(2,2.8,paste(sp.names[i]))
text(18,2.8,"Nominal")
plot(year,trueB,ylim=ylim,col=c(2),lty=c(1),lwd=c(2),type=("l"), ylab="rel. Abundance",xlab="Year",main="")
points(year,norm.pc_mu,cex=1.1)
lines(year,slope.pc,lty=5,lwd=2)
text(2,2.8,paste(sp.names[i]))
text(18,2.8,paste("PC",nPC,".R4",sep=""))
}
###### End of the 2 Habitat 4 Species simulation ###############################
###########################################################################
################### 4 Habitat 10 Species Scenario ###############################
###########################################################################
rm(list=ls())
require("mgcv")
require("nFactors")
# Define Habitat
# 1 DeepReef
# 2 ShallowReef
# 3 DeepSediment
# 4 ShallowSediment
# Select between control = 0 and test =1 scenario
scenario = 1
# set up
nspec <- 10
reps <- 500
sims <- 10
years <- 20
year <- c(1:years)
# tweedie model parameters
p=1.3
disp=10
# determine the probability for a habitat switch during trip t
pSwitch <- 0.25
# Set up data generation
# Species order
# KOB, GLBK, HAKE, PANG,CRPN,SNTR,ROMN,DRGD,RSTM,SHRK
sp.names = c("KOB","GLBK","HAKE","PANG","CRPN","SNTR","ROMN","DRGD","RSTM"
,"SHRK")
#random r
r <- runif(nspec)*0.1*(-1)^round((runif(nspec)*4),0)
#random biomass trends
B1 <- c(200*exp(rnorm(3,0,0.5)),50*exp(rnorm(1,0,0.5)),200*exp(rnorm(1,0,0.5)),50*
exp(rnorm(5,0,0.5)))
### get biomass trajectories for 20 years
B <- mat.or.vec(years,nspec)
for (i in 1:nspec)
{
B[,i] <- B1[i] * exp((year-1) * r[i])
}
biomass_dat <- data.frame(year,B)
ylim = c(0,max(biomass_dat*1.05))
windows(width=7, height=5)
matplot(year, B,type="b",ylim=ylim, ylab="Biomass",pch=1:4,main="Biomass dynamics")
# Catchabilties for species i and habitat h
# KOB, GLBK, HAKE, PANG,CRPN,SNTR,ROMN,DRGD,RSTM,SHRK
q1 <- c(1.0,0.9,0.05,0.1,0.01,0.01,0.01,0.01,0.01,0.5)
q2 <- c(0.02,0.1,1.0,0.8,0.01,0.01,0.01,0.01,0.01,0.05)
q3 <- c(0.01,0.01,0.01,0.01,1.0,0.05,0.2,0.2,0.3,0.05)
q4 <- c(0.01,0.01,0.01,0.01,0.1,1.0,0.8,0.7,0.9,0.3)
# Effort error
error <- 0.2
#dataset structure
sim_year <- rep(c(1:20),each = reps)
KOB <- rep(B[,1], each = reps)
GLBK <- rep(B[,2], each = reps)
HAKE <- rep(B[,3], each = reps)
PANG <- rep(B[,4], each = reps)
CRPN <- rep(B[,5], each = reps)
SNTR <- rep(B[,6], each = reps)
ROMN <- rep(B[,7], each = reps)
DGRD <- rep(B[,8], each = reps)
RSTM <- rep(B[,9], each = reps)
SHRK <- rep(B[,10], each = reps)
n <- length(sim_year)
half <- c(1:10)
habitats <- c(1:4)
if(scenario==0)
{
vecP1 <- c(rep(0.2,20))*exp(rnorm(years,0,error))
vecP2 <- c(rep(0.25,20))*exp(rnorm(years,0,error))
vecP3 <- c(rep(0.3,20))*exp(rnorm(years,0,error))
vecP4 <- 1 - (vecP1+vecP2+vecP3)
}
if(scenario==1)
{
vecP1 <- c(0.45+0*(c(1:10)-1),0.11+0*(c(1:10)-1))*exp(rnorm(years,0,error))
vecP2 <- c(0.1+0*(c(1:10)-1),0.45+0*(c(1:10)-1))*exp(rnorm(years,0,error))
vecP3 <- (0.3-0.008*(year-1))*exp(rnorm(years,0,error))
vecP4 <- 1 - (vecP1+vecP2+vecP3)
}
#effort trends pi(h,y)
effort<- (rbind(vecP1,vecP2,vecP3,vecP4))
rand.effort <- effort[sample(1:4,replace=F),]
windows(height=4,width=8)
par(mfrow=c(1,1),mex=0.75, cex=1, cex.lab=1.1,mar=c(5,5,1.5,1))
plot(year,vecP1,type="b",pch=0, xlab="Year",ylab=expression(pi[h][","][y]),main="Effort dynamics",ylim=c(0,1))
lines(year,vecP2,type="b",pch=1,lty=2)
lines(year,vecP3,type="b",pch=1,lty=3)
lines(year,vecP4,type="b",pch=1,lty=4)
#probabilities for habitat to trip assignments
prob1 <- rep((rand.effort[1,]), each=reps)
prob2 <- rep(rand.effort[2,]/(1-rand.effort[1,]), each=reps)
prob3 <- rep(rand.effort[3,]/(1-(rand.effort[1,]+rand.effort[2,])), each=reps)
# assigng up to 2 habitats to each fishing trip
habitatA <- ifelse(runif(n)<prob1,1,ifelse(runif(n)<prob2,2,ifelse(runif(n)<prob3,3,4)))
habitatB <- ifelse(runif(n)<prob1,1,ifelse(runif(n)<prob2,2,ifelse(runif(n)<prob3,3,4)))
###########################################################################
######################## Generate dataset ##################################### ###########################################################################
#Determine the fraction of effort (f(h,t))
propA <- 1/(1+exp(-(runif(n)-0.5)/0.5))
propB <- 1-propA
sim_dat0 <- data.frame(sim_year, habitatA, habitatB,propA,propB,KOB,GLBK,HAKE,PANG,CRPN,SNTR,ROMN,DGRD,RSTM,SHRK)
######## generate vectors with habitats 1 - 2
habA <- as.vector(sim_dat0$habitatA)
habB <- as.vector(sim_dat0$habitatB)
#######Generate CPUE records as a function of f(q,B,f)
Biomass <- sim_dat0[,6:(5+nspec)]
qA <- (if(habA[1] == 1) (rTweedie((q1*t(Biomass[1,])),p=p,disp)) else
if(habA[1] == 2)(rTweedie((q2*t(Biomass[1,])),p=p,disp)) else
if(habA[1] == 3)(rTweedie((q3*t(Biomass[1,])),p=p,disp)) else (rTweedie((q4*t(Biomass[1,])),p=p,disp)))
qB <- (if(runif(1)>pSwitch) (qA) else if(habB[1] == 1) (rTweedie((q1*t(Biomass[1,])),p=p,disp)) else
if(habB[1] == 2)(rTweedie((q2*t(Biomass[1,])),p=p,disp)) else
if(habB[1] == 3)(rTweedie((q3*t(Biomass[1,])),p=p,disp)) else (rTweedie((q4*t(Biomass[1,])),p=p,disp)))
for(i in 2:n)
{
rqA <- (if(habA[i] == 1) (rTweedie((q1*t(Biomass[i,])),p=p,disp)) else if(habA[i] == 2)(rTweedie((q2*t(Biomass[i,])),p=1.2,disp)) else if(habA[i] == 3)(rTweedie((q3*t(Biomass[i,])),p=1.2,disp)) else (rTweedie((q4*t(Biomass[i,])),p=p,disp)))
qA <- rbind(qA,rqA)
qB <- rbind(qB,if(runif(1)>pSwitch) (rqA) else
if(habB[i] == 1) (rTweedie((q1*t(Biomass[i,])),p=p,disp)) else
if(habB[i] == 2)(rTweedie((q2*t(Biomass[i,])),p=p,disp)) else
if(habB[i] == 3)(rTweedie((q3*t(Biomass[i,])),p=p,disp)) else (rTweedie((q4*t(Biomass[i,])),p=p,disp)))
}
Bname <- data.frame(ifelse(Biomass>0,1,1))
cpueA <- propA*qA*Bname
cpueB <- propB*qB*Bname
cpue <- (cpueA+cpueB)
### create simulationed dataset
dat <- data.frame(year=sim_year,habitatA,habitatB,cpue)
ncols <- 3+nspec
dat1 <- subset(dat,apply(dat[,4:ncols], 1,sum)>0)
###########################################################################
########################### PCA ###########################################
###########################################################################
#extract catch composition
spec_dat <- dat1[,4:ncols]
spec_dat$sum <- apply(spec_dat, 1,sum)
mmult_dat <- subset(spec_dat,sum>0)
## prepare species composition for PCA
pca_dat <- (sqrt(sqrt(mmult_dat[,1:nspec]/mmult_dat$sum)))
#run PCA
pca<-prcomp(pca_dat, scale=T)
# Predict PCA Loadings
pr_pca <- predict(pca,pca_dat)
pcs <- data.frame(pr_pca[,1:nspec])
dataset <- cbind(dat1,pcs)
# get Eigenvalue
eig = summary(pca)$sdev^2
# run Optimal Coordinates test
OCtest = nScree(eig)
windows()
plotnScree(OCtest)
# retain number of PCs in combination with Eigenvalue > 1
nPC = min(OCtest$Components$nkaiser,OCtest$Components$noc)
###########################################################################
#################### Run standardization models for each species ###################
###########################################################################
# make CPUE histogram
br=50
windows(width=10, height=12)
par(mfrow=c(5,2),mex=0.65, cex=0.6, cex.lab=1,mar=c(5,5,1,1))
for (i in 1:nspec)
{
xlim <- c(0,0.5*max(cpue[,i]))
br <- round((max(cpue[,i])/log(max(cpue[,i])/20)),0)
hist(cpue[,i],breaks=br,col="gray",xlab="CPUE",main="",xlim=xlim)
legend('topright',paste(names(cpue[i])),bty="n")
}
#SPECIES LIST: KOB, GLBK, HAKE, PANG, CRPN, SNTR, ROMN, DRGD, RSTM, SHRK
# number of note points for splines
kn = 6
# run standardization models for each species
for(i in 1:nspec)
{
B <- biomass_dat[,i+1]
dataset$y <- dataset[,i+3]
# fitting year-effect (Nominal)
m_base <- gam(y~as.factor(year),data=dataset,family=Tweedie(p=p))
# select DPC variant
if(nPC==1)
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn),data = dataset,family=Tweedie(p=p))
if(nPC==2)
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn)+s(PC2,bs="cr",k=kn),data = dataset,family=Tweedie(p=p))
if(nPC==3)
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn)+s(PC2,bs="cr",k=kn)+s(PC3,bs="cr",k=kn),data = dataset,family=Tweedie(p=p))
if(nPC>3)
m_dpc <- gam(y~as.factor(year)+s(PC1,bs="cr",k=kn)+s(PC2,bs="cr",k=kn)+s(PC3,bs="cr",k=kn)+s(PC3,bs="cr",k=kn),data = dataset,family=Tweedie(p=p))
# calc. normalized trend for known biomass
trueB <- B/(mean(B))
## create prediction dataset
sel_pcs<- apply(dataset[,(4+nspec):(7+nspec)],2,FUN=mean)
tn <-1
tactic <- tn
for( j in 2:20)
{
sel_pcs <- rbind(apply(dataset[,(4+nspec):(7+nspec)],2,FUN=mean))
}
pred_dat <- data.frame(year=year,sel_pcs)
# predict nominal CPUE
bc <- (predict(m_base,pred_dat,type="response"))
bc_mu <- as.matrix(bc)
norm.bc_mu <- bc_mu/mean(bc_mu)
# predict DPC CPUE
pc <- (predict(m_dpc,pred_dat,type="response"))
pc_mu <- as.matrix(pc)
norm.pc_mu <- pc_mu/mean(pc_mu)
lm.bc <- lm(log(norm.bc_mu)~year)
lm.pc <- lm(log(norm.pc_mu)~year)
# fit log-linear trends
slope.bc <- exp(predict(lm.bc,data.frame(year)))
slope.pc <- exp(predict(lm.pc,data.frame(year)))
if(i==1) windows(width=7., height=9)
if(i==1) par(mfrow=c(5,2),mex=0.65, cex=0.6, cex.lab=1,mar=c(5,5,1,1))
if(i==6) windows(width=7., height=9)
if(i==6) par(mfrow=c(5,2),mex=0.65, cex=0.6, cex.lab=1,mar=c(5,5,1,1))
ylim = c(0,3)
plot(year,trueB,ylim=ylim,col=c(2),lty=c(1),lwd=c(2),type=("l"), ylab="rel. Abundance",xlab="Year",main="")
points(year,norm.bc_mu,cex=1.1)
lines(year,slope.bc,lty=5,lwd=2)
text(2,2.8,paste(sp.names[i]))
text(18,2.8,"Nominal")
plot(year,trueB,ylim=ylim,col=c(2),lty=c(1),lwd=c(2),type=("l"), ylab="rel. Abundance",xlab="Year",main="")
points(year,norm.pc_mu,cex=1.1)
lines(year,slope.pc,lty=5,lwd=2)
text(2,2.8,paste(sp.names[i]))
text(18,2.8,paste("PC",nPC,".R4",sep=""))
}
###### End of the 4 Habitat 10 Species simulation ###############################