In eeholmes/VRAPCpp: VRAP 2nd edition with C++ for RER and Viability Computations
require(here) #package to intelligently figure out project base directory
#require(devtools); install_github("eeholmes/VRAP")
require(VRAP)
require(VRAPS)
knitr::opts_chunk$set(echo = TRUE)
vignetteFiles <- here("vignette_files")
if(!file.exists(vignetteFiles)) dir.create(vignetteFiles)
runSims <- TRUE # set this to true to run all simulations
clean <- FALSE # clean up the files made by VRAP 1.0
# 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/")
Introduction
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.
Model
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})$$
ER calibration
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
}
Reimplementation of VRAP in R
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 was only an option when covariates were included in the VRAP 1.0. In VRAP 2.0, this option is made available through the parameter errorType by specifying "LOGNORMAL".
First we define an input list that includes all of the necessary information to run the simulation.
inputs <- list(
Title = "Background Vignette",
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 ER targets
Then we use the code above to create a function to run the simulations for each target ER. This allows us to iterate over a range of exploitation rates. The function is RunSims2R() in the VRAPS package. To view the function type:
library(VRAPS)
RunSims2R
Finally we can use RunSims2R() to generate simulated data.
setwd(vignetteFiles)
# 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")
}
Compared to results from the VRAP package
To compare this to the results from the VRAP package we need to convert the input list into a rav file that VRAP 1.0 can read. This is accomplished the WriteRavFile2() function in the VRAPS package. After loading the VRAPS package with library(VRAPS), you can view the function by typing WriteRavFile2 on the command line.
library(VRAPS)
WriteRavFile
Then we can run the VRAP procedure from VRAP 1.0 using the Main() function from the VRAP package. We will run VRAP twice in order to compare results from two different runs of the simulation. We will save the output to Rdata files to use later for plotting.
setwd(vignetteFiles)
WriteRavFile(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 functions to make the plots
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 (r inputs$ECrit in our example) and compare this to what we get from the VRAP 1.0 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 (r inputs$ERecovery 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))
Faster with Rcpp and C++
The new simplifed VRAP 2.0 function produces results that are comparable to the original VRAP 1.0 code (in the VRAP package). However, the VRAP 2.0 run times are comparable (i.e. not that much faster than VRAP 1.0).
We can re-implement most of the simulation code in a C++ function using the Rcpp package. In VRAP 2.0, this function is called simFish. The function simFish operates on a single exploitation rate. The function code is attached at the end of this vignette.
The C++ function can be called from R as follows. For example's sake, we use a target ER of 0.2.
# 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=inputs$NRuns,NYears=inputs$NYears,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)
We can assemble this code into an R function that iterates through the different exploitation rates and calls the C++ function. This function is RunSims2C() in the VRAPS package.
library(VRAPS)
RunSims2C
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 confirmed that both the R and C++ implementations of the VRAP code are 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=""))
log normal residuals with autocorrelation
In VRAP 2.0 (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 VRAP 1.0. In VRAP 2.0, 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="ER Target",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.
bEcritC <- apply(resultsLNc$totEsc,1,function(x) mean(x<inputs$ECrit))
bEcritC2 <- apply(resultsLNc2$totEsc,1,function(x) mean(x<inputs$ECrit))
plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="ER target",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.
Next steps
- Allow for other SR functions
- Allow for covariates?
- Allow for a vector with NYears values describing the median of recruitment residuals. This could be included as a coraviate with coefficient fixed at 1.
- Consider other residual structures (e.g. auto-correlated sequence + random noise).
if(clean){
file.remove("tmp.rav","tmprav.rav","results.Rdat","vrapOut.tmp.esc","vrapOut.tmp.byr","vrapOut.tmp.sum","vrapOut1.Rdat","vrapOut2.Rdat")
}
simFish C++ 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 do not 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;
}
eeholmes/VRAPCpp documentation built on Feb. 6, 2024, 7:25 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
require(here) #package to intelligently figure out project base directory #require(devtools); install_github("eeholmes/VRAP") require(VRAP) require(VRAPS) knitr::opts_chunk$set(echo = TRUE) vignetteFiles <- here("vignette_files") if(!file.exists(vignetteFiles)) dir.create(vignetteFiles) runSims <- TRUE # set this to true to run all simulations clean <- FALSE # clean up the files made by VRAP 1.0 # 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/")
Introduction
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.
Model
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})$$
ER calibration
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 }
Reimplementation of VRAP in R
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 was only an option when covariates were included in the VRAP 1.0. In VRAP 2.0, this option is made available through the parameter errorType by specifying "LOGNORMAL".
First we define an input list that includes all of the necessary information to run the simulation.
inputs <- list( Title = "Background Vignette", 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 ER targets
Then we use the code above to create a function to run the simulations for each target ER. This allows us to iterate over a range of exploitation rates. The function is RunSims2R() in the VRAPS package. To view the function type:
library(VRAPS)
RunSims2R
Finally we can use RunSims2R() to generate simulated data.
setwd(vignetteFiles) # 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") }
Compared to results from the VRAP package
To compare this to the results from the VRAP package we need to convert the input list into a rav file that VRAP 1.0 can read. This is accomplished the WriteRavFile2() function in the VRAPS package. After loading the VRAPS package with library(VRAPS), you can view the function by typing WriteRavFile2 on the command line.
library(VRAPS)
WriteRavFile
Then we can run the VRAP procedure from VRAP 1.0 using the Main() function from the VRAP package. We will run VRAP twice in order to compare results from two different runs of the simulation. We will save the output to Rdata files to use later for plotting.
setwd(vignetteFiles) WriteRavFile(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 functions to make the plots 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 (r inputs$ECrit in our example) and compare this to what we get from the VRAP 1.0 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 (r inputs$ERecovery 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))
Faster with Rcpp and C++
The new simplifed VRAP 2.0 function produces results that are comparable to the original VRAP 1.0 code (in the VRAP package). However, the VRAP 2.0 run times are comparable (i.e. not that much faster than VRAP 1.0).
We can re-implement most of the simulation code in a C++ function using the Rcpp package. In VRAP 2.0, this function is called simFish. The function simFish operates on a single exploitation rate. The function code is attached at the end of this vignette.
The C++ function can be called from R as follows. For example's sake, we use a target ER of 0.2.
# 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=inputs$NRuns,NYears=inputs$NYears,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)
We can assemble this code into an R function that iterates through the different exploitation rates and calls the C++ function. This function is RunSims2C() in the VRAPS package.
library(VRAPS)
RunSims2C
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 confirmed that both the R and C++ implementations of the VRAP code are 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=""))
log normal residuals with autocorrelation
In VRAP 2.0 (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 VRAP 1.0. In VRAP 2.0, 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="ER Target",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.
bEcritC <- apply(resultsLNc$totEsc,1,function(x) mean(x<inputs$ECrit)) bEcritC2 <- apply(resultsLNc2$totEsc,1,function(x) mean(x<inputs$ECrit)) plot(1,1,xlim=c(0,0.8),ylim=c(0,1),xlab="ER target",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.
Next steps
- Allow for other SR functions
- Allow for covariates?
- Allow for a vector with NYears values describing the median of recruitment residuals. This could be included as a coraviate with coefficient fixed at 1.
- Consider other residual structures (e.g. auto-correlated sequence + random noise).
if(clean){ file.remove("tmp.rav","tmprav.rav","results.Rdat","vrapOut.tmp.esc","vrapOut.tmp.byr","vrapOut.tmp.sum","vrapOut1.Rdat","vrapOut2.Rdat") }
simFish C++ 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 do not 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; }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.