require(here) #package to intelligently figure out project base directory #require(devtools); install_github("eeholmes/VRAP") require(VRAP) knitr::opts_chunk$set(echo = TRUE) workingDirectory <- here() #here() is the base directory of project runSims <- TRUE # set this to true to run all simulations # note that using .at(i,j) for indexing in c++ is safer because it includes bounds checks. ## the following was necessary to get the c++ code to compile Sys.setenv(PATH = paste("C:/Rtools/bin", Sys.getenv("PATH"), sep=";")) Sys.setenv(BINPREF = "C:/Rtools/mingw_$(WIN)/bin/")
In this document we explore alternative implementations of the VRAP model. We reimplement a basic version VRAP first using R and then including C++ code (via Rcpp) to speed things up. VRAP is expanded slightly to allow for auto-correlated log-normal recruitment residuals as well as an alternative to the gamma distribution used in VRAP when there are no co-variates. Model details are described below.
Start with an initial population of fish in the ocean, $N_{premort:~y,a}$, ages 1 through 5. These fish are assumed to have not yet experienced natural and fishing morality for the year
Then, for each age, natural mortality, $mort$, and pre-terminal harvest, $hr$, are applied to generate the post mortality cohort:
$$N_{postmort: y,a} = N_{premort:~y,a}(1-mort_{y,a})(1-hr_{pt:~y,a})$$ A proportion of these mature,
$$N_{mature:~y,a} = N_{postmort:~y,a}mat_{y,a}$$
while the rest remain to form next years pre-mortality cohort. $$N_{premort:~y+1,a+1} = N_{postmort:~y,a}(1-mat_{y,a})$$ The vacant age 1, population is then refilled using the spawner-recruit relationship.
$$N_{premort:~y,1} = R_{y} = SRfunc(E_{y},prod,cap)e^{\beta X_{y}}e^{Z_{y}}$$ based on total adult escapement, $E_{y}$, the spawner-recruit parameters, $prod$ and $cap$, effects of environmental covariates $\beta X_{y}$, and lognormal error with autocorrelation, $Z_{y}$ defined as:
$$Z_{y} \sim Normal(\rho Z_{y-1},\sigma_{rec}\sqrt{1-\rho^2})$$ Note: The current VRAP assumes gamma distributed residuals if there are no covariates.
Total adult escapement is the sum of the mature age 3-5 fish after the terminal harvest.
$$E_y = \sum_{a=3}^5N_{postmort:~y,a}(1-hr_{t:~y,a})$$
Because the model is parameterized in terms of age-specific calendar year harvest rates, as opposed to exploitation rates, there needs to be a process of figuring out what harvest rates to use to achieve the desired exploitation rates. The calendar year exploitation rate is defined here as the total adult equivalent (AEQ) harvest in a year ($harv_{AEQ:~y}$) divided by escapement plus the harvest.
$$ER_y = {harv_{AEQ:~y} \over {harv_{AEQ:~y}+Esc_y}}$$ Here $harv_{AEQ:~y}$ is equal to the total terminal harvest plus the AEQ adjusted non-terminal harvest.
$$harv_{AEQ:~y} = \sum_{a=1}^{5} \left(harv_{PT:~y,a}AEQ_{a} + harv_{T:~y,a}\right)$$ where the preterminal (mixed maturity) harvest is
$$harv_{PT:~y,a} = N_{premort:~y,a}(1-mort_{y,a})hr_{pt:~y,a}$$
and the terminal (mature) harvest is
$$harv_{T:~y,a} = N_{premort:~y,a}(1-mort_{y,a})(1-hr_{pt:~y,a})hr_{t:~y,a}$$
The adult equivalents, $AEQ$, is defined with a recurrence relationship:
$$AEQ_5 = mat_5$$ $$AEQ_{a} = mat_{a} + (1 - mat_{a}) AEQ_{a+1} (1-mort_{a+1})$$
Notice that these $AEQ$ values are applied post natural mortality.
Lets look at some code that scales the harvest rate to achieve the target exploitation rate.
Let's start by defining the initial population and necessary parameteters
# initial calendar year population (ages 1 - 5) CohortStart <- c(8000,4000,2000,700,250) # harvest rates (these will evenutally be scaled) MatU <- c(0,0.05,0.2,0.4,0.1) # pre-termina or mixed-maturity FR PTU <- c(0,0.03, 0.2, 0.3, 0.2) # terminal or mature FR # maturation rates MatRate <- c(0,0.01,0.15,0.7,1) # natural mortality NatMort <- c(0.5,0.4,0.3,0.2,0.1)
We can calculate adult equivalents using the recurrence equation above.
AEQ = c(0,0,0,0,MatRate[5]) for(age in 4:2){ AEQ[age] = MatRate[age] + (1-NatMort[age+1]) * (1 - MatRate[age]) * AEQ[age+1] }
Now let's calculate the exploitation rate based on this information. Here I used code straight from the VRAP CompEscpmnt.R
function.
# initialize cohort Cohort <- CohortStart ### Code straight from the VRAP function CompEscpmnt.R ### # COMPUTE PRETERMINAL MORTALITY AND UPDATE COHORT PTMort = PTU * Cohort TempCohort = Cohort - PTMort # COMPUTE MATURE RUN AND UPDATE COHORT MatRun = TempCohort * MatRate TempCohort = TempCohort - MatRun # COMPUTE MATURE MORTALITY AND ESCAPEMENT MatMort = MatU * MatRun Escpmnt = MatRun - MatMort #spawners Escpmnt[Escpmnt < 1] = 0 # COMPUTE AEQ TOTAL MORTALITY EXPLOITATION RATE AEQMort = AEQ * PTMort + MatMort TotAEQMort = sum(AEQMort) TotEscpmnt = sum(Escpmnt) ### End code ### # exploitation rate ER <- TotAEQMort/(TotAEQMort+TotEscpmnt)
We can simplify this code as follows
# we can simplify the code above to: PTMort <- PTU * Cohort MatMort <- MatU * Cohort * (1 - PTU) * MatRate AEQmort <- PTMort * AEQ + MatMort Escapmnt <- Cohort * (1 - PTU) * (1 - MatU) * MatRate ER <- sum(AEQmort)/(sum(AEQmort)+sum(Escpmnt)) # or even simpler: AEQmort <- Cohort*(PTU*AEQ + (1-PTU)*MatRate*MatU) # AEQMort Escpmnt <- Cohort*(1-PTU)*(1-MatU)*MatRate # Escapmnt ER <- sum(AEQmort)/(sum(AEQmort)+sum(Escpmnt)) # Exploitation rate
Now that we can calculate the exploitation rate we can create an algorithm to adjust the harvest rates until the target exploitation rate is achieved. This is accomplished by starting with a scaling factor of 1, HRscale
, and then multiplying it by the ratio of the target exploitation rate divided by the actual exploitation rate. This process is repeated until the actual exploitation rate is sufficiently close to the target.
targetER <- 0.4 HRscale <- 1 # multiply this times the harvest rates. repeat{ # adjust preterminal and terminal fishing rates PTUAdj <- PTU*HRscale MatUAdj <- MatU*HRscale # calculate AEQ fishing mortality, escapement, and the exploitation rate AEQmort <- Cohort*(PTUAdj*AEQ + (1-PTUAdj)*MatRate*MatUAdj) Escpmnt <- Cohort*(1-PTUAdj)*(1-MatUAdj)*MatRate ER <- sum(AEQmort)/(sum(AEQmort)+sum(Escpmnt)) # calculate the error rate (how far the actual ER is from the target) ERerror <- abs(ER-targetER)/targetER # print the results cat(paste("actual ER = ",round(ER,3),", goal = ",targetER,", abs(actual-target)/target = ",round(ERerror,3),"\n",sep="")) # exit loop if you are close enough if(ERerror < 0.001) break else HRscale <- HRscale*targetER/ER }
Here we create a new VRAP function that takes an input list and returns the simulations. These can then be post-processed with a separate function to generate the desired RER statistics. This version implements a simplified version of VRAP without many of the options. However, it does allow for auto-correlated lognormal recruitment residuals, which is only an option when covariates are included in the original vRAP. This option is made available through the parameter errorType.
First we define an input list that includes all of the necessary information to run the simulation.
inputs <- list( RanSeed = 0, ConvergeCrit = 0.001, SRType = "Ric2", depen = "NO", EscChoice = "YES", # define the productivity and capacity parameters for the spawner-recruit function prod = 2, cap = 1000, # maturation rates MatRate = c(0,0.01,0.15,0.7,1), # natural mortality NatMort = c(0.5,0.4,0.3,0.2,0.1), # harvest rates (these will eventually be scaled) MatU = c(0,0.05,0.2,0.4,0.1), # terminal or mature FR PTU = c(0,0.03, 0.2, 0.3, 0.2), # pre-terminal or mixed-maturity FR # initial popultion CohortStart = c(8000,4000,2000,700,250), # the type of distribution used for the recruitment residuals (gamma or logNormal) errorType = "GAMMA", # If errorType == "GAMMA" SRErrorA and SRerror B are # - the shape and scale parameters of the gamma distribution. # If errorType == "logNormal" they are the # - log of the standard deviation and lag 1 autocorrelation of a normal distribution with mean zero # that is then exponentiated. SRErrorA = 5.7749, SRErrorB = 0.1875, # set management error gamma parameters MgmtError = TRUE, GammaMgmtA = 100, GammaMgmtB = 1/100, # set escapement thresholds ECrit = 200, ERecovery = 400, EndAv = 5, #years to average # parametere for setting exploitation rate start, stop and steps StepFunc = "ER", StepStart = 0, StepEnd = 0.8, StepSize = 0.02, # stepsize in terms of the exploitation rate, baseER*stepSize NYears = 25, # years in each forward projectiong NRuns = 1000 # the number of times the simulation should be repeated. ) inputs$StepNum = round((inputs$StepEnd - inputs$StepStart) / inputs$StepSize + 1) # number of steps
Then we create the function to run the simulations. This allows us to iterate over exploitation rates.
RunSims2R <- function(inputs,rngSeed=NULL,verbose=FALSE){ # initialize output array totEsc <- array(NA,dim=c(inputs$StepNum,inputs$NRuns,inputs$NYears)) HRscale <- 1 # multiply this times the harvest rates. # calcualte the target exploitation rates based on start, step size, and steps targetER <- inputs$StepStart + inputs$StepSize * (0:(inputs$StepNum - 1)) # calculate AEQ and recruitsFromAgeOneFish AEQ <- c(0,0,0,0,inputs$MatRate[5]) for(age in 4:2){ AEQ[age] <- inputs$MatRate[age] + (1-inputs$NatMort[age+1]) * (1 - inputs$MatRate[age]) * AEQ[age+1] } recruitsFromAgeOneFish <- (1-inputs$NatMort[1])*(1-inputs$NatMort[2])*AEQ[2] for(ERind in 1:inputs$StepNum){ print(paste("============= target ER =",targetER[ERind])) for(sim in 1:inputs$NRuns){ # loop through 1000 25 yr simulations logSRerror <- rnorm(1, 0, sd=inputs$SRErrorB) # not currently used Cohort <- inputs$CohortStart # initialize population for(year in 1:25){ # loop through 25 year simulation # apply natural mortality Cohort <- Cohort*(1-inputs$NatMort) # generate management error actualER <- targetER[ERind] if(inputs$MgmtError) actualER <- min(actualER * rgamma(1, inputs$GammaMgmtA, scale=inputs$GammaMgmtB),1) # loop to achieve target exploitation rate unless targetER=0 if(actualER==0){ PTUAdj <- 0 MatUAdj <- 0 }else{ numTrys <- 1 lastAEQmort <- 99 repeat{ # adjust preterminal and terminal fishing rates PTUAdj <- inputs$PTU*HRscale MatUAdj <- inputs$MatU*HRscale # can't be larger than 1 PTUAdj[PTUAdj>1] <- 1 MatUAdj[MatUAdj>1] <- 1 # calculate AEQ fishing mortality, escapement, and the exploitation rate AEQmort <- Cohort*(PTUAdj*AEQ + (1-PTUAdj)*inputs$MatRate*MatUAdj) Escpmnt <- Cohort*(1-PTUAdj)*(1-MatUAdj)*inputs$MatRate totAEQmort <- sum(AEQmort) totEscpmnt <- sum(Escpmnt) ER <- totAEQmort/(totAEQmort+totEscpmnt) # calculate the error rate (how far the actual ER is from the target) ERerror <- abs(ER-actualER)/actualER # exit loop if you are close enough OR other criteria are met. Otherwise adjust HRscale. if(totAEQmort+totEscpmnt < 1 | totAEQmort==0 | numTrys > 100 | totAEQmort==lastAEQmort){ if(verbose){ cat(paste("Target ER = ",targetER[ERind]," Sim = ",sim," Year = ",year, " goal - actual = ",round(actualER,3)," - ",round(ER,3), " HRscale = ",round(HRscale,3)," numTrys = ",numTrys, " totEsc = ",round(totEscpmnt,1)," totAEQmort = ",round(totAEQmort,1),"\n",sep="")) } break }else if(ERerror < 0.001) break else HRscale <- HRscale*actualER/ER numTrys <- numTrys+1 lastAEQmort <- totAEQmort } } # calculate new cohort newCohort <- Cohort*(1-PTUAdj)*(1-inputs$MatRate) Escpmnt <- Cohort*(1-PTUAdj)*(1-MatUAdj)*inputs$MatRate Escpmnt[Escpmnt < 1] <- 0 # calculate adult escapement adultEscapement <- sum(Escpmnt[3:5]) # age the cohort Cohort[2:5] <- newCohort[1:4] # now fill in age 1 fish using the spawner-recruit function. AEQrecruits <- inputs$prod * adultEscapement * exp(-adultEscapement / inputs$cap) if(inputs$errorType=="GAMMA"){ SRerror <- rgamma(1,inputs$SRErrorA,scale=inputs$SRErrorB) }else if(inputs$errorType=="LOGNORMAL"){ # SRErrorA = lognormal sd, SRErrorB = autocorrelation logSRerror <- inputs$SRErrorB*logSRerror + sqrt(1-inputs$SRErrorB^2)*rnorm(1, 0, inputs$SRErrorA) SRerror <- exp(logSRerror) } Cohort[1] <- AEQrecruits*SRerror/recruitsFromAgeOneFish totEsc[ERind,sim,year] <- sum(Escpmnt[2:5]) } } } list(inputs=inputs, totEsc=totEsc) }
Finally we can use this function to generate simulated data.
setwd(workingDirectory) # Only run if var runSims is TRUE, otherwise use saved results. if(runSims){ results <- RunSims2R(inputs) save(results,file="results.Rdat") }else{ load("results.Rdat") }
To compare this to the results from the VRAP package we need to convert the input list into a rav file. This can be accomplished with a simple function.
source(here("R/generateRavFile.R")) print(generateRavFile)
Then we can run the VRAP procedure using the Main function from the VRAP package.
setwd(workingDirectory) generateRavFile(inputs) # Either use the VRAP Main function to run VRAP (twice) and then save the results, or load the saved results. if(runSims){ library(VRAP) vrapOut1 <- Main(InFile="tmp.rav", OutFileBase="vrapOut.tmp", NRuns=1000, silent=TRUE, lcores=4) save(vrapOut1,file="vrapOut1.Rdat") vrapOut2 <- Main(InFile="tmp.rav", OutFileBase="vrapOut.tmp", NRuns=1000, silent=TRUE, lcores=4) save(vrapOut2,file="vrapOut2.Rdat") }else{ load("vrapOut1.Rdat") load("vrapOut2.Rdat") }
Some plots comparing the average escapement values.
# first create the functions plotSimAvg <- function(simDat,vrapOut){ avgs <- apply(simDat[,,],c(1,3),mean) plot(1,1,xlim=c(1,25),ylim=range(c(vrapOut$SummaryStats$AvgEscpmnt,avgs)),xlab="Year", ylab="Average escapement", type="n", bty="l") for(i in 1:inputs$StepNum){ lines(1:25,vrapOut$SummaryStats$AvgEscpmnt[i,]) lines(1:25,avgs[i,],lty=2) } legend(x=15,y=max(avgs),legend=c("VRAP","Simulations above"),lty=c(1,2)) } compPlot<- function(vrapOut1,vrapOut2){ plot(1,1,xlim=c(1,25),ylim=range(c(vrapOut1$SummaryStats$AvgEscpmnt,vrapOut2$SummaryStats$AvgEscpmnt)),xlab="Year", ylab="Average escapement", type="n", bty="l") for(i in 1:inputs$StepNum){ lines(1:25,vrapOut1$SummaryStats$AvgEscpmnt[i,]) lines(1:25,vrapOut2$SummaryStats$AvgEscpmnt[i,],lty=2) } legend(x=15,y=max(vrapOut1$SummaryStats$AvgEscpmnt),legend=c("VRAP 1","VRAP 2"),lty=c(1,2)) } # then use the functions to create the plots plotSimAvg(results$totEsc,vrapOut1) compPlot(vrapOut1,vrapOut2)
Finally, we calculate the proportion of years that escapement is below the lower critical escapement threshold (200 in our example) and compare this to what we get from the VRAP output
bEcrit <- apply(results$totEsc,1,function(x) mean(x<200)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="Year",ylab="% years below critical Esc threshold",type="n",bty="l") targetER <- inputs$StepStart + inputs$StepSize * (0:(inputs$StepNum-1)) lines(targetER,bEcrit) lines(targetER,vrapOut1$SummaryStats$AvgECrit,lty=3) lines(targetER,vrapOut2$SummaryStats$AvgECrit,lty=3) legend(x=0,y=0.8,legend=c("Simulations above","VRAP (2 sims)"),lty=c(1,3))
We can also look a the proportion of years that the geometric mean of the last 5 years of escapement is greater than or equal to the rebuilding escapement threshold (400 in this example).
n <- inputs$NYears meanVals <- apply(results$totEsc,c(1,2),function(x) exp(mean(log(x[(n-4):n])))) aRcrit <- apply(meanVals,1,function(x) mean(x>=400)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="Year",ylab="% years above rebuilding esc threshold",type="n",bty="l") targetER <- inputs$StepStart + inputs$StepSize * (0:(inputs$StepNum-1)) lines(targetER,aRcrit) lines(targetER,vrapOut1$SummaryStats$PropRec,lty=3) lines(targetER,vrapOut2$SummaryStats$PropRec,lty=3) legend(x=0.5,y=0.8,legend=c("Simulations above","VRAP (2 sims)"),lty=c(1,3))
The new simplifed VRAP function appears to produce results that are comparable to the original VRAP code (as defined by the VRAP package). However, the run times are comparable (i.e. not that much faster than the original).
Here, I re-implement most of the simulation code in a C++ function, simFish, using the Rcpp package. The function simFish operates on a single exploitation rate. The code that iterates over different exploitation rates is provided in a separate R function.
#include <Rcpp.h> using namespace Rcpp; // [[Rcpp::export]] NumericMatrix simFish(int NRuns, int NYears, double targetER, bool MgmtError, double GammaMgmtA, double GammaMgmtB, String errorType, double SRErrorA, double SRErrorB, NumericVector CohortStart, double prod, double cap, NumericVector MatRate, NumericVector NatMort, NumericVector PTU, NumericVector MatU, NumericVector AEQ){ NumericMatrix totEsc(NRuns, NYears); NumericVector Cohort(5); NumericVector newCohort(5); NumericVector PTUAdj(5); NumericVector MatUAdj(5); NumericVector AEQmort(5); NumericVector Escpmnt(5); double lastAEQmort, HRscale, totAEQmort, totEscpmnt, adultEscapement, AEQrecruits, ER, realizedER, ERerror, SRerror, logSRerror; bool converged; int numTrys; double recruitsFromAgeOneFish = (1-NatMort(0))*(1-NatMort(1))*AEQ(1); //Rcpp::Rcout << "targetER=" << targetER << "\n"; for(int sim = 0; sim < NRuns; sim++) { // loop through simulations Cohort = clone(CohortStart); // initialize population. Use clone to create a deep copy (i.e don't just copy a pointer) logSRerror = rnorm(1, 0, SRErrorA)[0]; for(int year = 0; year < NYears; year++) { // loop through years Cohort = Cohort*(1-NatMort); if(MgmtError) ER = std::min(targetER * rgamma(1, GammaMgmtA, GammaMgmtB)[0],1.0); else ER = targetER; //Rcpp::Rcout << ER << ","; if(ER==0){ std::fill(PTUAdj.begin(), PTUAdj.end(), 0); // set all elements to 0 std::fill(MatUAdj.begin(), MatUAdj.end(), 0); // set all elements to 0 }else{ numTrys = 1; lastAEQmort = 99; converged = false; HRscale = 1; // Rcpp::Rcout <<"year=" << year << "Cohort=" << Cohort << "\n"; while(!converged){ // adjust preterminal and terminal fishing rates PTUAdj = pmin(PTU*HRscale,1); MatUAdj = pmin(MatU*HRscale,1); // calculate AEQ fishing mortality, escapement, and the exploitation rate AEQmort = Cohort*(PTUAdj*AEQ + (1-PTUAdj)*MatRate*MatUAdj); Escpmnt = Cohort*(1-PTUAdj)*(1-MatUAdj)*MatRate; totAEQmort = sum(AEQmort); totEscpmnt = sum(Escpmnt); realizedER = totAEQmort/(totAEQmort+totEscpmnt); // calculate the error rate (how far the actual ER is from the target) // Rcpp::Rcout << "year=" << year << "ER=" << ER << "\n"; ERerror = std::abs(ER-realizedER)/ER; // exit loop if you are close enough OR other criteria are met //Rcpp::Rcout << "PTUAdj,MatUAdj=" << PTUAdj << "," << MatUAdj << "\n"; //Rcpp::Rcout << "numTrys=" << numTrys << " HRscale=" << HRscale << "\n"; //Rcpp::Rcout << "totAEQmort=" << totAEQmort << " totEscpmnt=" << totEscpmnt << "\n"; //Rcpp::Rcout << "ER=" << ER << " realizedER=" << realizedER << "\n"; //Rcpp::Rcout << "ERerror=" << ERerror << "\n"; if((totAEQmort+totEscpmnt < 1) || (totAEQmort==0) || (numTrys > 100) || (totAEQmort==lastAEQmort)){ converged = true; }else if(ERerror < 0.001){ converged = true; }else{ HRscale = HRscale*ER/realizedER; } numTrys = numTrys+1; lastAEQmort = totAEQmort; } } // calculate new cohort newCohort = Cohort*(1-PTUAdj)*(1-MatRate); // Rcpp::Rcout << "newCohort=" << newCohort << "\n"; Escpmnt = pmax(Cohort*(1-PTUAdj)*(1-MatUAdj)*MatRate,0); // calculate adult escapement adultEscapement = Escpmnt(2) + Escpmnt(3) + Escpmnt(4); // age the cohort for(int ageInd = 0; ageInd < 4; ageInd++) Cohort(ageInd+1) = newCohort(ageInd); // now fill in age 1 fish using the spawner-recruit function. AEQrecruits = prod * adultEscapement * exp(-adultEscapement / cap); if(errorType=="GAMMA"){ SRerror = rgamma(1,SRErrorA,SRErrorB)[0]; }else if(errorType=="LOGNORMAL"){ // SRErrorA = lognormal sd, SRErrorB = autocorrelation logSRerror = SRErrorB*logSRerror + sqrt(1-pow(SRErrorB,2.0))*rnorm(1, 0, SRErrorA)[0]; SRerror = exp(logSRerror); } Cohort(0) = AEQrecruits*SRerror/recruitsFromAgeOneFish; totEsc(sim,year) = Escpmnt(1) + Escpmnt(2) + Escpmnt(3) + Escpmnt(4); } } return totEsc; }
The C++ function can be called from R as follows.
# calculate AEQ AEQ <- c(0,0,0,0,inputs$MatRate[5]) for(age in 4:2){ AEQ[age] <- inputs$MatRate[age] + (1-inputs$NatMort[age+1]) * (1 - inputs$MatRate[age]) * AEQ[age+1] } # run the simulation xx <- simFish(NRuns=1000,NYears=25,targetER=0.2, MgmtError=inputs$MgmtError, GammaMgmtA=inputs$GammaMgmtA,GammaMgmtB=inputs$GammaMgmtB, errorType=inputs$errorType, SRErrorA=inputs$SRErrorA,SRErrorB=inputs$SRErrorB, CohortStart=inputs$CohortStart,prod=inputs$prod,cap=inputs$cap, MatRate=inputs$MatRate,NatMort=inputs$NatMort, PTU=inputs$PTU,MatU=inputs$MatU,AEQ=AEQ) # print out some of the results print(xx[1:10,1:10])
Now we can recreate the full R function that iterates through the different exploitation rates and calls the c++ function.
RunSims2C <- function(inputs,rngSeed=NULL){ # initialize output array totEsc <- array(NA,dim=c(inputs$StepNum, inputs$NRuns, inputs$NYears)) HRscale <- 1 # multiply this times the harvest rates. # calcualte the target exploitation rates based on start, step size, and steps targetER <- inputs$StepStart + inputs$StepSize * (0:(inputs$StepNum-1)) # calculate AEQ and recruitsFromAgeOneFish AEQ <- c(0,0,0,0,inputs$MatRate[5]) for(age in 4:2){ AEQ[age] <- inputs$MatRate[age] + (1-inputs$NatMort[age+1]) * (1 - inputs$MatRate[age]) * AEQ[age+1] } recruitsFromAgeOneFish <- (1-inputs$NatMort[1])*(1-inputs$NatMort[2])*AEQ[2] for(ERind in 1:inputs$StepNum){ print(paste("============= target ER =",targetER[ERind])) totEsc[ERind,,] <- simFish(NRuns=inputs$NRuns,NYears=inputs$NYears,targetER=targetER[ERind], MgmtError=inputs$MgmtError, GammaMgmtA=inputs$GammaMgmtA,GammaMgmtB=inputs$GammaMgmtB, errorType=inputs$errorType, SRErrorA=inputs$SRErrorA,SRErrorB=inputs$SRErrorB, CohortStart=inputs$CohortStart,prod=inputs$prod,cap=inputs$cap, MatRate=inputs$MatRate,NatMort=inputs$NatMort, PTU=inputs$PTU,MatU=inputs$MatU,AEQ=AEQ) } list(inputs=inputs, totEsc=totEsc) }
And then finally we can use this new R function to generate simulations for the range of exploitation rates.
cResults <- RunSims2C(inputs)
Now let's compare the results to the original VRAP (Main from VRAP library) by looking at the average total escapement.
plotSimAvg(cResults$totEsc,vrapOut1)
And the proportion of time below the critical threshold.
bEcrit <- apply(cResults$totEsc,1,function(x) mean(x<200)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="Year",ylab="% years below critical Esc threshold",type="n",bty="l") targetER <- inputs$StepStart + inputs$StepSize * (0:(inputs$StepNum-1)) lines(targetER,bEcrit) lines(targetER,vrapOut1$SummaryStats$AvgECrit,lty=3) lines(targetER,vrapOut2$SummaryStats$AvgECrit,lty=3) legend(x=0,y=0.8,legend=c("Simulations above","VRAP (2 sims)"),lty=c(1,3))
Now that we have comfirmed that both the R and C++ implementations of the VRAP code our producing results comparable to the original VRAP we can look at the difference in speed between the R an C++ implementation.
t1 <- Sys.time() cResults <- RunSims2C(inputs) t2 <- Sys.time() rResults <- RunSims2R(inputs) t3 <- Sys.time()
cat(paste("c-version ",format(t2-t1),"\n",sep="")) cat(paste("r-version ",format(t3-t2),"\n",sep="")) tRat <- round(100*as.numeric(difftime(t2,t1,units="secs"))/as.numeric(difftime(t3,t2,units="secs")),1) cat(paste("The c-version takes ",tRat,"% as much time as the r-version\n",sep="")) cat(paste("Or is ",round(100/tRat)," times faster.\n",sep=""))
In the new implementations of the VRAP (both R and C++) the errorType parameter can be set to "GAMMA" or "LOGNORMAL". This determines the distribution of the recruitment residuals. For the results above we used the "GAMMA" option which is what is assumed when there are no covariates in the original VRAP. In these new implementations, the logNormal distribution allows for autocorrelated logNormal recruitment residuals (as described in the model section above). Here we test this option.
First lets create a new input list with errorType set to "LOGNORMAL"
inputs2 <- inputs inputs2$errorType <- "LOGNORMAL" inputs2$SRErrorA <- 0.5 # lognormal stdev = 0.5 inputs2$SRErrorB <- 0.75 # autocorrelation = 0.75
Now run with the new r and c++ implementations.
resultsLNr <- RunSims2R(inputs2) resultsLNc <- RunSims2C(inputs2)
First lets make sure the R and C++ versions produce similar results.
bEcritR <- apply(resultsLNr$totEsc,1,function(x) mean(x<200)) bEcritC <- apply(resultsLNc$totEsc,1,function(x) mean(x<200)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="Year",ylab="% years below critical Esc threshold",type="n",bty="l") targetER <- inputs2$StepStart + inputs2$StepSize * (0:(inputs2$StepNum-1)) lines(targetER,bEcritR) lines(targetER,bEcritC,lty=3) legend(x=0,y=0.8,legend=c("R function","C++ function"),lty=c(1,3))
Now let's compare simulations with and without autocorrelation.
inputs3 <- inputs2 inputs3$SRErrorB <- 0 resultsLNc2 <- RunSims2C(inputs3)
We can look at a few individual time series of total escapement for runs with and without autocorrelation so see if there is any signs of the autocorrelation in the recruitment residuals translating to patterns in total escapement.
par(mfrow=c(3,1),mar=c(1,1,1,1),oma=c(5,4,5,1)) plot(1:25,resultsLNc$totEsc[1,1,],type="o",pch=16,xlab="",ylab="") plot(1:25,resultsLNc$totEsc[2,1,],type="o",pch=16,xlab="",ylab="") plot(1:25,resultsLNc$totEsc[3,1,],type="o",pch=16,xlab="",ylab="") mtext(side=1,outer=TRUE,text="Year",line=1) mtext(side=2,outer=TRUE,text="Total escapement",line=1) mtext(side=3,outer=TRUE,text="Autocor=0.75",line=1)
par(mfrow=c(3,1),mar=c(1,1,1,1),oma=c(5,4,5,1)) plot(1:25,resultsLNc2$totEsc[1,1,],type="o",pch=16,xlab="",ylab="") plot(1:25,resultsLNc2$totEsc[2,1,],type="o",pch=16,xlab="",ylab="") plot(1:25,resultsLNc2$totEsc[3,1,],type="o",pch=16,xlab="",ylab="") mtext(side=1,outer=TRUE,text="Year",line=1) mtext(side=2,outer=TRUE,text="Total escapement",line=1) mtext(side=3,outer=TRUE,text="Autocor=0",line=1)
And we can compare the proportion of years below the critical threshold for the runs with and without autocorrelation. Notice that the median and sd are the same for the two runs.
bEcritC <- apply(resultsLNc$totEsc,1,function(x) mean(x<200)) bEcritC2 <- apply(resultsLNc2$totEsc,1,function(x) mean(x<200)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="Year",ylab="% years below critical Esc threshold",type="n",bty="l") targetER <- inputs2$StepStart + inputs2$StepSize * (0:(inputs2$StepNum-1)) lines(targetER,bEcritC) lines(targetER,bEcritC2,lty=3) legend(x=0,y=0.8,legend=c("autocor=0.75","autocor=0"),lty=c(1,3))
As might be expected, strong autocorrelation (0.75) leads to a higher proportion of years below the threshold.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.