| chapter3 {MQMF} | R Documentation |
chapter3 The 27 R-code chunks from Simple Population Models
Description
chapter3 is not an active function but rather acts as a repository for the various example code chunks found in chapter3. There are 27 r-code chunks in chapter3.
Usage
chapter3(verbose = TRUE)
Arguments
verbose |
Should instructions to written to the console, default=TRUE |
Examples
## Not run:
### The Discrete Logistic Model
# R-chunk 1 page 36
# Code to produce Figure 3.1. Note the two one-line functions
surprod <- function(Nt,r,K) return((r*Nt)*(1-(Nt/K)))
densdep <- function(Nt,K) return((1-(Nt/K)))
r <- 1.2; K <- 1000.0; Nt <- seq(10,1000,10)
oldpar <- par(no.readonly=TRUE) # this line not in book
# plotprep(width=7, height=5, newdev=FALSE)
par(mfrow=c(2,1),mai=c(0.4,0.4,0.05,0.05),oma=c(0.0,0,0.0,0.0))
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plot1(Nt,surprod(Nt,r,K),xlab="Population Nt",defpar=FALSE,
ylab="Production")
plot1(Nt,densdep(Nt,K),xlab="Population Nt",defpar=FALSE,
ylab="Density-Dependence")
par(oldpar) # this line not in book
### Dynamic Behaviour
# R-chunk 2 page 38
#Code for Figure 3.2. Try varying the value of rv from 0.5-2.8
yrs <- 100; rv=2.8; Kv <- 1000.0; Nz=100; catch=0.0; p=1.0
ans <- discretelogistic(r=rv,K=Kv,N0=Nz,Ct=catch,Yrs=yrs,p=p)
avcatch <- mean(ans[(yrs-50):yrs,"nt"],na.rm=TRUE) #used in text
label <- paste0("r=",rv," K=",Kv," Ct=",catch, " N0=",Nz," p=",p=p)
oldpar <- par(no.readonly=TRUE) # this line not in book
plot(ans, main=label, cex=0.9, font=7) #Schaefer dynamics
par(oldpar) # this line not in book
# R-chunk 3 page 39
#run discrete logistic dynamics for 600 years
yrs=600
ans <- discretelogistic(r=2.55,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)
# R-chunk 4 page 40, code not in the book
#tabulate the last 30 years of the dynamics needs knitr
library(knitr)
kable(halftable(ans[(yrs-29):yrs,],yearcol="year",subdiv=3),digits=c(0,1,1,0,1,1,0,1,1))
### Finding Boundaries between Behaviours.
# R-chunk 5 page 40
#run discretelogistic and search for repeated values of Nt
yrs <- 600
ans <- discretelogistic(r=2.55,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)
avt <- round(apply(ans[(yrs-100):(yrs-1),2:3],1,mean),2)
count <- table(avt)
count[count > 1] # with r=2.55 you should find an 8-cycle limit
# R-chunk 6 page 41
#searches for unique solutions given an r value see Table 3.2
testseq <- seq(1.9,2.59,0.01)
nseq <- length(testseq)
result <- matrix(0,nrow=nseq,ncol=2,
dimnames=list(testseq,c("r","Unique")))
yrs <- 600
for (i in 1:nseq) { # i = 31
rval <- testseq[i]
ans <- discretelogistic(r=rval,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)
ans <- ans[-yrs,] # remove last year, see str(ans) for why
ans[,"nt1"] <- round(ans[,"nt1"],3) #try hashing this out
result[i,] <- c(rval,length(unique(tail(ans[,"nt1"],100))))
}
# R-chunk 7 page 41 - 42, Table 3.2. Code not in the book.
#unique repeated Nt values 100 = non-equilibrium or chaos
kable(halftable(result,yearcol = "r"),)
### Classical Bifurcation Diagram of Chaos
# R-chunk 8 pages 42 - 43
#the R code for the bifurcation function
bifurcation <- function(testseq,taill=100,yrs=1000,limy=0,incx=0.001){
nseq <- length(testseq)
result <- matrix(0,nrow=nseq,ncol=2,
dimnames=list(testseq,c("r","Unique Values")))
result2 <- matrix(NA,nrow=nseq,ncol=taill)
for (i in 1:nseq) {
rval <- testseq[i]
ans <- discretelogistic(r=rval,K=1000.0,N0=100,Ct=0.0,Yrs=yrs)
ans[,"nt1"] <- round(ans[,"nt1"],4)
result[i,] <- c(rval,length(unique(tail(ans[,"nt1"],taill))))
result2[i,] <- tail(ans[,"nt1"],taill)
}
if (limy[1] == 0) limy <- c(0,getmax(result2,mult=1.02))
oldpar <- parset() #plot taill values against taill of each r value
on.exit(par(oldpar)) # this line not in book
plot(rep(testseq[1],taill),result2[1,],type="p",pch=16,cex=0.1,
ylim=limy,xlim=c(min(testseq)*(1-incx),max(testseq)*(1+incx)),
xlab="r value",yaxs="i",xaxs="i",ylab="Equilibrium Numbers",
panel.first=grid())
for (i in 2:nseq)
points(rep(testseq[i],taill),result2[i,],pch=16,cex=0.1)
return(invisible(list(result=result,result2=result2)))
} # end of bifurcation
# R-chunk 9 page 43
#Alternative r value arrangements for you to try; Fig 3.3
#testseq <- seq(2.847,2.855,0.00001) #hash/unhash as needed
#bifurcation(testseq,limy=c(600,740),incx=0.0001) # t
#testseq <- seq(2.6225,2.6375,0.00001) # then explore
#bifurcation(testseq,limy=c(660,730),incx=0.0001)
testseq <- seq(1.9,2.975,0.0005) # modify to explore
bifurcation(testseq,limy=0)
### The Effect of Fishing on Dynamics
# R-chunk 10 page 43 - 44.
#Effect of catches on stability properties of discretelogistic
yrs=50; Kval=1000.0
nocatch <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=0,Yrs=yrs)
catch50 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=50,Yrs=yrs)
catch200 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=200,Yrs=yrs)
catch300 <- discretelogistic(r=2.56,K=Kval,N0=500,Ct=300,Yrs=yrs)
# R-chunk 11 page 45
#Effect of different catches on n-cyclic behaviour Fig3.4
plottime <- function(x,ylab) {
yrs <- nrow(x)
plot1(x[,"year"],x[,"nt"],ylab=ylab,defpar=FALSE)
avB <- round(mean(x[(yrs-40):yrs,"nt"],na.rm=TRUE),3)
mtext(avB,side=1,outer=F,line=-1.1,font=7,cex=1.0)
} # end of plottime
#the oma argument is used to adjust the space around the graph
oldpar <- par(no.readonly=TRUE) # this line not in book
par(mfrow=c(2,2),mai=c(0.25,0.4,0.05,0.05),oma=c(1.0,0,0.25,0))
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plottime(nocatch,"Catch = 0")
plottime(catch50,"Catch = 50")
plottime(catch200,"Catch = 200")
plottime(catch300,"Catch = 300")
mtext("years",side=1,outer=TRUE,line=-0.2,font=7,cex=1.0)
par(oldpar)
# R-chunk 12 page 46
#Phase plot for Schaefer model Fig 3.5
plotphase <- function(x,label,ymax=0) { #x from discretelogistic
yrs <- nrow(x)
colnames(x) <- tolower(colnames(x))
if (ymax[1] == 0) ymax <- getmax(x[,c(2:3)])
plot(x[,"nt"],x[,"nt1"],type="p",pch=16,cex=1.0,ylim=c(0,ymax),
yaxs="i",xlim=c(0,ymax),xaxs="i",ylab="nt1",xlab="",
panel.first=grid(),col="darkgrey")
begin <- trunc(yrs * 0.6) #last 40% of yrs = 20, when yrs=50
points(x[begin:yrs,"nt"],x[begin:yrs,"nt1"],pch=18,col=1,cex=1.2)
mtext(label,side=1,outer=F,line=-1.1,font=7,cex=1.2)
} # end of plotphase
oldpar <- par(no.readonly=TRUE) # this line not in book
par(mfrow=c(2,2),mai=c(0.25,0.25,0.05,0.05),oma=c(1.0,1.0,0,0))
par(cex=0.75, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7)
plotphase(nocatch,"Catch = 0",ymax=1300)
plotphase(catch50,"Catch = 50",ymax=1300)
plotphase(catch200,"Catch = 200",ymax=1300)
plotphase(catch300,"Catch = 300",ymax=1300)
mtext("nt",side=1,outer=T,line=0.0,font=7,cex=1.0)
mtext("nt+1",side=2,outer=T,line=0.0,font=7,cex=1.0)
par(oldpar) # this line not in book
### Determinism
## Age-Structured Modelling Concepts
### Survivorship in a Cohort
# R-chunk 13 pages 48 - 49
#Exponential population declines under different Z. Fig 3.6
yrs <- 50; yrs1 <- yrs + 1 # to leave room for B[0]
years <- seq(0,yrs,1)
B0 <- 1000 # now alternative total mortality rates
Z <- c(0.05,0.1,0.2,0.4,0.55)
nZ <- length(Z)
Bt <- matrix(0,nrow=yrs1,ncol=nZ,dimnames=list(years,Z))
Bt[1,] <- B0
for (j in 1:nZ) for (i in 2:yrs1) Bt[i,j] <- Bt[(i-1),j]*exp(-Z[j])
oldp <- plot1(years,Bt[,1],xlab="Years",ylab="Population Size",lwd=2)
if (nZ > 1) for (j in 2:nZ) lines(years,Bt[,j],lwd=2,col=j,lty=j)
legend("topright",legend=paste0("Z = ",Z),col=1:nZ,lwd=3,
bty="n",cex=1,lty=1:5)
par(oldp) # this line not in book
### Instantaneous vs Annual Mortality Rates
# R-chunk 14 page 51
#Prepare matrix of harvest rate vs time to appoximate F
Z <- -log(0.5)
timediv <- c(2,4,12,52,365,730,2920,8760,525600)
yrfrac <- 1/timediv
names(yrfrac) <- c("6mth","3mth","1mth","1wk","1d","12h","3h","1h","1m")
nfrac <- length(yrfrac)
columns <- c("yrfrac","divisor","yrfracH","Remain")
result <- matrix(0,nrow=nfrac,ncol=length(columns),
dimnames=list(names(yrfrac),columns))
for (i in 1:nfrac) {
timestepmort <- Z/timediv[i]
N <- 1000
for (j in 1:timediv[i]) N <- N * (1-timestepmort)
result[i,] <- c(yrfrac[i],timediv[i],timestepmort,N)
}
# R-chunk 15 page 51 Table 3.3, code not shown in book
#output of constant Z for shorter and shorter periods
kable(result,digits=c(10,0,8,4))
# R-chunk 16 page 51
#Annual harvest rate against instantaneous F, Fig 3.7
Fi <- seq(0.001,2,0.001)
H <- 1 - exp(-Fi)
oldpar <- parset() # a wrapper for simplifying defining the par values
plot(Fi,H,type="l",lwd=2,panel.first=grid(),xlab="Instantaneous Fishing Mortality F",
ylab="Annual Proportion Mortality H")
lines(c(0,1),c(0,1),lwd=2,lty=2,col=2)
par(oldpar) # this line not in book
## Simple Yield per Recruit
# R-chunk 17 page 53
# Simple Yield-per-Recruit see Russell (1942)
age <- 1:11; nage <- length(age); N0 <- 1000 # some definitions
# weight-at-age values
WaA <- c(NA,0.082,0.175,0.283,0.4,0.523,0.7,0.85,0.925,0.99,1.0)
# now the harvest rates
H <- c(0.01,0.06,0.11,0.16,0.21,0.26,0.31,0.36,0.55,0.8)
nH <- length(H)
NaA <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,H)) # storage
CatchN <- NaA; CatchW <- NaA # define some storage matrices
for (i in 1:nH) { # loop through the harvest rates
NaA[1,i] <- N0 # start each harvest rate with initial numbers
for (age in 2:nage) { # loop through over-simplified dynamics
NaA[age,i] <- NaA[(age-1),i] * (1 - H[i])
CatchN[age,i] <- NaA[(age-1),i] - NaA[age,i]
}
CatchW[,i] <- CatchN[,i] * WaA
} # transpose the vector of total catches to
totC <- t(colSums(CatchW,na.rm=TRUE)) # simplify later printing
# R-chunk 18 page 54 Table 3.4 code not shown in book
#Tabulate numbers-at-age for different harvest rates needs knitr
kable(NaA,digits=c(0,0,0,0,0,0,0,0,1,1),row.names=TRUE)
# R-chunk 19 page 54, Table 3.5, code not shown in book.
#Tabulate Weight-at-age for different harvest rates
kable(CatchW[2:11,],digits=c(2,2,2,2,2,2,2,2,2,2),row.names=TRUE)
# R-chunk 20 page 54, Table 3.6, code not shown in book.
#Total weights vs Harvest rate
kable(totC,digits=c(1,1,1,1,1,1,1,1,1,1))
# R-chunk 21 page 55
#Use MQMF::plot1 for a quick plot of the total catches. Figure 3.8
oldpar <- plot1(H,totC,xlab="Harvest Rate",ylab="Total Yield",lwd=2)
par(oldpar) # to reset the par values if desired
### Selectivity in Yield-per-Recruit
# R-chunk 22 Page 56
#Logistic S shaped cureve for maturity
ages <- seq(0,50,1)
sel1 <- mature(-3.650425,0.146017,sizeage=ages) #-3.65/0.146=25
sel2 <- mature(-6,0.2,ages)
sel3 <- mature(-6,0.24,ages)
oldp <- plot1(ages,sel1,xlab="Age Yrs",ylab="Selectivity",cex=0.75,lwd=2)
lines(ages,sel2,col=2,lwd=2,lty=2)
lines(ages,sel3,col=3,lwd=2,lty=3)
abline(v=25,col="grey",lty=2)
abline(h=c(0.25,0.5,0.75),col="grey",lty=2)
legend("topleft",c("25_15.04","30_10.986","25_9.155"),col=c(1,2,3),
lwd=3,cex=1.1,bty="n",lty=1:3)
par(oldp)
### The Baranov Catch Equation
# R-chunk 23 Page 58
# Baranov catch equation
age <- 0:12; nage <- length(age)
sa <-mature(-4,2,age) #selectivity-at-age
H <- 0.2; M <- 0.35
FF <- -log(1 - H)#Fully selected instantaneous fishing mortality
Ft <- sa * FF # instantaneous Fishing mortality-at-age
N0 <- 1000
out <- cbind(bce(M,Ft,N0,age),"Select"=sa) # out becomes Table 3.7
# R-chunk 24 page 59, Table 3.7, code not shown in book.
#tabulate output from Baranov Catch Equations
kable(out,digits=c(3,3,3,3))
### Growth and Weight-at-Age
## Full Yield-per-Recruit
# R-chunk 25 Page 60 - 61
# A more complete YPR analysis
age <- 0:20; nage <- length(age) #storage vectors and matrices
laa <- vB(c(50.0,0.25,-1.5),age) # length-at-age
WaA <- (0.015 * laa ^ 3.0)/1000 # weight-at-age as kg
H <- seq(0.01,0.65,0.05); nH <- length(H)
FF <- round(-log(1 - H),5) # Fully selected fishing mortality
N0 <- 1000
M <- 0.1
numt <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,FF))
catchN <- matrix(0,nrow=nage,ncol=nH,dimnames=list(age,FF))
as50 <- c(1,2,3)
yield <- matrix(0,nrow=nH,ncol=length(as50),dimnames=list(H,as50))
for (sel in 1:length(as50)) {
sa <- logist(as50[sel],1.0,age) # selectivity-at-age
for (harv in 1:nH) {
Ft <- sa * FF[harv] # Fishing mortality-at-age
out <- bce(M,Ft,N0,age)
numt[,harv] <- out[,"Nt"]
catchN[,harv] <- out[,"Catch"]
yield[harv,sel] <- sum(out[,"Catch"] * WaA,na.rm=TRUE)
} # end of harv loop
} # end of sel loop
# R-chunk 26 Page 61
#A full YPR analysis Figure 3.10
oldp <- plot1(H,yield[,3],xlab="Harvest Rate",ylab="Yield",cex=0.75,lwd=2)
lines(H,yield[,2],lwd=2,col=2,lty=2)
lines(H,yield[,1],lwd=2,col=3,lty=3)
legend("bottomright",legend=as50,col=c(3,2,1),lwd=3,bty="n",
cex=1.0,lty=c(3,2,1))
par(oldp)
# R-chunk 27 page 62, Table 3.8, code not shown in book.
#Tabulate yield-per-recruit using Baranoc catch equation
kable(yield,digits=c(2,3,3,3))
## End(Not run)
[Package MQMF version 0.1.5 Index]