# chapter6 The 53 R-code chunks from On Uncertainty chapter6 is not an active function but rather acts as a repository for the various example code chunks found in chapter6. There are 53 r-code chunks in chapter6. You should, of course, feel free to use and modify any of these example chunks in your own work. ```r chapter6(verbose = TRUE) ``` ## Arguments - `verbose`: Should instructions to written to the console, default=TRUE ## Examples ```r ## Not run: # All the example code from # On Uncertainty # On Uncertainty ## Introduction ### Types of Uncertainty ### The Example Model # R-chunk 1 Page 189 #Fit a surplus production model to abdat fisheries data data(abdat); logce <- log(abdat$cpue) param <- log(c(0.42,9400,3400,0.05)) label=c("r","K","Binit","sigma") # simpspm returns bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat,logobs=logce) outfit(bestmod,title="SP-Model",parnames=label) #backtransforms # R-chunk 2 Page 190 #plot the abdat data and the optimum sp-model fit Fig 6.1 predce <- exp(simpspm(bestmod$estimate,abdat)) optresid <- abdat[,"cpue"]/predce #multiply by predce for obsce ymax <- getmax(c(predce,abdat$cpue)) oldp <- plot1(abdat$year,(predce*optresid),type="l",maxy=ymax,cex=0.9, ylab="CPUE",xlab="Year",lwd=3,col="grey",lty=1) points(abdat$year,abdat$cpue,pch=1,col=1,cex=1.1) lines(abdat$year,predce,lwd=2,col=1) # best fit line par(oldp) # return par to old settings; this line not in book ## Bootstrapping ### Empirical Probability Density Distributions ## A Simple Bootstrap Example # R-chunk 3 Page 193 #regression between catches of NPF prawn species Fig 6.2 data(npf) model <- lm(endeavour ~ tiger,data=npf) oldp <- plot1(npf$tiger,npf$endeavour,type="p",xlab="Tiger Prawn (t)", ylab="Endeavour Prawn (t)",cex=0.9) abline(model,col=1,lwd=2) correl <- sqrt(summary(model)$r.squared) pval <- summary(model)$coefficients[2,4] label <- paste0("Correlation ",round(correl,5)," P = ",round(pval,8)) text(2700,180,label,cex=1.0,font=7,pos=4) par(oldp) # return par to old settings; this line not in book # R-chunk 4 Page 194 # 5000 bootstrap estimates of correlation coefficient Fig 6.3 set.seed(12321) # better to use a less obvious seed, if at all N <- 5000 # number of bootstrap samples result <- numeric(N) #a vector to store 5000 correlations for (i in 1:N) { #sample index from 1:23 with replacement pick <- sample(1:23,23,replace=TRUE) #sample is an R function result[i] <- cor(npf$tiger[pick],npf$endeavour[pick]) } rge <- range(result) # store the range of results CI <- quants(result) # calculate quantiles; 90%CI = 5% and 95% restrim <- result[result > 0] #remove possible -ve values for plot oldp <- parset(cex=1.0) #set up a plot window and draw a histogram bins <- seq(trunc(range(restrim)[1]*10)/10,1.0,0.01) outh <- hist(restrim,breaks=bins,main="",col=0,xlab="Correlation") abline(v=c(correl,mean(result)),col=c(4,3),lwd=c(3,2),lty=c(1,2)) abline(v=CI[c(2,4)],col=4,lwd=2) # and 90% confidence intervals text(0.48,400,makelabel("Range ",rge,sep=" ",sigdig=4),font=7,pos=4) label <- makelabel("90%CI ",CI[c(2,4)],sep=" ",sigdig=4) text(0.48,300,label,cex=1.0,font=7,pos=4) par(oldp) # return par to old settings; this line not in book ## Bootstrapping Time-Series Data # R-chunk 5 Page 196 # fitting Schaefer model with log-normal residuals with 24 years data(abdat); logce <- log(abdat$cpue) # of abalone fisheries data param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) #log values bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat,logobs=logce) optpar <- bestmod$estimate # these are still log-transformed predce <- exp(simpspm(optpar,abdat)) #linear-scale pred cpue optres <- abdat[,"cpue"]/predce # optimum log-normal residual optmsy <- exp(optpar[1])*exp(optpar[2])/4 sampn <- length(optres) # number of residuals and of years # R-chunk 6 Page 196 Table 6.1 code not included in the book outtab <- halftable(cbind(abdat,predce,optres),subdiv=2) kable(outtab, digits=c(0,0,3,3,3,0,0,3,3,3), caption='(ref:tab601)') # R-chunk 7 Pages 196 - 197 # 1000 bootstrap Schaefer model fits; takes a few seconds start <- Sys.time() # use of as.matrix faster than using data.frame bootfish <- as.matrix(abdat) # and avoid altering original data N <- 1000; years <- abdat[,"year"] # need N x years matrices columns <- c("r","K","Binit","sigma") results <- matrix(0,nrow=N,ncol=sampn,dimnames=list(1:N,years)) bootcpue <- matrix(0,nrow=N,ncol=sampn,dimnames=list(1:N,years)) parboot <- matrix(0,nrow=N,ncol=4,dimnames=list(1:N,columns)) for (i in 1:N) { # fit the models and save solutions bootcpue[i,] <- predce * sample(optres, sampn, replace=TRUE) bootfish[,"cpue"] <- bootcpue[i,] #calc and save bootcpue bootmod <- nlm(f=negLL,p=optpar,funk=simpspm,indat=bootfish, logobs=log(bootfish[,"cpue"])) parboot[i,] <- exp(bootmod$estimate) #now save parameters results[i,] <- exp(simpspm(bootmod$estimate,abdat)) #and predce } cat("total time = ",Sys.time()-start, "seconds \n") # R-chunk 8 Page 197 # bootstrap replicates in grey behind main plot Fig 6.4 oldp <- plot1(abdat[,"year"],abdat[,"cpue"],type="n",xlab="Year", ylab="CPUE") # type="n" just lays out an empty plot for (i in 1:N) # ready to add the separate components lines(abdat[,"year"],results[i,],lwd=1,col="grey") points(abdat[,"year"],abdat[,"cpue"],pch=16,cex=1.0,col=1) lines(abdat[,"year"],predce,lwd=2,col=1) par(oldp) # return par to old settings; this line not in book # R-chunk 9 Pages 198 - 199 #histograms of bootstrap parameters and model outputs Fig 6.5 dohist <- function(invect,nmvar,bins=30,bootres,avpar) { #adhoc hist(invect[,nmvar],breaks=bins,main="",xlab=nmvar,col=0) abline(v=c(exp(avpar),bootres[pick,nmvar]),lwd=c(3,2,3,2), col=c(3,4,4,4)) } msy <- parboot[,"r"]*parboot[,"K"]/4 #calculate bootstrap MSY msyB <- quants(msy) #from optimum bootstrap parameters oldp <- parset(plots=c(2,2),cex=0.9) bootres <- apply(parboot,2,quants); pick <- c(2,3,4) #quantiles dohist(parboot,nmvar="r",bootres=bootres,avpar=optpar[1]) dohist(parboot,nmvar="K",bootres=bootres,avpar=optpar[2]) dohist(parboot,nmvar="Binit",bootres=bootres,avpar=optpar[3]) hist(msy,breaks=30,main="",xlab="MSY",col=0) abline(v=c(optmsy,msyB[pick]),lwd=c(3,2,3,2),col=c(3,4,4,4)) par(oldp) # return par to old settings; this line not in book ### Parameter Correlation # R-chunk 10 Page 200 #relationships between parameters and MSY Fig 6.6 parboot1 <- cbind(parboot,msy) # note rgb use, alpha allows for shading, try 1/15 or 1/10 pairs(parboot1,pch=16,col=rgb(red=1,green=0,blue=0,alpha = 1/20)) ## Asymptotic Errors # R-chunk 11 Page 203 #Fit Schaefer model and generate the Hessian data(abdat) param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) # Note inclusion of the option hessian=TRUE in nlm function bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat, logobs=log(abdat[,"cpue"]),hessian=TRUE) outfit(bestmod,backtran = TRUE) #try typing bestmod in console # Now generate the confidence intervals vcov <- solve(bestmod$hessian) # solve inverts matrices sterr <- sqrt(diag(vcov)) #diag extracts diagonal from a matrix optpar <- bestmod$estimate #use qt for t-distrib quantiles U95 <- optpar + qt(0.975,20)*sterr # 4 parameters hence L95 <- optpar - qt(0.975,20)*sterr # (24 - 4) df cat("\n r K Binit sigma \n") cat("Upper 95% ",round(exp(U95),5),"\n") # backtransform cat("Optimum ",round(exp(optpar),5),"\n")#\n =linefeed in cat cat("Lower 95% ",round(exp(L95),5),"\n") ### Uncertainty about the Model Outputs ### Sampling from a Multi-Variate Normal Distribution # R-chunk 12 Page 204 # Use multi-variate normal to generate percentile CI Fig 6.7 library(mvtnorm) # use RStudio, or install.packages("mvtnorm") N <- 1000 # number of multi-variate normal parameter vectors years <- abdat[,"year"]; sampn <- length(years) # 24 years mvncpue <- matrix(0,nrow=N,ncol=sampn,dimnames=list(1:N,years)) columns <- c("r","K","Binit","sigma") # Fill parameter vectors with N vectors from rmvnorm mvnpar <- matrix(exp(rmvnorm(N,mean=optpar,sigma=vcov)), nrow=N,ncol=4,dimnames=list(1:N,columns)) # Calculate N cpue trajectories using simpspm for (i in 1:N) mvncpue[i,] <- exp(simpspm(log(mvnpar[i,]),abdat)) msy <- mvnpar[,"r"]*mvnpar[,"K"]/4 #N MSY estimates # plot data and trajectories from the N parameter vectors oldp <- plot1(abdat[,"year"],abdat[,"cpue"],type="p",xlab="Year", ylab="CPUE",cex=0.9) for (i in 1:N) lines(abdat[,"year"],mvncpue[i,],col="grey",lwd=1) points(abdat[,"year"],abdat[,"cpue"],pch=16,cex=1.0)#orig data lines(abdat[,"year"],exp(simpspm(optpar,abdat)),lwd=2,col=1) par(oldp) # return par to old settings; this line not in book # R-chunk 13 Page 205 #correlations between parameters when using mvtnorm Fig 6.8 pairs(cbind(mvnpar,msy),pch=16,col=rgb(red=1,0,0,alpha = 1/10)) # R-chunk 14 Pages 206 - 207 #N parameter vectors from the multivariate normal Fig 6.9 mvnres <- apply(mvnpar,2,quants) # table of quantiles pick <- c(2,3,4) # select rows for 5%, 50%, and 95% meanmsy <- mean(msy) # optimum bootstrap parameters msymvn <- quants(msy) # msy from mult-variate normal estimates plothist <- function(x,optp,label,resmvn) { hist(x,breaks=30,main="",xlab=label,col=0) abline(v=c(exp(optp),resmvn),lwd=c(3,2,3,2),col=c(3,4,4,4)) } # repeated 4 times, so worthwhile writing a short function oldp <- par(no.readonly=TRUE) par(mfrow=c(2,2),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) plothist(mvnpar[,"r"],optpar[1],"r",mvnres[pick,"r"]) plothist(mvnpar[,"K"],optpar[2],"K",mvnres[pick,"K"]) plothist(mvnpar[,"Binit"],optpar[3],"Binit",mvnres[pick,"Binit"]) plothist(msy,meanmsy,"MSY",msymvn[pick]) par(oldp) # return par to old settings; this line not in book # R-chunk 15 Page 208 Table 6.2 code not included in the book #Tabulate percentile CI from bootstrap (B) and multi-variate (mvn) kable(cbind(bootres,msyB),digits=c(4,3,3,4,3), caption='(ref:tab602)') kable(cbind(mvnres,msymvn),digits=c(4,3,3,4,3)) ## Likelihood Profiles # R-chunk 16 Page 209 #Fit the Schaefer surplus production model to abdat data(abdat); logce <- log(abdat$cpue) # using negLL param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) optmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat,logobs=logce) outfit(optmod,parnames=c("r","K","Binit","sigma")) # R-chunk 17 Page 210 #the code for MQMF's negLLP function negLLP <- function(pars, funk, indat, logobs, initpar=pars, notfixed=c(1:length(pars)),...) { usepar <- initpar #copy the original parameters into usepar usepar[notfixed] <- pars[notfixed] #change 'notfixed' values npar <- length(usepar) logpred <- funk(usepar,indat,...) #funk uses the usepar values pick <- which(is.na(logobs)) # proceed as in negLL if (length(pick) > 0) { LL <- -sum(dnorm(logobs[-pick],logpred[-pick],exp(pars[npar]), log=T)) } else { LL <- -sum(dnorm(logobs,logpred,exp(pars[npar]),log=T)) } return(LL) } # end of negLLP # R-chunk 18 Page 211 #does negLLP give same answers as negLL when no parameters fixed? param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) bestmod <- nlm(f=negLLP,p=param,funk=simpspm,indat=abdat,logobs=logce) outfit(bestmod,parnames=c("r","K","Binit","sigma")) # R-chunk 19 Page 211 #Likelihood profile for r values 0.325 to 0.45 rval <- seq(0.325,0.45,0.001) # set up the test sequence ntrial <- length(rval) # create storage for the results columns <- c("r","K","Binit","sigma","-veLL") result <- matrix(0,nrow=ntrial,ncol=length(columns), dimnames=list(rval,columns))# close to optimum bestest <- c(r= 0.32,K=11000,Binit=4000,sigma=0.05) for (i in 1:ntrial) { #i <- 1 param <- log(c(rval[i],bestest[2:4])) #recycle bestest values parinit <- param #to improve the stability of nlm as r changes bestmodP <- nlm(f=negLLP,p=param,funk=simpspm,initpar=parinit, indat=abdat,logobs=log(abdat$cpue),notfixed=c(2:4), typsize=magnitude(param),iterlim=1000) bestest <- exp(bestmodP$estimate) result[i,] <- c(bestest,bestmodP$minimum) # store each result } minLL <- min(result[,"-veLL"]) #minimum across r values used. # R-chunk 20 Page 212 Table 6.3 code not included in the book #tabulate first 12 records from likelihood profile kable(head(result,12),digits=c(3,3,3,4,5), caption='(ref:tab603)') ### Likelihood Ratio Based Confidence Intervals # R-chunk 21 Page 213 #likelihood profile on r from the Schaefer model Fig 6.10 plotprofile(result,var="r",lwd=2) # review the code # R-chunk 22 Page 214 #Likelihood profile for K values 7200 to 12000 Kval <- seq(7200,12000,10) ntrial <- length(Kval) columns <- c("r","K","Binit","sigma","-veLL") resultK <- matrix(0,nrow=ntrial,ncol=length(columns), dimnames=list(Kval,columns)) bestest <- c(r= 0.45,K=7500,Binit=2800,sigma=0.05) for (i in 1:ntrial) { param <- log(c(bestest[1],Kval[i],bestest[c(3,4)])) parinit <- param bestmodP <- nlm(f=negLLP,p=param,funk=simpspm,initpar=parinit, indat=abdat,logobs=log(abdat$cpue), notfixed=c(1,3,4),iterlim=1000) bestest <- exp(bestmodP$estimate) resultK[i,] <- c(bestest,bestmodP$minimum) } minLLK <- min(resultK[,"-veLL"]) #kable(head(result,12),digits=c(4,3,3,4,5)) # if wanted. # R-chunk 23 Page 214 #likelihood profile on K from the Schaefer model Fig 6.11 plotprofile(resultK,var="K",lwd=2) ### -ve Log-Likelihoods or Likelihoods # R-chunk 24 Page 215 #translate -velog-likelihoods into likelihoods likes <- exp(-resultK[,"-veLL"])/sum(exp(-resultK[,"-veLL"]),na.rm=TRUE) resK <- cbind(resultK,likes,cumlike=cumsum(likes)) # R-chunk 25 Page 216 Table 6.4 code not included in the book #tabulate head of likelihood profile matrix for K kable(head(resK,8),digits=c(4,0,3,4,5,9,7),caption='(ref:tab604)') # R-chunk 26 Page 216 Figure 6.12 code not in the book #K parameter likelihood profile Fig 6.12 oldp <- plot1(resK[,"K"],resK[,"likes"],xlab="K value", ylab="Likelihood",lwd=2) lower <- which.closest(0.025,resK[,"cumlike"]) mid <- which(resK[,"likes"] == max(resK[,"likes"])) upper <- which.closest(0.975,resK[,"cumlike"]) abline(v=c(resK[c(lower,mid,upper),"K"]),col=1,lwd=c(1,2,1)) label <- makelabel("",resK[c(lower,mid,upper),"K"],sep=" ") text(9500,0.005,label,cex=1.2,pos=4) par(oldp) # return par to old settings; this line not in book ### Percentile Likelihood Profiles for Model Outputs # R-chunk 27 Page 217 - 218 #examine effect on -veLL of MSY values from 740 - 1050t #need a different negLLP() function, negLLO(): O for output. #now optvar=888.831 (rK/4), the optimum MSY, varval ranges 740-1050 #and wght is the weighting to give to the penalty negLLO <- function(pars,funk,indat,logobs,wght,optvar,varval) { logpred <- funk(pars,indat) LL <- -sum(dnorm(logobs,logpred,exp(tail(pars,1)),log=T)) + wght*((varval - optvar)/optvar)^2 #compare with negLL return(LL) } # end of negLLO msyP <- seq(740,1020,2.5); optmsy <- exp(optmod$estimate[1])*exp(optmod$estimate[2])/4 ntrial <- length(msyP) wait <- 400 columns <- c("r","K","Binit","sigma","-veLL","MSY","pen","TrialMSY") resultO <- matrix(0,nrow=ntrial,ncol=length(columns),dimnames=list(msyP,columns)) bestest <- c(r= 0.47,K=7300,Binit=2700,sigma=0.05) for (i in 1:ntrial) { # i <- 1 param <- log(bestest) bestmodO <- nlm(f=negLLO,p=param,funk=simpspm,indat=abdat, logobs=log(abdat$cpue),wght=wait, optvar=optmsy,varval=msyP[i],iterlim=1000) bestest <- exp(bestmodO$estimate) ans <- c(bestest,bestmodO$minimum,bestest[1]*bestest[2]/4, wait *((msyP[i] - optmsy)/optmsy)^2,msyP[i]) resultO[i,] <- ans } minLLO <- min(resultO[,"-veLL"]) # R-chunk 28 Page 218 Table 6.5 code not included in the book #tabulate first and last few records of profile on MSY kable(head(resultO[,1:7],4),digits=c(3,3,3,4,2,3,2),caption='(ref:tab605)') kable(tail(resultO[,1:7],4),digits=c(3,3,3,4,2,3,2)) # R-chunk 29 Page 219 Figure 6.13 code not included in the book #likelihood profile on MSY from the Schaefer model Fig 6.13 plotprofile(resultO,var="TrialMSY",lwd=2) ## Bayesian Posterior Distributions ### Generating the Markov Chain ### The Starting Point ### The Burn-in Period ### Convergence to the Stationary Distribution ### The Jumping Distribution ### Application of MCMC to the Example # R-chunk 30 Page 225 #activate and plot the fisheries data in abdat Fig 6.14 data(abdat) # type abdat in the console to see contents plotspmdat(abdat) #use helper function to plot fishery stats vs year ### Markov Chain Monte Carlo ### A First Example of an MCMC # R-chunk 31 Pages 228 - 229 # Conduct MCMC analysis to illustrate burn-in. Fig 6.15 data(abdat); logce <- log(abdat$cpue) fish <- as.matrix(abdat) # faster to use a matrix than a data.frame! begin <- Sys.time() # enable time taken to be calculated chains <- 1 # 1 chain per run; normally do more burnin <- 0 # no burn-in for first three chains N <- 100 # Number of MCMC steps to keep step <- 4 # equals one step per parameter so no thinning priorcalc <- calcprior # define the prior probability function scales <- c(0.065,0.055,0.065,0.425) #found by trial and error set.seed(128900) #gives repeatable results in book; usually omitted inpar <- log(c(r= 0.4,K=11000,Binit=3600,sigma=0.05)) result1 <- do_MCMC(chains,burnin,N,step,inpar,negLL,calcpred=simpspm, calcdat=fish,obsdat=logce,priorcalc,scales) inpar <- log(c(r= 0.35,K=8500,Binit=3400,sigma=0.05)) result2 <- do_MCMC(chains,burnin,N,step,inpar,negLL,calcpred=simpspm, calcdat=fish,obsdat=logce,priorcalc,scales) inpar <- log(c(r= 0.45,K=9500,Binit=3200,sigma=0.05)) result3 <- do_MCMC(chains,burnin,N,step,inpar,negLL,calcpred=simpspm, calcdat=fish,obsdat=logce,priorcalc,scales) burnin <- 50 # strictly a low thinning rate of 4; not enough step <- 16 # 16 thinstep rate = 4 parameters x 4 = 16 N <- 10000 # 16 x 10000 = 160,000 steps + 50 burnin inpar <- log(c(r= 0.4,K=9400,Binit=3400,sigma=0.05)) result4 <- do_MCMC(chains,burnin,N,step,inpar,negLL,calcpred=simpspm, calcdat=fish,obsdat=logce,priorcalc,scales) post1 <- result1[[1]][[1]] post2 <- result2[[1]][[1]] post3 <- result3[[1]][[1]] postY <- result4[[1]][[1]] cat("time = ",Sys.time() - begin,"\n") cat("Accept = ",result4[[2]],"\n") # R-chunk 32 Pages 229 - 230 #first example and start of 3 initial chains for MCMC Fig6.15 oldp <- parset(cex=0.85) P <- 75 # the first 75 steps only start to explore parameter space plot(postY[,"K"],postY[,"r"],type="p",cex=0.2,xlim=c(7000,13000), ylim=c(0.28,0.47),col=8,xlab="K",ylab="r",panel.first=grid()) lines(post2[1:P,"K"],post2[1:P,"r"],lwd=1,col=1) points(post2[1:P,"K"],post2[1:P,"r"],pch=15,cex=1.0) lines(post1[1:P,"K"],post1[1:P,"r"],lwd=1,col=1) points(post1[1:P,"K"],post1[1:P,"r"],pch=1,cex=1.2,col=1) lines(post3[1:P,"K"],post3[1:P,"r"],lwd=1,col=1) points(post3[1:P,"K"],post3[1:P,"r"],pch=2,cex=1.2,col=1) par(oldp) # return par to old settings; this line not in book # R-chunk 33 Pages 230 - 231 #pairs plot of parameters from the first MCMC Fig 6.16 posterior <- result4[[1]][[1]] msy <-posterior[,1]*posterior[,2]/4 pairs(cbind(posterior[,1:4],msy),pch=16,col=rgb(1,0,0,1/50),font=7) # R-chunk 34 Pages 231 - 232 #plot the traces from the first MCMC example Fig 6.17 posterior <- result4[[1]][[1]] oldp <- par(no.readonly=TRUE) # this line not in book par(mfrow=c(4,2),mai=c(0.4,0.4,0.05,0.05),oma=c(0.0,0,0.0,0.0)) par(cex=0.8, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) label <- colnames(posterior) N <- dim(posterior)[1] for (i in 1:4) { ymax <- getmax(posterior[,i]); ymin <- getmin(posterior[,i]) plot(1:N,posterior[,i],type="l",lwd=1,ylim=c(ymin,ymax), panel.first=grid(),ylab=label[i],xlab="Step") plot(density(posterior[,i]),lwd=2,col=2,panel.first=grid(),main="") } par(oldp) # return par to old settings; this line not in book # R-chunk 35 Page 233 #Use acf to examine auto-correlation with thinstep = 16 Fig 6.18 posterior <- result4[[1]][[1]] label <- colnames(posterior)[1:4] oldp <- parset(plots=c(2,2),cex=0.85) for (i in 1:4) auto <- acf(posterior[,i],type="correlation",lwd=2, plot=TRUE,ylab=label[i],lag.max=20) par(oldp) # return par to old settings; this line not in book # R-chunk 36 Pages 233 - 234 #setup MCMC with thinstep of 128 per parameter Fig 6.19 begin=gettime() scales <- c(0.06,0.05,0.06,0.4) inpar <- log(c(r= 0.4,K=9400,Binit=3400,sigma=0.05)) result <- do_MCMC(chains=1,burnin=100,N=1000,thinstep=512,inpar, negLL,calcpred=simpspm,calcdat=fish, obsdat=logce,calcprior,scales,schaefer=TRUE) posterior <- result[[1]][[1]] label <- colnames(posterior)[1:4] oldp <- parset(plots=c(2,2),cex=0.85) for (i in 1:4) auto <- acf(posterior[,i],type="correlation",lwd=2, plot=TRUE,ylab=label[i],lag.max=20) par(oldp) # return par to old settings; this line not in book cat(gettime() - begin) ### Marginal Distributions # R-chunk 37 Pages 235 - 236 # plot marginal distributions from the MCMC Fig 6.20 dohist <- function(x,xlab) { # to save a little space return(hist(x,main="",breaks=50,col=0,xlab=xlab,ylab="", panel.first=grid())) } # ensure we have the optimum solution available param <- log(c(r= 0.42,K=9400,Binit=3400,sigma=0.05)) bestmod <- nlm(f=negLL,p=param,funk=simpspm,indat=abdat, logobs=log(abdat$cpue)) optval <- exp(bestmod$estimate) posterior <- result[[1]][[1]] #example above N=1000, thin=512 oldp <- par(no.readonly=TRUE) par(mfrow=c(5,1),mai=c(0.4,0.3,0.025,0.05),oma=c(0,1,0,0)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) np <- length(param) for (i in 1:np) { #store invisible output from hist for later use outH <- dohist(posterior[,i],xlab=colnames(posterior)[i]) abline(v=optval[i],lwd=3,col=4) tmp <- density(posterior[,i]) scaler <- sum(outH$counts)*(outH$mids[2]-outH$mids[1]) tmp$y <- tmp$y * scaler lines(tmp,lwd=2,col=2) } msy <- posterior[,"r"]*posterior[,"K"]/4 mout <- dohist(msy,xlab="MSY") tmp <- density(msy) tmp$y <- tmp$y * (sum(mout$counts)*(mout$mids[2]-mout$mids[1])) lines(tmp,lwd=2,col=2) abline(v=(optval[1]*optval[2]/4),lwd=3,col=4) mtext("Frequency",side=2,outer=T,line=0.0,font=7,cex=1.0) par(oldp) # return par to old settings; this line not in book ## The Use of Rcpp # R-chunk 38 Pages 236 - 237 #profile the running of do_MCMC using the now well known abdat data(abdat); logce <- log(abdat$cpue); fish <- as.matrix(abdat) param <- log(c(r=0.39,K=9200,Binit=3400,sigma=0.05)) Rprof(append=TRUE) # note the use of negLL1() result <- do_MCMC(chains=1,burnin=100,N=20000,thinstep=16,inpar=param, infunk=negLL1,calcpred=simpspm,calcdat=fish, obsdat=logce,priorcalc=calcprior, scales=c(0.07,0.06,0.07,0.45)) Rprof(NULL) outprof <- summaryRprof() # R-chunk 39 Page 238 Table 6.6 code not included in the book #tabulate output of Rprof on do_MCMC function kable(head(outprof$by.self,12),caption='(ref:tab606)') ### Addressing Vectors and Matrices ### Replacement for simpspm() # R-chunk 40 Page 240 library(Rcpp) #Send a text string containing the C++ code to cppFunction this will #take a few seconds to compile, then the function simpspmC will #continue to be available during the rest of your R session. The #code in this chunk could be included into its own R file, and then #the R source() function can be used to include the C++ into a #session. indat must have catch in col2 (col1 in C++), and cpue in #col3 (col2 in C++). Note the use of ; at the end of each line. #Like simpspm(), this returns only the log(predicted cpue). 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; }') # R-chunk 41 Page 241 #Ensure results obtained from simpspm and simpspmC are same library(microbenchmark) data(abdat) fishC <- as.matrix(abdat) # Use a matrix rather than a data.frame inpar <- log(c(r= 0.389,K=9200,Binit=3300,sigma=0.05)) spmR <- exp(simpspm(inpar,fishC)) # demonstrate equivalence #need to declare all arguments in simpspmC, no default values spmC <- exp(simpspmC(inpar,fishC,schaefer=TRUE)) out <- microbenchmark( # verything identical calling function simpspm(inpar,fishC,schaefer=TRUE), simpspmC(inpar,fishC,schaefer=TRUE), times=1000 ) out2 <- summary(out)[,2:8] out2 <- rbind(out2,out2[2,]/out2[1,]) rownames(out2) <- c("simpspm","simpspmC","TimeRatio") # R-chunk 42 Page 241 Table 6.7 code not included in the book #compare results from simpspm and simpspmC kable(halftable(cbind(spmR,spmC)),row.names=TRUE,digits=c(4,4,4,4,4,4),caption='(ref:tab607)') # R-chunk 43 Page 242 Table 6.8 code not included in the book #output from microbenchmark comparison of simpspm and simpspmC kable(out2,row.names=TRUE,digits=c(3,3,3,3,3,3,3,0),caption='(ref:tab608)') # R-chunk 44 Pages 242 - 243 #How much does using simpspmC in do_MCMC speed the run time? #Assumes Rcpp code has run, eg source("Rcpp_functions.R") set.seed(167423) #Can use getseed() to generate a suitable seed beginR <- gettime() #to enable estimate of time taken setscale <- c(0.07,0.06,0.07,0.45) reps <- 2000 #Not enough but sufficient for demonstration param <- log(c(r=0.39,K=9200,Binit=3400,sigma=0.05)) resultR <- do_MCMC(chains=1,burnin=100,N=reps,thinstep=128, inpar=param,infunk=negLL1,calcpred=simpspm, calcdat=fishC,obsdat=log(abdat$cpue),schaefer=TRUE, priorcalc=calcprior,scales=setscale) timeR <- gettime() - beginR cat("time = ",timeR,"\n") cat("acceptance rate = ",resultR$arate," \n") postR <- resultR[[1]][[1]] set.seed(167423) # Use the same pseudo-random numbers and the beginC <- gettime() # same starting point to make the comparsion param <- log(c(r=0.39,K=9200,Binit=3400,sigma=0.05)) resultC <- do_MCMC(chains=1,burnin=100,N=reps,thinstep=128, inpar=param,infunk=negLL1,calcpred=simpspmC, calcdat=fishC,obsdat=log(abdat$cpue),schaefer=TRUE, priorcalc=calcprior,scales=setscale) timeC <- gettime() - beginC cat("time = ",timeC,"\n") # note the same acceptance rates cat("acceptance rate = ",resultC$arate," \n") postC <- resultC[[1]][[1]] cat("Time Ratio = ",timeC/timeR) # R-chunk 45 Page 243 #compare marginal distributions of the 2 chains Fig 6.21 oldp <- par(no.readonly=TRUE) # this line not in the book par(mfrow=c(1,1),mai=c(0.45,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) maxy <- getmax(c(density(postR[,"K"])$y,density(postC[,"K"])$y)) plot(density(postR[,"K"]),lwd=2,col=1,xlab="K",ylab="Density", main="",ylim=c(0,maxy),panel.first=grid()) lines(density(postC[,"K"]),lwd=3,col=5,lty=2) par(oldp) # return par to old settings; this line not in book ### Multiple Independent Chains # R-chunk 46 Page 244 #run multiple = 3 chains setscale <- c(0.07,0.06,0.07,0.45) # I only use a seed for set.seed(9393074) # reproducibility within this book reps <- 10000 # reset the timer beginC <- gettime() # remember a thinstep=256 is insufficient resultC <- do_MCMC(chains=3,burnin=100,N=reps,thinstep=256, inpar=param,infunk=negLL1,calcpred=simpspmC, calcdat=fishC,obsdat=log(fishC[,"cpue"]), priorcalc=calcprior,scales=setscale,schaefer=TRUE) cat("time = ",gettime() - beginC," secs \n") # R-chunk 47 Pages 244 - 245 #3 chain run using simpspmC, 10000 reps, thinstep=256 Fig 6.22 oldp <- par(no.readonly=TRUE) par(mfrow=c(2,2),mai=c(0.4,0.45,0.05,0.05),oma=c(0.0,0,0.0,0.0)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) label <- c("r","K","Binit","sigma") for (i in 1:4) { plot(density(resultC$result[[2]][,i]),lwd=2,col=1, xlab=label[i],ylab="Density",main="",panel.first=grid()) lines(density(resultC$result[[1]][,i]),lwd=2,col=2) lines(density(resultC$result[[3]][,i]),lwd=2,col=3) } par(oldp) # return par to old settings; this line not in book # R-chunk 48 Pages 245 - 246 #generate summary stats from the 3 MCMC chains av <- matrix(0,nrow=3,ncol=4,dimnames=list(1:3,label)) sig2 <- av # do the variance relsig <- av # relative to mean of all chains for (i in 1:3) { tmp <- resultC$result[[i]] av[i,] <- apply(tmp[,1:4],2,mean) sig2[i,] <- apply(tmp[,1:4],2,var) } cat("Average \n") av cat("\nVariance per chain \n") sig2 cat("\n") for (i in 1:4) relsig[,i] <- sig2[,i]/mean(sig2[,i]) cat("Variance Relative to Mean Variance of Chains \n") relsig # R-chunk 49 Pages 246 - 247 #compare quantile from the 2 most widely separate MCMC chains tmp <- resultC$result[[2]] # the 10000 values of each parameter cat("Chain 2 \n") msy1 <- tmp[,"r"]*tmp[,"K"]/4 ch1 <- apply(cbind(tmp[,1:4],msy1),2,quants) round(ch1,4) tmp <- resultC$result[[3]] cat("Chain 3 \n") msy2 <- tmp[,"r"]*tmp[,"K"]/4 ch2 <- apply(cbind(tmp[,1:4],msy2),2,quants) round(ch2,4) cat("Percent difference ") cat("\n2.5% ",round(100*(ch1[1,] - ch2[1,])/ch1[1,],4),"\n") cat("50% ",round(100*(ch1[3,] - ch2[3,])/ch1[3,],4),"\n") cat("97.5% ",round(100*(ch1[5,] - ch2[5,])/ch1[5,],4),"\n") ### Replicates Required to Avoid Serial Correlation # R-chunk 50 Page 248 #compare two higher thinning rates per parameter in MCMC param <- log(c(r=0.39,K=9200,Binit=3400,sigma=0.05)) setscale <- c(0.07,0.06,0.07,0.45) result1 <- do_MCMC(chains=1,burnin=100,N=2000,thinstep=1024, inpar=param,infunk=negLL1,calcpred=simpspmC, calcdat=fishC,obsdat=log(abdat$cpue), priorcalc=calcprior,scales=setscale,schaefer=TRUE) result2 <- do_MCMC(chains=1,burnin=50,N=1000,thinstep=2048, inpar=param,infunk=negLL1,calcpred=simpspmC, calcdat=fishC,obsdat=log(abdat$cpue), priorcalc=calcprior,scales=setscale,schaefer=TRUE) # R-chunk 51 Page 248 #autocorrelation of 2 different thinning rate chains Fig6.23 posterior1 <- result1$result[[1]] posterior2 <- result2$result[[1]] label <- colnames(posterior1)[1:4] oldp <- par(no.readonly=TRUE) par(mfrow=c(4,2),mai=c(0.25,0.45,0.05,0.05),oma=c(1.0,0,1.0,0.0)) par(cex=0.85, mgp=c(1.35,0.35,0), font.axis=7,font=7,font.lab=7) for (i in 1:4) { auto <- acf(posterior1[,i],type="correlation",plot=TRUE, ylab=label[i],lag.max=20,xlab="",ylim=c(0,0.3),lwd=2) if (i == 1) mtext(1024,side=3,line=-0.1,outer=FALSE,cex=1.2) auto <- acf(posterior2[,i],type="correlation",plot=TRUE, ylab=label[i],lag.max=20,xlab="",ylim=c(0,0.3),lwd=2) if (i == 1) mtext(2048,side=3,line=-0.1,outer=FALSE,cex=1.2) } mtext("Lag",side=1,line=-0.1,outer=TRUE,cex=1.2) par(oldp) # return par to old settings; this line not in book # R-chunk 52 Page 249 #visual comparison of 2 chains marginal densities Fig 6.24 oldp <- parset(plots=c(2,2),cex=0.85) label <- c("r","K","Binit","sigma") for (i in 1:4) { plot(density(result1$result[[1]][,i]),lwd=4,col=1,xlab=label[i], ylab="Density",main="",panel.first=grid()) lines(density(result2$result[[1]][,i]),lwd=2,col=5,lty=2) } par(oldp) # return par to old settings; this line not in book # R-chunk 53 Pages 250 - 251 #tablulate a summary of the two different thinning rates. cat("1024 thinning rate \n") posterior <- result1$result[[1]] msy <-posterior[,1]*posterior[,2]/4 tmp1 <- apply(cbind(posterior[,1:4],msy),2,quants) rge <- apply(cbind(posterior[,1:4],msy),2,range) tmp1 <- rbind(tmp1,rge[2,] - rge[1,]) rownames(tmp1)[6] <- "Range" print(round(tmp1,4)) posterior2 <- result2$result[[1]] msy2 <-posterior2[,1]*posterior2[,2]/4 cat("2048 thinning rate \n") tmp2 <- apply(cbind(posterior2[,1:4],msy2),2,quants) rge2 <- apply(cbind(posterior2[,1:4],msy2),2,range) tmp2 <- rbind(tmp2,rge2[2,] - rge2[1,]) rownames(tmp2)[6] <- "Range" print(round(tmp2,4)) cat("Inner 95% ranges and Differences between total ranges \n") cat("95% 1 ",round((tmp1[5,] - tmp1[1,]),4),"\n") cat("95% 2 ",round((tmp2[5,] - tmp2[1,]),4),"\n") cat("Diff ",round((tmp2[6,] - tmp1[6,]),4),"\n") ## End(Not run) ```

chapter6 function