| chapter7 {MQMF} | R Documentation |
chapter7 The 67 R-code chunks from Surplus Production Models
Description
chapter7 is not an active function but rather acts as a repository for the various example code chunks found in chapter7. There are 67 r-code chunks in chapter7 You should, of course, feel free to use and modify any of these example chunks in your own work.
Usage
chapter7(verbose = TRUE)
Arguments
verbose |
Should instructions to written to the console, default=TRUE |
Examples
## Not run:
# All the example code from # Surplus Production Models
# Surplus Production Models
## Introduction
### Data Needs
### The Need for Contrast
### When are Catch-Rates Informative
# R-chunk 1 Page 256
#Yellowfin-tuna data from Schaefer 12957
# R-chunk 2 Page 256 Table 7.1 code not in the book
data(schaef)
kable(halftable(schaef,subdiv=2),digits=c(0,0,0,4))
# R-chunk 3 Page 256
#schaef fishery data and regress cpue and catch Fig 7.1
oldp <- parset(plots=c(3,1),margin=c(0.35,0.4,0.05,0.05))
plot1(schaef[,"year"],schaef[,"catch"],ylab="Catch",xlab="Year",
defpar=FALSE,lwd=2)
plot1(schaef[,"year"],schaef[,"cpue"],ylab="CPUE",xlab="Year",
defpar=FALSE,lwd=2)
plot1(schaef[,"catch"],schaef[,"cpue"],type="p",ylab="CPUE",
xlab="Catch",defpar=FALSE,pch=16,cex=1.0)
model <- lm(schaef[,"cpue"] ~ schaef[,"catch"])
abline(model,lwd=2,col=2) # summary(model)
par(oldp) # return par to old settings; this line not in book
# R-chunk 4 Page 257
#cross correlation between cpue and catch in schaef Fig 7.2
oldp <- parset(cex=0.85) #sets par values for a tidy base graphic
ccf(x=schaef[,"catch"],y=schaef[,"cpue"],type="correlation",
ylab="Correlation",plot=TRUE)
par(oldp) # return par to old settings; this line not in book
# R-chunk 5 Page 257
#now plot schaef data with timelag of 2 years on cpue Fig 7.3
oldp <- parset(plots=c(3,1),margin=c(0.35,0.4,0.05,0.05))
plot1(schaef[1:20,"year"],schaef[1:20,"catch"],ylab="Catch",
xlab="Year",defpar=FALSE,lwd=2)
plot1(schaef[3:22,"year"],schaef[3:22,"cpue"],ylab="CPUE",
xlab="Year",defpar=FALSE,lwd=2)
plot1(schaef[1:20,"catch"],schaef[3:22,"cpue"],type="p",
ylab="CPUE",xlab="Catch",defpar=FALSE,cex=1.0,pch=16)
model2 <- lm(schaef[3:22,"cpue"] ~ schaef[1:20,"catch"])
abline(model2,lwd=2,col=2)
par(oldp) # return par to old settings; this line not in book
# R-chunk 6 Page 259
#write out a summary of he regression model2
summary(model2)
## Some Equations
### Production Functions
# R-chunk 7 Page 262
#plot productivity and density-dependence functions Fig7.4
prodfun <- function(r,Bt,K,p) return((r*Bt/p)*(1-(Bt/K)^p))
densdep <- function(Bt,K,p) return((1/p)*(1-(Bt/K)^p))
r <- 0.75; K <- 1000.0; Bt <- 1:1000
sp <- prodfun(r,Bt,K,1.0) # Schaefer equivalent
sp0 <- prodfun(r,Bt,K,p=1e-08) # Fox equivalent
sp3 <- prodfun(r,Bt,K,3) #left skewed production, marine mammal?
oldp <- parset(plots=c(2,1),margin=c(0.35,0.4,0.1,0.05))
plot1(Bt,sp,type="l",lwd=2,xlab="Stock Size",
ylab="Surplus Production",maxy=200,defpar=FALSE)
lines(Bt,sp0 * (max(sp)/max(sp0)),lwd=2,col=2,lty=2) # rescale
lines(Bt,sp3*(max(sp)/max(sp3)),lwd=3,col=3,lty=3) # production
legend(275,100,cex=1.1,lty=1:3,c("p = 1.0 Schaefer","p = 1e-08 Fox",
"p = 3 LeftSkewed"),col=c(1,2,3),lwd=3,bty="n")
plot1(Bt,densdep(Bt,K,p=1),xlab="Stock Size",defpar=FALSE,
ylab="Density-Dependence",maxy=2.5,lwd=2)
lines(Bt,densdep(Bt,K,1e-08),lwd=2,col=2,lty=2)
lines(Bt,densdep(Bt,K,3),lwd=3,col=3,lty=3)
par(oldp) # return par to old settings; this line not in book
### The Schaefer Model
### Sum of Squared Residuals
### Estimating Management Statistics
# R-chunk 8 Page 266
#compare Schaefer and Fox MSY estimates for same parameters
param <- c(r=1.1,K=1000.0,Binit=800.0,sigma=0.075)
cat("MSY Schaefer = ",getMSY(param,p=1.0),"\n") # p=1 is default
cat("MSY Fox = ",getMSY(param,p=1e-08),"\n")
### The Trouble with Equilibria
## Model Fitting
### A Possible Workflow for Stock Assessment
# R-chunk 9 Page 269
#Initial model 'fit' to the initial parameter guess Fig 7.5
data(schaef); schaef <- as.matrix(schaef)
param <- log(c(r=0.1,K=2250000,Binit=2250000,sigma=0.5))
negatL <- negLL(param,simpspm,schaef,logobs=log(schaef[,"cpue"]))
ans <- plotspmmod(inp=param,indat=schaef,schaefer=TRUE,
addrmse=TRUE,plotprod=FALSE)
# R-chunk 10 Pages 270 - 271
#Fit the model first using optim then nlm in sequence
param <- log(c(0.1,2250000,2250000,0.5))
pnams <- c("r","K","Binit","sigma")
best <- optim(par=param,fn=negLL,funk=simpspm,indat=schaef,
logobs=log(schaef[,"cpue"]),method="BFGS")
outfit(best,digits=4,title="Optim",parnames = pnams)
cat("\n")
best2 <- nlm(negLL,best$par,funk=simpspm,indat=schaef,
logobs=log(schaef[,"cpue"]))
outfit(best2,digits=4,title="nlm",parnames = pnams)
# R-chunk 11 Page 271
#optimum fit. Defaults used in plotprod and schaefer Fig 7.6
ans <- plotspmmod(inp=best2$estimate,indat=schaef,addrmse=TRUE,
plotprod=TRUE)
# R-chunk 12 Page 272
#the high-level structure of ans; try str(ans$Dynamics)
str(ans, width=65, strict.width="cut",max.level=1)
# R-chunk 13 Page 273
#compare the parameteric MSY with the numerical MSY
round(ans$Dynamics$sumout,3)
cat("\n Productivity Statistics \n")
summspm(ans) # the q parameter needs more significantr digits
### Is the Analysis Robust?
# R-chunk 14 Page 274
#conduct a robustness test on the Schaefer model fit
data(schaef); schaef <- as.matrix(schaef); reps <- 12
param <- log(c(r=0.15,K=2250000,Binit=2250000,sigma=0.5))
ansS <- fitSPM(pars=param,fish=schaef,schaefer=TRUE, #use
maxiter=1000,funk=simpspm,funkone=FALSE) #fitSPM
#getseed() #generates random seed for repeatable results
set.seed(777852) #sets random number generator with a known seed
robout <- robustSPM(inpar=ansS$estimate,fish=schaef,N=reps,
scaler=40,verbose=FALSE,schaefer=TRUE,
funk=simpspm,funkone=FALSE)
#use str(robout) to see the components included in the output
# R-chunk 15 Page 275 Table 7.2 code not in the book
#outcome of robustness tests
kable(robout$results[,1:5],digits=c(3,4,3,4,3))
kable(robout$results[,6:11],digits=c(3,4,3,4,5,0))
# R-chunk 16 Pages 275 - 276
#Repeat robustness test on fit to schaef data 100 times
set.seed(777854)
robout2 <- robustSPM(inpar=ansS$estimate,fish=schaef,N=100,
scaler=25,verbose=FALSE,schaefer=TRUE,
funk=simpspm,funkone=TRUE,steptol=1e-06)
lastbits <- tail(robout2$results[,6:11],10)
# R-chunk 17 Page 276 Table 7.3 code not in the book
#last 10 rows of robustness test showing deviations
kable(lastbits,digits=c(5,1,1,4,5,0))
# R-chunk 18 Page 276
# replicates from the robustness test Fig 7.7
result <- robout2$results
oldp <- parset(plots=c(2,2),margin=c(0.35,0.45,0.05,0.05))
hist(result[,"r"],breaks=15,col=2,main="",xlab="r")
hist(result[,"K"],breaks=15,col=2,main="",xlab="K")
hist(result[,"Binit"],breaks=15,col=2,main="",xlab="Binit")
hist(result[,"MSY"],breaks=15,col=2,main="",xlab="MSY")
par(oldp) # return par to old settings; this line not in book
# R-chunk 19 Page 277
#robustSPM parameters against each other Fig 7.8
pairs(result[,c("r","K","Binit","MSY")],upper.panel=NULL,pch=1)
### Using Different Data?
# R-chunk 20 Page 278
#Now use the dataspm data-set, which is noisier
set.seed(777854) #other random seeds give different results
data(dataspm); fish <- dataspm #to generalize the code
param <- log(c(r=0.24,K=5174,Binit=2846,sigma=0.164))
ans <- fitSPM(pars=param,fish=fish,schaefer=TRUE,maxiter=1000,
funkone=TRUE)
out <- robustSPM(ans$estimate,fish,N=100,scaler=15, #making
verbose=FALSE,funkone=TRUE) #scaler=10 gives
result <- tail(out$results[,6:11],10) #16 sub-optimal results
# R-chunk 21 Page 279 Table 7.4 code not in the book
#last 10 trials of robustness on dataspm fit
kable(result,digits=c(4,2,2,4,4,3))
## Uncertainty
### Likelihood Profiles
# R-chunk 22 Page 280
# Fig 7.9 Fit of optimum to the abdat data-set
data(abdat); fish <- as.matrix(abdat)
colnames(fish) <- tolower(colnames(fish)) # just in case
pars <- log(c(r=0.4,K=9400,Binit=3400,sigma=0.05))
ans <- fitSPM(pars,fish,schaefer=TRUE) #Schaefer
answer <- plotspmmod(ans$estimate,abdat,schaefer=TRUE,addrmse=TRUE)
# R-chunk 23 Pages 280 - 282
# likelihood profiles for r and K for fit to abdat Fig 7.10
#doprofile input terms are vector of values, fixed parameter
#location, starting parameters, and free parameter locations.
#all other input are assumed to be in the calling environment
doprofile <- function(val,loc,startest,indat,notfix=c(2:4)) {
pname <- c("r","K","Binit","sigma","-veLL")
numv <- length(val)
outpar <- matrix(NA,nrow=numv,ncol=5,dimnames=list(val,pname))
for (i in 1:numv) { #
param <- log(startest) # reset the parameters
param[loc] <- log(val[i]) #insert new fixed value
parinit <- param # copy revised parameter vector
bestmod <- nlm(f=negLLP,p=param,funk=simpspm,initpar=parinit,
indat=indat,logobs=log(indat[,"cpue"]),notfixed=notfix)
outpar[i,] <- c(exp(bestmod$estimate),bestmod$minimum)
}
return(outpar)
}
rval <- seq(0.32,0.46,0.001)
outr <- doprofile(rval,loc=1,startest=c(rval[1],11500,5000,0.25),
indat=fish,notfix=c(2:4))
Kval <- seq(7200,11500,200)
outk <- doprofile(Kval,loc=2,c(0.4,7200,6500,0.3),indat=fish,notfix=c(1,3,4))
oldp <- parset(plots=c(2,1),cex=0.85,outmargin=c(0.5,0.5,0,0))
plotprofile(outr,var="r",defpar=FALSE,lwd=2) #MQMF function
plotprofile(outk,var="K",defpar=FALSE,lwd=2)
par(oldp) # return par to old settings; this line not in book
### Bootstrap Confidence Intervals
# R-chunk 24 Page 283
#find optimum Schaefer model fit to dataspm data-set Fig 7.11
data(dataspm)
fish <- as.matrix(dataspm)
colnames(fish) <- tolower(colnames(fish))
pars <- log(c(r=0.25,K=5500,Binit=3000,sigma=0.25))
ans <- fitSPM(pars,fish,schaefer=TRUE,maxiter=1000) #Schaefer
answer <- plotspmmod(ans$estimate,fish,schaefer=TRUE,addrmse=TRUE)
# R-chunk 25 Page 284
#bootstrap the log-normal residuals from optimum model fit
set.seed(210368)
reps <- 1000 # can take 10 sec on a large Desktop. Be patient
#startime <- Sys.time() # schaefer=TRUE is the default
boots <- spmboot(ans$estimate,fishery=fish,iter=reps)
#print(Sys.time() - startime) # how long did it take?
str(boots,max.level=1)
# R-chunk 26 Page 285
#Summarize bootstrapped parameter estimates as quantiles seen in Table 7.5
bootpar <- boots$bootpar
rows <- colnames(bootpar)
columns <- c(c(0.025,0.05,0.5,0.95,0.975),"Mean")
bootCI <- matrix(NA,nrow=length(rows),ncol=length(columns),
dimnames=list(rows,columns))
for (i in 1:length(rows)) {
tmp <- bootpar[,i]
qtil <- quantile(tmp,probs=c(0.025,0.05,0.5,0.95,0.975),na.rm=TRUE)
bootCI[i,] <- c(qtil,mean(tmp,na.rm=TRUE))
}
# R-chunk 27 page 285 # not visible in the book but this generates Table 7.5
kable(bootCI,digits=c(4,4,4,4,4,4))
# R-chunk 28 Page 286
#boostrap CI. Note use of uphist to expand scale Fig 7.12
colf <- c(1,1,1,4); lwdf <- c(1,3,1,3); ltyf <- c(1,1,1,2)
colsf <- c(2,3,4,6)
oldp <- parset(plots=c(3,2))
hist(bootpar[,"r"],breaks=25,main="",xlab="r")
abline(v=c(bootCI["r",colsf]),col=colf,lwd=lwdf,lty=ltyf)
uphist(bootpar[,"K"],maxval=14000,breaks=25,main="",xlab="K")
abline(v=c(bootCI["K",colsf]),col=colf,lwd=lwdf,lty=ltyf)
hist(bootpar[,"Binit"],breaks=25,main="",xlab="Binit")
abline(v=c(bootCI["Binit",colsf]),col=colf,lwd=lwdf,lty=ltyf)
uphist(bootpar[,"MSY"],breaks=25,main="",xlab="MSY",maxval=450)
abline(v=c(bootCI["MSY",colsf]),col=colf,lwd=lwdf,lty=ltyf)
hist(bootpar[,"Depl"],breaks=25,main="",xlab="Final Depletion")
abline(v=c(bootCI["Depl",colsf]),col=colf,lwd=lwdf,lty=ltyf)
hist(bootpar[,"Harv"],breaks=25,main="",xlab="End Harvest Rate")
abline(v=c(bootCI["Harv",colsf]),col=colf,lwd=lwdf,lty=ltyf)
par(oldp) # return par to old settings; this line not in book
# R-chunk 29 Page 286
#Fig7.13 1000 bootstrap trajectories for dataspm model fit
dynam <- boots$dynam
years <- fish[,"year"]
nyrs <- length(years)
oldp <- parset()
ymax <- getmax(c(dynam[,,"predCE"],fish[,"cpue"]))
plot(fish[,"year"],fish[,"cpue"],type="n",ylim=c(0,ymax),
xlab="Year",ylab="CPUE",yaxs="i",panel.first = grid())
for (i in 1:reps) lines(years,dynam[i,,"predCE"],lwd=1,col=8)
lines(years,answer$Dynamics$outmat[1:nyrs,"predCE"],lwd=2,col=0)
points(years,fish[,"cpue"],cex=1.2,pch=16,col=1)
percs <- apply(dynam[,,"predCE"],2,quants)
arrows(x0=years,y0=percs["5%",],y1=percs["95%",],length=0.03,
angle=90,code=3,col=0)
par(oldp) # return par to old settings; this line not in book
# R-chunk 30 Page 288
#Fit the Fox model to dataspm; note different parameters
pars <- log(c(r=0.15,K=6500,Binit=3000,sigma=0.20))
ansF <- fitSPM(pars,fish,schaefer=FALSE,maxiter=1000) #Fox version
bootsF <- spmboot(ansF$estimate,fishery=fish,iter=reps,schaefer=FALSE)
dynamF <- bootsF$dynam
# R-chunk 31 Pages 288 - 289
# bootstrap trajectories from both model fits Fig 7.14
oldp <- parset()
ymax <- getmax(c(dynam[,,"predCE"],fish[,"cpue"]))
plot(fish[,"year"],fish[,"cpue"],type="n",ylim=c(0,ymax),
xlab="Year",ylab="CPUE",yaxs="i",panel.first = grid())
for (i in 1:reps) lines(years,dynamF[i,,"predCE"],lwd=1,col=1,lty=1)
for (i in 1:reps) lines(years,dynam[i,,"predCE"],lwd=1,col=8)
lines(years,answer$Dynamics$outmat[1:nyrs,"predCE"],lwd=2,col=0)
points(years,fish[,"cpue"],cex=1.1,pch=16,col=1)
percs <- apply(dynam[,,"predCE"],2,quants)
arrows(x0=years,y0=percs["5%",],y1=percs["95%",],length=0.03,
angle=90,code=3,col=0)
legend(1985,0.35,c("Schaefer","Fox"),col=c(8,1),bty="n",lwd=3)
par(oldp) # return par to old settings; this line not in book
### Parameter Correlations
# R-chunk 32 Page 290
# plot variables against each other, use MQMF panel.cor Fig 7.15
pairs(boots$bootpar[,c(1:4,6,7)],lower.panel=panel.smooth,
upper.panel=panel.cor,gap=0,lwd=2,cex=0.5)
### Asymptotic Errors
# R-chunk 33 Page 290
#Start the SPM analysis using asymptotic errors.
data(dataspm) # Note the use of hess=TRUE in call to fitSPM
fish <- as.matrix(dataspm) # using as.matrix for more speed
colnames(fish) <- tolower(colnames(fish)) # just in case
pars <- log(c(r=0.25,K=5200,Binit=2900,sigma=0.20))
ans <- fitSPM(pars,fish,schaefer=TRUE,maxiter=1000,hess=TRUE)
# R-chunk 34 page 291
#The hessian matrix from the Schaefer fit to the dataspm data
outfit(ans)
# R-chunk 35 Page 292
#calculate the var-covar matrix and the st errors
vcov <- solve(ans$hessian) # calculate variance-covariance matrix
label <- c("r","K", "Binit","sigma")
colnames(vcov) <- label; rownames(vcov) <- label
outvcov <- rbind(vcov,sqrt(diag(vcov)))
rownames(outvcov) <- c(label,"StErr")
# R-chunk 36 Page 290 Table 7.6 code not in the book
# tabulate the variance covariance matrix and StErrs
kable(outvcov,digits=c(5,5,5,5))
# R-chunk 37 Pages 292 - 293
#generate 1000 parameter vectors from multi-variate normal
library(mvtnorm) # use RStudio, or install.packages("mvtnorm")
N <- 1000 # number of parameter vectors, use vcov from above
mvn <- length(fish[,"year"]) #matrix to store cpue trajectories
mvncpue <- matrix(0,nrow=N,ncol=mvn,dimnames=list(1:N,fish[,"year"]))
columns <- c("r","K","Binit","sigma")
optpar <- ans$estimate # Fill matrix with mvn parameter vectors
mvnpar <- matrix(exp(rmvnorm(N,mean=optpar,sigma=vcov)),nrow=N,
ncol=4,dimnames=list(1:N,columns))
msy <- mvnpar[,"r"]*mvnpar[,"K"]/4
nyr <- length(fish[,"year"])
depletion <- numeric(N) #now calculate N cpue series in linear space
for (i in 1:N) { # calculate dynamics for each parameter set
dynamA <- spm(log(mvnpar[i,1:4]),fish)
mvncpue[i,] <- dynamA$outmat[1:nyr,"predCE"]
depletion[i] <- dynamA$outmat["2016","Depletion"]
}
mvnpar <- cbind(mvnpar,msy,depletion) # try head(mvnpar,10)
# R-chunk 38 Page 293
#data and trajectories from 1000 MVN parameter vectors Fig 7.16
oldp <- plot1(fish[,"year"],fish[,"cpue"],type="p",xlab="Year",
ylab="CPUE",maxy=2.0)
for (i in 1:N) lines(fish[,"year"],mvncpue[i,],col="grey",lwd=1)
points(fish[,"year"],fish[,"cpue"],pch=1,cex=1.3,col=1,lwd=2) # data
lines(fish[,"year"],exp(simpspm(optpar,fish)),lwd=2,col=1)# pred
percs <- apply(mvncpue,2,quants) # obtain the quantiles
arrows(x0=fish[,"year"],y0=percs["5%",],y1=percs["95%",],length=0.03,
angle=90,code=3,col=1) #add 90% quantiles
msy <- mvnpar[,"r"]*mvnpar[,"K"]/4 # 1000 MSY estimates
text(2010,1.75,paste0("MSY ",round(mean(msy),3)),cex=1.25,font=7)
par(oldp) # return par to old settings; this line not in book
# R-chunk 39 Pages 293 - 294
#Isolate errant cpue trajectories Fig 7.17
pickd <- which(mvncpue[,"2016"] < 0.40)
oldp <- plot1(fish[,"year"],fish[,"cpue"],type="n",xlab="Year",
ylab="CPUE",maxy=6.25)
for (i in 1:length(pickd))
lines(fish[,"year"],mvncpue[pickd[i],],col=1,lwd=1)
points(fish[,"year"],fish[,"cpue"],pch=16,cex=1.25,col=4)
lines(fish[,"year"],exp(simpspm(optpar,fish)),lwd=3,col=2,lty=2)
par(oldp) # return par to old settings; this line not in book
# R-chunk 40 Page 294
#Use adhoc function to plot errant parameters Fig 7.18
oldp <- parset(plots=c(2,2),cex=0.85)
outplot <- function(var1,var2,pickdev) {
plot1(mvnpar[,var1],mvnpar[,var2],type="p",pch=16,cex=1.0,
defpar=FALSE,xlab=var1,ylab=var2,col=8)
points(mvnpar[pickdev,var1],mvnpar[pickdev,var2],pch=16,cex=1.0)
}
outplot("r","K",pickd) # assumes mvnpar in working environment
outplot("sigma","Binit",pickd)
outplot("r","Binit",pickd)
outplot("K","Binit",pickd)
par(oldp) # return par to old settings; this line not in book
# R-chunk 41 Page 296
#asymptotically sampled parameter vectors Fig 7.19
pairs(mvnpar,lower.panel=panel.smooth, upper.panel=panel.cor,
gap=0,cex=0.25,lwd=2)
# R-chunk 42 Page 297
# Get the ranges of parameters from bootstrap and asymptotic
bt <- apply(bootpar,2,range)[,c(1:4,6,7)]
ay <- apply(mvnpar,2,range)
out <- rbind(bt,ay)
rownames(out) <- c("MinBoot","MaxBoot","MinAsym","MaxAsym")
# R-chunk 43 Page 297 Table 7.7 code not in the book
#tabulate ranges from two approsches
kable(out,digits=c(4,3,3,4,3,4))
### Sometimes Asymptotic Errors Work
# R-chunk 44 Pages 297 - 298
#repeat asymptotice errors using abdat data-set Figure 7.20
data(abdat)
fish <- as.matrix(abdat)
pars <- log(c(r=0.4,K=9400,Binit=3400,sigma=0.05))
ansA <- fitSPM(pars,fish,schaefer=TRUE,maxiter=1000,hess=TRUE)
vcovA <- solve(ansA$hessian) # calculate var-covar matrix
mvn <- length(fish[,"year"])
N <- 1000 # replicates
mvncpueA <- matrix(0,nrow=N,ncol=mvn,dimnames=list(1:N,fish[,"year"]))
columns <- c("r","K","Binit","sigma")
optparA <- ansA$estimate # Fill matrix of parameter vectors
mvnparA <- matrix(exp(rmvnorm(N,mean=optparA,sigma=vcovA)),
nrow=N,ncol=4,dimnames=list(1:N,columns))
msy <- mvnparA[,"r"]*mvnparA[,"K"]/4
for (i in 1:N) mvncpueA[i,]<-exp(simpspm(log(mvnparA[i,]),fish))
mvnparA <- cbind(mvnparA,msy)
oldp <- plot1(fish[,"year"],fish[,"cpue"],type="p",xlab="Year",
ylab="CPUE",maxy=2.5)
for (i in 1:N) lines(fish[,"year"],mvncpueA[i,],col=8,lwd=1)
points(fish[,"year"],fish[,"cpue"],pch=16,cex=1.0) #orig data
lines(fish[,"year"],exp(simpspm(optparA,fish)),lwd=2,col=0)
par(oldp) # return par to old settings; this line not in book
# R-chunk 45 Page 298
#plot asymptotically sampled parameter vectors Figure 7.21
pairs(mvnparA,lower.panel=panel.smooth, upper.panel=panel.cor,
gap=0,pch=16,col=rgb(red=0,green=0,blue=0,alpha = 1/10))
### Bayesian Posteriors
# R-chunk 46 Page 299
#Fit the Fox Model to the abdat data Figure 7.22
data(abdat); fish <- as.matrix(abdat)
param <- log(c(r=0.3,K=11500,Binit=3300,sigma=0.05))
foxmod <- nlm(f=negLL1,p=param,funk=simpspm,indat=fish,
logobs=log(fish[,"cpue"]),iterlim=1000,schaefer=FALSE)
optpar <- exp(foxmod$estimate)
ans <- plotspmmod(inp=foxmod$estimate,indat=fish,schaefer=FALSE,
addrmse=TRUE, plotprod=TRUE)
# R-chunk 47 Page 301
# Conduct an MCMC using simpspmC on the abdat Fox SPM
# This means you will need to compile simpspmC from appendix
set.seed(698381) #for repeatability, possibly only on Windows10
begin <- gettime() # to enable the time taken to be calculated
inscale <- c(0.07,0.05,0.09,0.45) #note large value for sigma
pars <- log(c(r=0.205,K=11300,Binit=3200,sigma=0.044))
result <- do_MCMC(chains=1,burnin=50,N=2000,thinstep=512,
inpar=pars,infunk=negLL,calcpred=simpspmC,
obsdat=log(fish[,"cpue"]),calcdat=fish,
priorcalc=calcprior,scales=inscale,schaefer=FALSE)
# alternatively, use simpspm, but that will take longer.
cat("acceptance rate = ",result$arate," \n")
cat("time = ",gettime() - begin,"\n")
post1 <- result[[1]][[1]]
p <- 1e-08
msy <- post1[,"r"]*post1[,"K"]/((p + 1)^((p+1)/p))
# R-chunk 48 Page 302
#pairwise comparison for MCMC of Fox model on abdat Fig 7.23
pairs(cbind(post1[,1:4],msy),upper.panel = panel.cor,lwd=2,cex=0.2,
lower.panel=panel.smooth,col=1,gap=0.1)
# R-chunk 49 Page 302
# marginal distributions of 3 parameters and msy Figure 7.24
oldp <- parset(plots=c(2,2), cex=0.85)
plot(density(post1[,"r"]),lwd=2,main="",xlab="r") #plot has a method
plot(density(post1[,"K"]),lwd=2,main="",xlab="K") #for output from
plot(density(post1[,"Binit"]),lwd=2,main="",xlab="Binit") # density
plot(density(msy),lwd=2,main="",xlab="MSY") #try str(density(msy))
par(oldp) # return par to old settings; this line not in book
# R-chunk 50 Page 304
#MCMC r and K parameters, approx 50 + 90% contours. Fig7.25
puttxt <- function(xs,xvar,ys,yvar,lvar,lab="",sigd=0) {
text(xs*xvar[2],ys*yvar[2],makelabel(lab,lvar,sep=" ",
sigdig=sigd),cex=1.2,font=7,pos=4)
} # end of puttxt - a quick utility function
kran <- range(post1[,"K"]); rran <- range(post1[,"r"])
mran <- range(msy) #ranges used in the plots
oldp <- parset(plots=c(1,2),margin=c(0.35,0.35,0.05,0.1)) #plot r vs K
plot(post1[,"K"],post1[,"r"],type="p",cex=0.5,xlim=kran,
ylim=rran,col="grey",xlab="K",ylab="r",panel.first=grid())
points(optpar[2],optpar[1],pch=16,col=1,cex=1.75) # center
addcontours(post1[,"K"],post1[,"r"],kran,rran, #if fails make
contval=c(0.5,0.9),lwd=2,col=1) #contval smaller
puttxt(0.7,kran,0.97,rran,kran,"K= ",sigd=0)
puttxt(0.7,kran,0.94,rran,rran,"r= ",sigd=4)
plot(post1[,"K"],msy,type="p",cex=0.5,xlim=kran, # K vs msy
ylim=mran,col="grey",xlab="K",ylab="MSY",panel.first=grid())
points(optpar[2],getMSY(optpar,p),pch=16,col=1,cex=1.75)#center
addcontours(post1[,"K"],msy,kran,mran,contval=c(0.5,0.9),lwd=2,col=1)
puttxt(0.6,kran,0.99,mran,kran,"K= ",sigd=0)
puttxt(0.6,kran,0.97,mran,mran,"MSY= ",sigd=3)
par(oldp) # return par to old settings; this line not in book
# R-chunk 51 Page 305
#Traces for the Fox model parameters from the MCMC Fig7.26
oldp <- parset(plots=c(4,1),margin=c(0.3,0.45,0.05,0.05),
outmargin = c(1,0,0,0),cex=0.85)
label <- colnames(post1)
N <- dim(post1)[1]
for (i in 1:3) {
plot(1:N,post1[,i],type="l",lwd=1,ylab=label[i],xlab="")
abline(h=median(post1[,i]),col=2)
}
msy <- post1[,1]*post1[,2]/4
plot(1:N,msy,type="l",lwd=1,ylab="MSY",xlab="")
abline(h=median(msy),col=2)
mtext("Step",side=1,outer=T,line=0.0,font=7,cex=1.1)
par(oldp) # return par to old settings; this line not in book
# R-chunk 52 Page 306
#Do five chains of the same length for the Fox model
set.seed(6396679) # Note all chains start from same place, which is
inscale <- c(0.07,0.05,0.09,0.45) # suboptimal, but still the chains
pars <- log(c(r=0.205,K=11300,Binit=3220,sigma=0.044)) # differ
result <- do_MCMC(chains=5,burnin=50,N=2000,thinstep=512,
inpar=pars,infunk=negLL1,calcpred=simpspmC,
obsdat=log(fish[,"cpue"]),calcdat=fish,
priorcalc=calcprior,scales=inscale,
schaefer=FALSE)
cat("acceptance rate = ",result$arate," \n") # always check this
# R-chunk 53 Page 306
#Now plot marginal posteriors from 5 Fox model chains Fig7.27
oldp <- parset(plots=c(2,1),cex=0.85,margin=c(0.4,0.4,0.05,0.05))
post <- result[[1]][[1]]
plot(density(post[,"K"]),lwd=2,col=1,main="",xlab="K",
ylim=c(0,4.4e-04),panel.first=grid())
for (i in 2:5) lines(density(result$result[[i]][,"K"]),lwd=2,col=i)
p <- 1e-08
post <- result$result[[1]]
msy <- post[,"r"]*post[,"K"]/((p + 1)^((p+1)/p))
plot(density(msy),lwd=2,col=1,main="",xlab="MSY",type="l",
ylim=c(0,0.0175),panel.first=grid())
for (i in 2:5) {
post <- result$result[[i]]
msy <- post[,"r"]*post[,"K"]/((p + 1)^((p+1)/p))
lines(density(msy),lwd=2,col=i)
}
par(oldp) # return par to old settings; this line not in book
# R-chunk 54 Page 307
# get quantiles of each chain
probs <- c(0.025,0.05,0.5,0.95,0.975)
storeQ <- matrix(0,nrow=6,ncol=5,dimnames=list(1:6,probs))
for (i in 1:5) storeQ[i,] <- quants(result$result[[i]][,"K"])
x <- apply(storeQ[1:5,],2,range)
storeQ[6,] <- 100*(x[2,] - x[1,])/x[2,]
# R-chunk 55 Page 308 Table 7.8 code not in the book
#tabulate qunatiles of the five chains
kable(storeQ,digits=c(3,3,3,3,3))
## Management Advice
### Two Views of Risk
### Harvest Strategies
## Risk Assessment Projections
### Deterministic Projections
# R-chunk 56 Pages 310 - 311
#Prepare Fox model on abdat data for future projections Fig7.28
data(abdat); fish <- as.matrix(abdat)
param <- log(c(r=0.3,K=11500,Binit=3300,sigma=0.05))
bestmod <- nlm(f=negLL1,p=param,funk=simpspm,schaefer=FALSE,
logobs=log(fish[,"cpue"]),indat=fish,hessian=TRUE)
optpar <- exp(bestmod$estimate)
ans <- plotspmmod(inp=bestmod$estimate,indat=fish,schaefer=FALSE,
target=0.4,addrmse=TRUE, plotprod=FALSE)
# R-chunk 57 Page 312
out <- spm(bestmod$estimate,indat=fish,schaefer=FALSE)
str(out, width=65, strict.width="cut")
# R-chunk 58 Page 312 Table 7.9 code not in the book
#
kable(out$outmat[1:10,],digits=c(0,4,4,4,4,4,4))
# R-chunk 59 Page 313
# Fig 7.29
catches <- seq(700,1000,50) # projyr=10 is the default
projans <- spmprojDet(spmobj=out,projcatch=catches,plotout=TRUE)
### Accounting for Uncertainty
### Using Asymptotic Errors
# R-chunk 60 Page 315
# generate parameter vectors from a multivariate normal
# project dynamics under a constant catch of 900t
library(mvtnorm)
matpar <- parasympt(bestmod,N=1000) #generate parameter vectors
projs <- spmproj(matpar,fish,projyr=10,constC=900)#do dynamics
# R-chunk 61 Page 315
# Fig 7.30 1000 replicate projections asymptotic errors
outp <- plotproj(projs,out,qprob=c(0.1,0.5),refpts=c(0.2,0.4))
### Using Bootstrap Parameter Vectors
# R-chunk 62 Page 316
#bootstrap generation of plausible parameter vectors for Fox
reps <- 1000
boots <- spmboot(bestmod$estimate,fishery=fish,iter=reps,schaefer=FALSE)
matparb <- boots$bootpar[,1:4] #examine using head(matparb,20)
# R-chunk 63 Page 316
#bootstrap projections. Lower case b for boostrap Fig7.31
projb <- spmproj(matparb,fish,projyr=10,constC=900)
outb <- plotproj(projb,out,qprob=c(0.1,0.5),refpts=c(0.2,0.4))
### Using Samples from a Bayesian Posterior
# R-chunk 64 Pages 317 - 318
#Generate 1000 parameter vectors from Bayesian posterior
param <- log(c(r=0.3,K=11500,Binit=3300,sigma=0.05))
set.seed(444608)
N <- 1000
result <- do_MCMC(chains=1,burnin=100,N=N,thinstep=2048,
inpar=param,infunk=negLL,calcpred=simpspmC,
calcdat=fish,obsdat=log(fish[,"cpue"]),
priorcalc=calcprior,schaefer=FALSE,
scales=c(0.065,0.055,0.1,0.475))
parB <- result[[1]][[1]] #capital B for Bayesian
cat("Acceptance Rate = ",result[[2]],"\n")
# R-chunk 65 Page 318
# auto-correlation, or lack of, and the K trace Fig 7.32
oldp <- parset(plots=c(2,1),cex=0.85)
acf(parB[,2],lwd=2)
plot(1:N,parB[,2],type="l",ylab="K",ylim=c(8000,19000),xlab="")
par(oldp) # return par to old settings; this line not in book
# R-chunk 66 Page 318
# Fig 7.33
matparB <- as.matrix(parB[,1:4]) # B for Bayesian
projs <- spmproj(matparB,fish,constC=900,projyr=10) # project them
plotproj(projs,out,qprob=c(0.1,0.5),refpts=c(0.2,0.4)) #projections
## Concluding Remarks
## Appendix: The Use of Rcpp to Replace simpspm
# R-chunk 67 Page 321
library(Rcpp)
cppFunction('NumericVector simpspmC(NumericVector pars,
NumericMatrix indat, LogicalVector schaefer) {
int nyrs = indat.nrow();
NumericVector predce(nyrs);
NumericVector biom(nyrs+1);
double Bt, qval;
double sumq = 0.0;
double p = 0.00000001;
if (schaefer(0) == TRUE) {
p = 1.0;
}
NumericVector ep = exp(pars);
biom[0] = ep[2];
for (int i = 0; i < nyrs; i++) {
Bt = biom[i];
biom[(i+1)] = Bt + (ep[0]/p)*Bt*(1 - pow((Bt/ep[1]),p)) -
indat(i,1);
if (biom[(i+1)] < 40.0) biom[(i+1)] = 40.0;
sumq += log(indat(i,2)/biom[i]);
}
qval = exp(sumq/nyrs);
for (int i = 0; i < nyrs; i++) {
predce[i] = log(biom[i] * qval);
}
return predce;
}')
## End(Not run)
[Package MQMF version 0.1.5 Index]