# chapter5 The 23 R-code chunks from Static Models 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. ```r chapter5(verbose = TRUE) ``` ## Arguments - `verbose`: Should instructions to written to the console, default=TRUE ## Examples ```r ## 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) ```

chapter5 function