| chapter5 {MQMF} | R Documentation |
chapter5 The 23 R-code chunks from Static Models
Description
chapter5 is not an active function but rather acts as a repository for the various example code chunks found in chapter5. There are 23 r-code chunks in chapter5. You should, of course, feel free to use and modify any of these example chunks in your own work.
Usage
chapter5(verbose = TRUE)
Arguments
verbose |
Should instructions to written to the console, default=TRUE |
Examples
## Not run:
# All the example code from # Static Models
# Static Models
## Introduction
## Productivity Parameters
## Growth
### Seasonal Growth Curves
# R-chunk 1 Page 138
#vB growth curve fit to Pitcher and Macdonald derived seasonal data
data(minnow); week <- minnow$week; length <- minnow$length
pars <- c(75,0.1,-10.0,3.5); label=c("Linf","K","t0","sigma")
bestvB <- nlm(f=negNLL,p=pars,funk=vB,ages=week,observed=length,
typsize=magnitude(pars))
predL <- vB(bestvB$estimate,0:160)
outfit(bestvB,backtran = FALSE,title="Non-Seasonal vB",parnames=label)
# R-chunk 2 Page 139
#plot the non-seasonal fit and its residuals. Figure 5.1
oldp <- parset(plots=c(2,1),margin=c(0.35,0.45,0.02,0.05))
plot1(week,length,type="p",cex=1.0,col=2,xlab="Weeks",pch=16,
ylab="Length (mm)",defpar=FALSE)
lines(0:160,predL,lwd=2,col=1)
# calculate and plot the residuals
resids <- length - vB(bestvB$estimate,week)
plot1(week,resids,type="l",col="darkgrey",cex=0.9,lwd=2,
xlab="Weeks",lty=3,ylab="Normal Residuals",defpar=FALSE)
points(week,resids,pch=16,cex=1.1,col="red")
abline(h=0,col=1,lwd=1)
par(oldp) # return par to old settings; this line not in book
# R-chunk 3 Pages 139 - 140
# Fit seasonal vB curve, parameters = Linf, K, t0, C, s, sigma
svb <- function(p,ages,inc=52) {
return(p[1]*(1 - exp(-(p[4] * sin(2*pi*(ages - p[5])/inc) +
p[2] * (ages - p[3])))))
} # end of svB
spars <- c(bestvB$estimate[1:3],0.1,5,2.0) # keep sigma at end
bestsvb <- nlm(f=negNLL,p=spars,funk=svb,ages=week,observed=length,
typsize=magnitude(spars))
predLs <- svb(bestsvb$estimate,0:160)
outfit(bestsvb,backtran = FALSE,title="Seasonal Growth",
parnames=c("Linf","K","t0","C","s","sigma"))
# R-chunk 4 Page 140
#Plot seasonal growth curve and residuals Figure 5.2
oldp <- parset(plots=c(2,1)) # MQMF utility wrapper function
plot1(week,length,type="p",cex=0.9,col=2,xlab="Weeks",pch=16,
ylab="Length (mm)",defpar=FALSE)
lines(0:160,predLs,lwd=2,col=1)
# calculate and plot the residuals
resids <- length - svb(bestsvb$estimate,week)
plot1(week,resids,type="l",col="darkgrey",cex=0.9,xlab="Weeks",
lty=3,ylab="Normal Residuals",defpar=FALSE)
points(week,resids,pch=16,cex=1.1,col="red")
abline(h=0,col=1,lwd=1)
par(oldp) # return par to old settings; this line not in book
### Fabens Method with Tagging Data
# R-chunk 5 Pages 142 - 143
# tagging growth increment data from Black Island, Tasmania
data(blackisland); bi <- blackisland # just to keep things brief
oldp <- parset()
plot(bi$l1,bi$dl,type="p",pch=16,cex=1.0,col=2,ylim=c(-1,33),
ylab="Growth Increment mm",xlab="Initial Length mm",
panel.first = grid())
abline(h=0,col=1)
par(oldp) # return par to old settings; this line not in book
### Fitting Models to Tagging Data
# R-chunk 6 Page 144
# Fit the vB and Inverse Logistic to the tagging data
linm <- lm(bi$dl ~ bi$l1) # simple linear regression
param <- c(170.0,0.3,4.0); label <- c("Linf","K","sigma")
modelvb <- nlm(f=negNLL,p=param,funk=fabens,observed=bi$dl,indat=bi,
initL="l1",delT="dt") # could have used the defaults
outfit(modelvb,backtran = FALSE,title="vB",parnames=label)
predvB <- fabens(modelvb$estimate,bi)
cat("\n")
param2 <- c(25.0,130.0,35.0,3.0)
label2=c("MaxDL","L50","delta","sigma")
modelil <- nlm(f=negNLL,p=param2,funk=invl,observed=bi$dl,indat=bi,
initL="l1",delT="dt")
outfit(modelil,backtran = FALSE,title="IL",parnames=label2)
predil <- invl(modelil$estimate,bi)
# R-chunk 7 Page 145
#growth curves and regression fitted to tagging data Fig 5.4
oldp <- parset(margin=c(0.4,0.4,0.05,0.05))
plot(bi$l1,bi$dl,type="p",pch=16,cex=1.0,col=3,ylim=c(-2,31),
ylab="Growth Increment mm",xlab="Length mm",panel.first=grid())
abline(h=0,col=1)
lines(bi$l1,predvB,pch=16,col=1,lwd=3,lty=1) # vB
lines(bi$l1,predil,pch=16,col=2,lwd=3,lty=2) # IL
abline(linm,lwd=3,col=7,lty=2) # add dashed linear regression
legend("topright",c("vB","LinReg","IL"),lwd=3,bty="n",cex=1.2,
col=c(1,7,2),lty=c(1,2,2))
par(oldp) # return par to old settings; this line not in book
# R-chunk 8 Pages 145 - 146
#residuals for vB and inverse logistic for tagging data Fig 5.5
oldp <- parset(plots=c(1,2),outmargin=c(1,1,0,0),margin=c(.25,.25,.05,.05))
plot(bi$l1,(bi$dl - predvB),type="p",pch=16,col=1,ylab="",
xlab="",panel.first=grid(),ylim=c(-8,11))
abline(h=0,col=1)
mtext("vB",side=1,outer=FALSE,line=-1.1,cex=1.2,font=7)
plot(bi$l1,(bi$dl - predil),type="p",pch=16,col=1,ylab="",
xlab="",panel.first=grid(),ylim=c(-8,11))
abline(h=0,col=1)
mtext("IL",side=3,outer=FALSE,line=-1.2,cex=1.2,font=7)
mtext("Length mm",side=1,line=-0.1,cex=1.0,font=7,outer=TRUE)
mtext("Residual",side=2,line=-0.1,cex=1.0,font=7,outer=TRUE)
par(oldp) # return par to old settings; this line not in book
### A Closer Look at the Fabens Methods
### Implementation of Non-Constant Variances
# R-chunk 9 Page 149
# fit the Fabens tag growth curve with and without the option to
# modify variation with predicted length. See the MQMF function
# negnormL. So first no variation and then linear variation.
sigfunk <- function(pars,predobs) return(tail(pars,1)) #no effect
data(blackisland)
bi <- blackisland # just to keep things brief
param <- c(170.0,0.3,4.0); label=c("Linf","K","sigma")
modelvb <- nlm(f=negnormL,p=param,funk=fabens,funksig=sigfunk,
indat=bi,initL="l1",delT="dt")
outfit(modelvb,backtran = FALSE,title="vB constant sigma",parnames = label)
sigfunk2 <- function(pars,predo) { # linear with predicted length
sig <- tail(pars,1) * predo # sigma x predDL, see negnormL
pick <- which(sig <= 0) # ensure no negative sigmas from
sig[pick] <- 0.01 # possible negative predicted lengths
return(sig)
} # end of sigfunk2
param <- c(170.0,0.3,1.0); label=c("Linf","K","sigma")
modelvb2 <- nlm(f=negnormL,p=param,funk=fabens,funksig=sigfunk2,
indat=bi,initL="l1",delT="dt",
typsize=magnitude(param),iterlim=200)
outfit(modelvb2,backtran = FALSE,parnames = label,title="vB inverse DeltaL, sigma < 1")
# R-chunk 10 Page 150
#plot to two Faben's lines with constant and varying sigma Fig 5.6
predvB <- fabens(modelvb$estimate,bi)
predvB2 <- fabens(modelvb2$estimate,bi)
oldp <- parset(margin=c(0.4,0.4,0.05,0.05))
plot(bi$l1,bi$dl,type="p",pch=1,cex=1.0,col=1,ylim=c(-2,31),
ylab="Growth Increment mm",xlab="Length mm",panel.first=grid())
abline(h=0,col=1)
lines(bi$l1,predvB,col=1,lwd=2) # vB
lines(bi$l1,predvB2,col=2,lwd=2,lty=2) # IL
legend("topright",c("Constant sigma","Changing sigma"),lwd=3,
col=c(1,2),bty="n",cex=1.1,lty=c(1,2))
par(oldp) # return par to old settings; this line not in book
## Objective Model Selection
### Akiake's Information Criterion
# R-chunk 11 Page 152
#compare the relative model fits of Vb and IL
cat("von Bertalanffy \n")
aicbic(modelvb,bi)
cat("inverse-logistic \n")
aicbic(modelil,bi)
### Likelihood Ratio Test
# R-chunk 12 Page 154
# Likelihood ratio comparison of two growth models see Fig 5.4
vb <- modelvb$minimum # their respective -ve log-likelihoods
il <- modelil$minimum
dof <- 1
round(likeratio(vb,il,dof),8)
### Caveats on Likelihood Ratio Tests
## Remarks on Growth
## Maturity
### Introduction
### Alternative Maturity Ogives
# R-chunk 13 Page 158
# The Maturity data from tasab data-set
data(tasab) # see ?tasab for a list of the codes used
properties(tasab) # summarize properties of columns in tasab
table(tasab$site,tasab$sex) # sites 1 & 2 vs F, I, and M
# R-chunk 14 Page 158
#plot the proportion mature vs shell length Fig 5.7
propm <- tapply(tasab$mature,tasab$length,mean) #mean maturity at L
lens <- as.numeric(names(propm)) # lengths in the data
oldp <- plot1(lens,propm,type="p",cex=0.9,xlab="Length mm",
ylab="Proportion Mature")
par(oldp) # return par to old settings; this line not in book
# R-chunk 15 Pages 159 - 160
#Use glm to estimate mature logistic
binglm <- function(x,digits=6) { #function to simplify printing
out <- summary(x)
print(out$call)
print(round(out$coefficients,digits))
cat("\nNull Deviance ",out$null.deviance,"df",out$df.null,"\n")
cat("Resid.Deviance ",out$deviance,"df",out$df.residual,"\n")
cat("AIC = ",out$aic,"\n\n")
return(invisible(out)) # retain the full summary
} #end of binglm
tasab$site <- as.factor(tasab$site) # site as a factor
smodel <- glm(mature ~ site + length,family=binomial,data=tasab)
outs <- binglm(smodel) #outs contains the whole summary object
model <- glm(mature ~ length, family=binomial, data=tasab)
outm <- binglm(model)
cof <- outm$coefficients
cat("Lm50 = ",-cof[1,1]/cof[2,1],"\n")
cat("IQ = ",2*log(3)/cof[2,1],"\n")
# R-chunk 16 Page 161
#Add maturity logistics to the maturity data plot Fig 5.8
propm <- tapply(tasab$mature,tasab$length,mean) #prop mature
lens <- as.numeric(names(propm)) # lengths in the data
pick <- which((lens > 79) & (lens < 146))
oldp <- parset()
plot(lens[pick],propm[pick],type="p",cex=0.9, #the data points
xlab="Length mm",ylab="Proportion Mature",pch=1)
L <- seq(80,145,1) # for increased curve separation
pars <- coef(smodel)
lines(L,mature(pars[1],pars[3],L),lwd=3,col=3,lty=2)
lines(L,mature(pars[1]+pars[2],pars[3],L),lwd=3,col=2,lty=4)
lines(L,mature(coef(model)[1],coef(model)[2],L),lwd=2,col=1,lty=1)
abline(h=c(0.25,0.5,0.75),lty=3,col="grey")
legend("topleft",c("site1","both","site2"),col=c(3,1,2),lty=c(2,1,4),
lwd=3,bty="n")
par(oldp) # return par to old settings; this line not in book
### The Assumption of Symmetry
# R-chunk 17 Page 163
#Asymmetrical maturity curve from Schnute and Richard's curve Fig5.9
L = seq(50,160,1)
p=c(a=0.07,b=0.2,c=1.0,alpha=100)
asym <- srug(p=p,sizeage=L)
L25 <- linter(bracket(0.25,asym,L))
L50 <- linter(bracket(0.5,asym,L))
L75 <- linter(bracket(0.75,asym,L))
oldp <- parset()
plot(L,asym,type="l",lwd=2,xlab="Length mm",ylab="Proportion Mature")
abline(h=c(0.25,0.5,0.75),lty=3,col="grey")
abline(v=c(L25,L50,L75),lwd=c(1,2,1),col=c(1,2,1))
par(oldp) # return par to old settings; this line not in book
# R-chunk 18 Page 164 code not printed in the book
#Variation possible using the Schnute and Richard's Curve fig 5.10
# This code not printed in the book
tmplot <- function(vals,label) {
text(170,0.6,paste0(" ",label),font=7,cex=1.5)
legend("bottomright",legend=vals,col=1:nvals,lwd=3,bty="n",
cex=1.25,lty=c(1:nvals))
}
L = seq(50,180,1)
vals <- seq(0.05,0.09,0.01) # a value
nvals <- length(vals)
asym <- srug(p=c(a=vals[1],b=0.2,c=1.0,alpha=100),sizeage=L)
oldp <- parset(plots=c(2,2))
plot(L,asym,type="l",lwd=2,xlab="Length mm",ylab="Proportion Mature",
ylim=c(0,1.05))
abline(h=c(0.25,0.5,0.75),lty=3,col="darkgrey")
for (i in 2:nvals) {
asym <- srug(p=c(a=vals[i],b=0.2,c=1.0,alpha=100),sizeage=L)
lines(L,asym,lwd=2,col=i,lty=i)
}
tmplot(vals,"a")
vals <- seq(0.02,0.34,0.08) # b value
nvals <- length(vals)
asym <- srug(p=c(a=0.07,b=vals[1],c=1.0,alpha=100),sizeage=L)
plot(L,asym,type="l",lwd=2,xlab="Length mm",ylab="Proportion Mature",
ylim=c(0,1.05))
abline(h=c(0.25,0.5,0.75),lty=3,col="darkgrey")
for (i in 2:nvals) {
asym <- srug(p=c(a=0.07,b=vals[i],c=1.0,alpha=100),sizeage=L)
lines(L,asym,lwd=2,col=i,lty=i)
}
tmplot(vals,"b")
vals <- seq(0.95,1.05,0.025) # c value
nvals <- length(vals)
asym <- srug(p=c(a=0.07,b=0.2,c=vals[1],alpha=100),sizeage=L)
plot(L,asym,type="l",lwd=2,xlab="Length mm",ylab="Proportion Mature",
ylim=c(0,1.05))
abline(h=c(0.25,0.5,0.75),lty=3,col="darkgrey")
for (i in 2:nvals) {
asym <- srug(p=c(a=0.07,b=0.2,c=vals[i],alpha=100),sizeage=L)
lines(L,asym,lwd=2,col=i,lty=i)
}
tmplot(vals,"c")
vals <- seq(25,225,50) # alpha value
nvals <- length(vals)
asym <- srug(p=c(a=0.07,b=0.2,c=1.0,alpha=vals[1]),sizeage=L)
plot(L,asym,type="l",lwd=2,xlab="Length mm",ylab="Proportion Mature",
ylim=c(0,1.05))
abline(h=c(0.25,0.5,0.75),lty=3,col="darkgrey")
for (i in 2:nvals) {
asym <- srug(p=c(a=0.07,b=0.2,c=1.0,alpha=vals[i]),sizeage=L)
lines(L,asym,lwd=2,col=i,lty=i)
}
tmplot(vals,"alpha")
par(oldp) # return par to old settings; this line not in book
## Recruitment
### Introduction
### Properties of Good Stock Recruitment Relationships
### Recruitment Overfishing
### Beverton and Holt Recruitment
# R-chunk 19 Page 169
#plot the MQMF bh function for Beverton-Holt recruitment Fig 5.11
B <- 1:3000
bhb <- c(1000,500,250,150,50)
oldp <- parset()
plot(B,bh(c(1000,bhb[1]),B),type="l",ylim=c(0,1050),
xlab="Spawning Biomass",ylab="Recruitment")
for (i in 2:5) lines(B,bh(c(1000,bhb[i]),B),lwd=2,col=i,lty=i)
legend("bottomright",legend=bhb,col=c(1:5),lwd=3,bty="n",lty=c(1:5))
abline(h=c(500,1000),col=1,lty=2)
par(oldp) # return par to old settings; this line not in book
### Ricker Recruitment
# R-chunk 20 Page 170
#plot the MQMF ricker function for Ricker recruitment Fig 5.12
B <- 1:20000
rickb <- c(0.0002,0.0003,0.0004)
oldp <- parset()
plot(B,ricker(c(10,rickb[1]),B),type="l",xlab="Spawning Biomass",ylab="Recruitment")
for (i in 2:3)
lines(B,ricker(c(10,rickb[i]),B),lwd=2,col=i,lty=i)
legend("topright",legend=rickb,col=1:3,lty=1:3,bty="n",lwd=2)
par(oldp) # return par to old settings; this line not in book
### Deriso's Generalized Model
# R-chunk 21 Page 172
# plot of three special cases from Deriso-Schnute curve Fig. 5.13
deriso <- function(p,B) return(p[1] * B *(1 - p[2]*p[3]*B)^(1/p[3]))
B <- 1:10000
oldp <- plot1(B,deriso(c(10,0.001,-1),B),lwd=2,xlab="Spawning Biomass",
ylab="Recruitment")
lines(B,deriso(c(10,0.0004,0.25),B),lwd=2,col=2,lty=2) # DS
lines(B,deriso(c(10,0.0004,1e-06),B),lwd=2,col=3,lty=3) # Ricker
lines(B,deriso(c(10,0.0004,0.5),B),lwd=2,col=1,lty=3) # odd line
legend(x=7000,y=8500,legend=c("BH","DS","Ricker","odd line"),
col=c(1,2,3,1),lty=c(1,2,3,3),bty="n",lwd=3)
par(oldp) # return par to old settings; this line not in book
### Re-Parameterized Beverton-Holt Equation
### Re-Parameterized Ricker Equation
## Selectivity
### Introduction
### Logistic Selection
# R-chunk 22 Page 177
#Selectivity curves from logist and mature functions See Fig 5.14
ages <- seq(0,50,1); in50 <- 25.0
sel1 <- logist(in50,12,ages) #-3.65/0.146=L50=25.0
sel2 <- mature(-3.650425,0.146017,sizeage=ages)
sel3 <- mature(-6,0.2,ages)
sel4 <- logist(22.0,14,ages,knifeedge = TRUE)
oldp <- plot1(ages,sel1,xlab="Age Years",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)
lines(ages,sel4,col=4,lwd=2,lty=4)
abline(v=in50,col=1,lty=2); abline(h=0.5,col=1,lty=2)
legend("topleft",c("25_eq5.30","25_eq5.31","30_eq5.31","22_eq5.30N"),
col=c(1,2,3,4),lwd=3,cex=1.1,bty="n",lty=c(1:4))
par(oldp) # return par to old settings; this line not in book
### Dome Shaped Selection
# R-chunk 23 Page 179
#Examples of domed-shaped selectivity curves from domed. Fig.5.15
L <- seq(1,30,1)
p <- c(10,11,16,33,-5,-2)
oldp <- plot1(L,domed(p,L),type="l",lwd=2,ylab="Selectivity",xlab="Age Years")
p1 <- c(8,12,16,33,-5,-1)
lines(L,domed(p1,L),lwd=2,col=2,lty=2)
p2 <- c(9,10,16,33,-5,-4)
lines(L,domed(p2,L),lwd=2,col=4,lty=4)
par(oldp) # return par to old settings; this line not in book
## End(Not run)
[Package MQMF version 0.1.5 Index]