| Title: | Dynamic Model |
|---|---|
| Description: | This package uses comma-delimited data exported from the A&P Excel files to estimate the SR parameters for a user-specified SR function. |
| Authors: | Eli Holmes, Martin Liermann, Howard Coleman, NOAA, Seattle, USA |
| Maintainer: | Elizabeth Holmes - NOAA Federal <[email protected]> |
| License: | GPL-2 |
| Version: | 1.7.5 |
| Built: | 2025-01-27 02:40:40 UTC |
| Source: | https://github.com/nwfsc-math-bio/DM |
Calculate the harvest rate by year and age of return
calcHarvest(bdat)calcHarvest(bdat)
bdat |
data for the bayesian specification of a DM run |
Harvest rate (scalar).
calcInitVals computes initial values for MCMC algorithm used to compute the posterior distribution for the model parameters.
calcInitVals(mlEst, dat, input, bdat, sims = NULL)calcInitVals(mlEst, dat, input, bdat, sims = NULL)
mlEst |
The maximum likelihood estimate for the four SR parameters: log(prod), log(capacity), log(msCoef), log(flowCoef). Output from findOptimum(). |
dat |
data from the A & P file |
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM", SRfunction, includeMarineSurvival, includeFlow |
bdat |
The in a format appropriate for JAGS implementation of the model. |
sims |
MCMC simulations from a previous fit. If sims=NULL all initial values are based on estimates derived from the data. |
A list with initial values for prod (production), logCap (log capacity), msCoef (the coefficient for marine survival), flowCoef (the coeficient for stream survival), rMu and maturationRate_logit.
Calculate harvested recruits and escapement
calculateHarvested(r1)calculateHarvested(r1)
r1 |
a single run list output by the runModel function |
A list with harvested recruits and escapement numbers.
This function takes input and dat produced from the A & P files along with x (the posteriors) from getPosteriors(). It computes the process error estimates (MSE and autoCorr or mean and sterr) along with the SMSY for each posterior draw. It then assembles all these into a dataframe tDat and returns that.
calculatePEandSmsy(input, dat, postdraws)calculatePEandSmsy(input, dat, postdraws)
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM", SRfunction, includeMarineSurvival, includeFlow |
dat |
data from the A & P file |
postdraws |
posterior list from getPostDraws() call |
A dataframe tDat with a column for posterior draws from each parameter and the smsy for each draw.
Center capacity and productivity so that parameters can be compared to estimates from a Bayesian analysis where the covariates were centered.
centerProdCap(input, dat, params)centerProdCap(input, dat, params)
dat |
data from the A & P file |
params |
a vector of parameters with names flowCoef, msCoef, prod, cap. Assumed to be uncentered, e.g. the ML estimates from the SR model with covariates. |
a vector of centered parameters with names flowCoef, msCoef, prod, cap.
Center capacity and productivity so that parameters can be compared to estimates from an analysis where the covariates were (possibly) centered. VRAP uses the uncentered parameters. The input$centerFlow and input$centerMS variables specify whether and which covariates were centered for the Bayesian analysis.
conditionalcenterProdCap(input, dat, params)conditionalcenterProdCap(input, dat, params)
dat |
data from the A & P file |
params |
a vector of parameters with names flowCoef, msCoef, prod, cap. Assumed to be uncentered, e.g. the ML estimates from the SR model with covariates. |
a vector of centered parameters with names flowCoef, msCoef, prod, cap.
This function creates bugs data (bdat).
createBUGSdata(dat, input, priors = NULL, allData = FALSE)createBUGSdata(dat, input, priors = NULL, allData = FALSE)
dat |
data from the A & P file |
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM" or "SS", SRfunction, includeMarineSurvival, includeFlow |
priors |
the parameters for the prior distributions used in the JAGS models. |
allData |
if allData=TRUE the data for both model types is returned. If allData=FALSE, just the data necessary for the specified model is returned. |
dataframe returns a list with the data and priors formatted as input for the JAGS model.
createPriors(pMode = 20, pSig = 7, pMin = 0, pMax = 40, cMu = 9, cSig = 50, cMin = 3, cMax = 12, msMu = 1e-04, msSig = 10, msMin = 1e-04, msMax = 1000, flowMu = -1e-04, flowSig = 10, flowMin = -1000, flowMax = -1e-04)createPriors(pMode = 20, pSig = 7, pMin = 0, pMax = 40, cMu = 9, cSig = 50, cMin = 3, cMax = 12, msMu = 1e-04, msSig = 10, msMin = 1e-04, msMax = 1000, flowMu = -1e-04, flowSig = 10, flowMin = -1000, flowMax = -1e-04)
pMode |
Mode of the production prior (dlnorm) |
pSig |
Sigma of the production prior |
pMin |
Lowerbound of the production prior |
pMax |
Upperbound of the production prior |
cMu |
Mean of the log capacity prior (dnorm) |
cSig |
Sigma of the log capacity prior |
cMin |
Lowerbound of the log capacity prior |
cMax |
Upperbound of the log capacity prior |
msMu |
Mean of the marine survival coef prior (dnorm) |
msSig |
Sigma of the marine survival coef prior |
msMin |
Lowerbound of the marine survival coefficient prior |
msMax |
Upperbound of the marine survival coefficient prior |
flowMu |
Mean of the marine survival coef prior (dnorm) |
flowSig |
Sigma of the marine survival coef prior |
flowMin |
Lowerbound of the flow coefficient prior; the coefficient is always negative |
flowMax |
Upperbound of the flow coefficient prior; the coefficient is always negative |
a list with parameters for each prior.
This function takes the output from the DM model and creates a rav file that can be used forrunning the VRAP simulations. The parameter simInd indicates which MCMC sim to use. This takes the information from the input file in the DM spreadsheet (read in by readDMData) and the results of the MCMC simulations. Current version has no uncertainty in management error and base exploitation rate must be included as a parameter. This could be calculated with the SR params, maturation rates, natMort, and harvest.
createRAVfile(bdat, input, tDat, dat, filename = "temp_rav.rav", estType = "median", sim = 1, rav.options = list())createRAVfile(bdat, input, tDat, dat, filename = "temp_rav.rav", estType = "median", sim = 1, rav.options = list())
bdat |
data for the bayesian specification of a DM run |
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM", SRfunction, includeMarineSurvival, includeFlow |
tDat |
a data frame of the transformed posteriors and calculated process error and Smsy for each draw. From calculatePEandSmsy(). This is specifically for VRAP input. The productivity (a) and capacity (b) posteriors are not adjusted (centered) if the covariates do not have 0 mean. VRAP needs the SR parameters in this raw form using non-centered (de-meaned) covariates. |
dat |
data from the A & P file |
filename |
name to give the rav file |
estType |
("median") |
sim |
indicates which MCMC sim from the posteriors to use if estType is not median |
rav.options |
list of the rav file options to use that do not come from the posteriors.
|
nothing is returned but the rav file is written.
This function takes the input from an A & P files, runs DM to estimate SR parameters and saves a .rav file and .RData file. The .rav file can be downloaded or directly used for running the VRAP simulations by clicking the VRAP tab.
DM()DM()
Because the likelihood surfaces for these problems often has a severe bananna shape with a poorly defined maximum finding an optimum is often non-trivial If an optimum is given in the input file I use that as an initial value for nlm otherwise, I found that nlm often gets stuck in a local min. So here I have iterated between a genetic optimization algorithm and nlm. This works more often but still at times misses the optimum. This is adhoc.
The covariates are log-transformed and centered if input$centerMS=TRUE or input$centerFlow=TRUE.
findOptimum(dat, input, silent = FALSE)findOptimum(dat, input, silent = FALSE)
dat |
data from the A & P file |
input |
a list with the other values needed for a DM run. |
silent |
(TRUE/FALSE) |
In SRFunctions(), bev-holt is defined as S/(S*exp(-p[2])+exp(-p[1]))*exp(p[3]*logMS)*exp(p[4]*logFlow) In DM (writeBUGSmodel.R), R = [S/( (S/exp(logCap)) + (1/prod) )] exp(marineInd*logMS) exp(flowCoef*logFlow) so p[1] = log(prod), constrained to be positive p[2] = log(cap), constrained to be positive p[3] = msCoef, p[4] = flowCoef
A list. $estimate parameters at the minimum sum of squared residuals. The parameters are prod, cap, msCoef, flowCoef. $value is the sum of squared residuals at the minimum.
getPostDraws(dmObj)getPostDraws(dmObj)
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", "bdat" and "tDat" |
a list of vectors of the posterior draws for each parameter.
This uses data in the RData file from a saved DM run to create text describing the model fit.
info(RData.file)info(RData.file)
RData.file |
a saved DM run as a RData file with the objects "input","dat","result", and "tDat" |
a list with text strings.
This function is called by plotResults(). Age composition data versus calendar year. The years with age composition data are represented with a stacked barplot. With age 2 on the bottom (black) and age 5 on the top (white). The numbers above each bar represent the age composition sample size.
plotAgeComp(dmObj, proportions = TRUE)plotAgeComp(dmObj, proportions = TRUE)
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", "bdat" and "tDat" |
proportions |
If proportions=TRUE, the age composition proportions are plotted. If proportions=FALSE, the number of samples in each age are plotted. |
The plot is made.
Capacity posterior distibution. This is the posterior and prior distributions for the spawner-recruit capacity parameter.
plotCapacity(dmObj, xLab = "Capacity parameter", plotType = "histogram", binSize = NULL, xLims = NULL)plotCapacity(dmObj, xLab = "Capacity parameter", plotType = "histogram", binSize = NULL, xLims = NULL)
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", "bdat" and "tDat" |
xLab |
x-axis label |
plotType |
("histogram") |
binSize |
The width of the bin if the plot type is histogram. |
xLims |
x-axis limits |
The plot is made.
plotCoefs(dmObj, xLim = c(-3, 3), xLabs = c("Marine Survival coefficient", "Flow coefficient"))plotCoefs(dmObj, xLim = c(-3, 3), xLabs = c("Marine Survival coefficient", "Flow coefficient"))
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", "bdat" and "tDat" |
xLim |
x-axis limits |
xLabs |
x-axis labels |
The plot is made.
This function is used by plotResults()
plotPriors(priors)plotPriors(priors)
priors |
List of priors parameters. Created with createPriors(). |
The plot is made.
plotProductivity(dmObj, xLab = "Productivity parameter", plotType = "histogram", binSize = NULL, xLims = NULL)plotProductivity(dmObj, xLab = "Productivity parameter", plotType = "histogram", binSize = NULL, xLims = NULL)
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", "bdat" and "tDat" |
xLab |
x-axis label |
plotType |
("histogram") |
binSize |
The width of the bin if the plot type is histogram. |
xLims |
x-axis limits |
The plot is made.
Create several plots based on the DM model runs. This calls the functions for various plots: plotPriors, plotPosteriors, plotAgeComp, plotCoefs
plotResults(dmObj, plotType = c("a.and.p.data", "priors", "SR", "recruits", "covariates", "harvest", "EscData", "EscAndAge", "escapement", "barPlot", "ageComp"), plotDest = "default", plotName = "tmp")plotResults(dmObj, plotType = c("a.and.p.data", "priors", "SR", "recruits", "covariates", "harvest", "EscData", "EscAndAge", "escapement", "barPlot", "ageComp"), plotDest = "default", plotName = "tmp")
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", and "tDat" |
plotType |
The plots and order in which to make them
|
plotDest |
Plot destination. Can be "screen", "png", "default", "none". If "screen", separate graphics windows are opened. If "none", nothing is plotted but the figure captions are returned. |
plotName |
prefix to apply to filename for generated plots if plotDest="png" |
The plot is saved to a file or plotted to screen or default device depending on plotDest. A vector of plot captions is returned invisibly.
The DM tab of the A & P excel file has all the necessary inputs to run DM. This reads that .csv file and makes a list needed for the DM functions. creates input which has the other user specified values for the rav file.
readData(a.and.p.file, input, folder = "./", silent = FALSE)readData(a.and.p.file, input, folder = "./", silent = FALSE)
a.and.p.file |
a csv file saved from the DM tab of an A & P excel file |
input |
a list with the other values needed for a DM run. The following are the defaults used, but can be passed in to specify something different.
|
silent |
Whether to print progress messages. |
The result is a list with dat and input for the DM functions.
This function calls jags to get the posterior draws.
runJAGS(bdat, calcInits, sims = 10000, numChains = 3, details = FALSE)runJAGS(bdat, calcInits, sims = 10000, numChains = 3, details = FALSE)
bdat |
a list of containing the data necessary for the JAGS model |
calcInits |
a list or function that generates a list of intial values for the JAGS model |
sims |
the number of MCMC simulations used for each chain after a burnin of 10 percent of sims |
numChains |
number of MCMC chains to use |
details |
if details=TRUE, more detailed results are returned. |
If details=FALSE, a list of MCMC chains for the monitored model parameters is returned. If details=TRUE, a model object generated by the jags function is returned (this contains the simulations along with additional information)
This function takes the list created by readDMData() and runs the Bayesian model to estimate the SR parameters. The output is written to a .csv file and an .Rdata file. It writes a file, mod1.txt to the user's working directory
runModel(filename, input = NULL, priors = NULL, run = TRUE, sims = 10000, numChains = 3, oldSims = NULL, silent = FALSE)runModel(filename, input = NULL, priors = NULL, run = TRUE, sims = 10000, numChains = 3, oldSims = NULL, silent = FALSE)
filename |
csv file saved from the DM tab of A & P file. Must have age data. |
input |
A list that includes parameters for running the model (see readData for details) of particular interest are the parameters SRfunction, analysisType, and escapementObsSD. SRfunction sets the spawner-recruit function to be used. Current options are "ricker", "bevertonHolt" or "hockeyStick". analysisType can be either "DM" or "SS". DM corresponds to the traditional Dynamic Model which uses the maturation rates from the DynamicsInput tab in the A&P table. SS corresponds to the State Space model which estimates the maturation rates and includes observation error in the spawner counts. If escapementObsSD is a positive real number, observation error is fixed at that value. If it is NULL, observation error is a free parameter that is estimated. |
priors |
parameters for prior distributions in the JAGS model. See |
run |
if run=TRUE, the model is run and the results returned, otherwise all of the information necessary to run the model is returned without running the model |
sims |
length of each MCMC chain; default is 10000 |
numChains |
number of MCMC chains. default is 3 |
oldSims |
? |
silent |
Whether to print progress messages. |
a dmObj: a list of output for a single model fit.
selectSR(x)selectSR(x)
x |
name of SR function as a text string: "ricker", "bevertonHolt", "hockeyStick" |
a function
simplePlot(modelList, modNames = names(modelList)[modOrder], paramList = c("obsSD", "procSD", "prod", "logCap"), funcList = NULL, cols = 4, leftMargin = 12, modOrder = 1:length(modelList))simplePlot(modelList, modNames = names(modelList)[modOrder], paramList = c("obsSD", "procSD", "prod", "logCap"), funcList = NULL, cols = 4, leftMargin = 12, modOrder = 1:length(modelList))
modelList |
a list of multiple run objects (dmObj) based on different models. Each run object in the list is from one call of runModel |
modNames |
more detailed names of models. If the runs are named in the list, it uses those by default. |
paramList |
The parameters you want to compare from the model output. |
funcList |
functions that take the dmObj as their only parameter return a posterior distribution you want to compare across models. the user can make their own functions or utilize one of the available functions (procVarAutoCorr, and SmsyFunc) |
cols |
The number of columns of figures in the plot. |
leftMargin |
the size of the left figure margin (passed to oma in par). |
modOrder |
the order of the models. defaults to the the order in the list. |
The plot is made.
Takes the S, LOG flow, LOG marine survival covariates and p parameters and returns the recruits based on the specified SR function: ricker(), bevertonHolt(), hockeyStick(). Note that it assumes that the log of the covariates are input. All covariates appear in the form exp(p*cov). p[1]=log a; p[2]=log b. The other p are the covariate coefficients.
ricker(S, covariates, p) bevertonHolt(S, covariates, p) hockeyStick(S, covariates, p)ricker(S, covariates, p) bevertonHolt(S, covariates, p) hockeyStick(S, covariates, p)
S |
spawner count |
covariates |
A list with the covariates. The ms covariate is normally the first and the flow covariate is the second. Though order does not matter as the p's are matched based on order the covariate appears in the list. |
p |
the SR function a, b, and the coefficients for each covariate in the covariate list. All covariates appear in the form exp(exp(p)*cov). p[1]=a; p[2]=b. The other p match the covariates. |
recruits (scalar)
recruits from ricker: S*exp(p[1])*exp(-S/exp(p[2]))*exp(p[4]*flow)*exp(p[3]*marineInd)
recruits from beverton Holt: S/(S*exp(-p[2])+exp(-p[1]))*exp(p[4]*flow)*exp(p[3]*marineInd)
recruits from hockey stick: ifelse(exp(p[1])*S<exp(p[2]),exp(p[1])*S,exp(p[2]))*exp(p[4]*flow)*exp(p[3]*marineInd)
Implements the DM model from the original RER process. Estimates escapement by year using escapement data from previous years, maturation rates, harvest rates, and rates of natural mortality. The estimated escapement is log transformed and subtracted from the log transformed values of observed escapement. These differences are then squared and summed to get the sum of squares residual which is returned by the function. This function can be used to estimated the spawner-recruit parameters through the use of an optimization function.
SSTfuncDM(p, SRfunc, dat, input, returnAll = FALSE)SSTfuncDM(p, SRfunc, dat, input, returnAll = FALSE)
p |
The values of the SR function parameters. p[1] and p[2] are logged to keep them positive. p[1]=log(prod), p[2]=log(cap). The other p depend on what covariates are in the model. These are unlogged as they could be positive or negative depending on the covariate.. The length of p is 2 + number of covariates. If there are 2 covariates, then p[3]= msCoef and p[4]=flowCoef. |
SRfunc |
a SR function: ricker, bevertonHolt, hockeyStick, as returned by selectSR() |
dat |
data from the A & P file |
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM", SRfunction, covariates, includeMarineSurvival, includeFlow |
returnAll |
(TRUE/FALSE) |
Returns the sum of squares residuals between the observed and estimated escapement. If returnAll is true, the estimated and observed escapement are also returned.
Summarize DM runs
summarizeRuns(runList)summarizeRuns(runList)
runList |
a list of multiple run objects (dmObj). Each run object in the list is from one call of runModel |
a data frame with summary statistics for each run
Transform capacity and productivity estimates to what they would be if the covariates were centered during estimation (Bayesian or ML). VRAP uses the uncentered parameters. Note, if the coefficients were centered during estimation, input$centerMS=TRUE and/or input$centerFlow=TRUE. These values are used to determine whether to apply the correction. Thus input is assumed to match the parameter values in params.
uncenterProdCap(input, dat, params)uncenterProdCap(input, dat, params)
dat |
data from the A & P file |
params |
a vector of parameters with names flowCoef, msCoef, prod, cap. Assumed to be from the runModel() Bayesian analysis or findOptimum() call, so centered (or not) based on input$centerFlow and input$centerMS |
a vector of uncentered parameters with names flowCoef, msCoef, prod, cap.
centerProdCap
create a mod1.txt which contains the model in the BUGS language
writeBUGScode(input = NULL, outputText = FALSE)writeBUGScode(input = NULL, outputText = FALSE)
input |
a list with the other values needed for a DM run. The following are examples naturalMort, analysisType = "DM", SRfunction, covariates, includeMarineSurvival, includeFlow |
outputText |
(TRUE/FALSE) |
Note that the flow coefficient is always negative and prior set as a normal with upper (negative bound)
nothing but writes the file mod1.txt which is read in to run the BUGS model
writeReport(dmObj = NULL, dmObj.RData.file = NULL, output.file = "report1", rav.options = list(), output.format = c("pdf_document"), output.dir = getwd(), input.type = "xRnw")writeReport(dmObj = NULL, dmObj.RData.file = NULL, output.file = "report1", rav.options = list(), output.format = c("pdf_document"), output.dir = getwd(), input.type = "xRnw")
dmObj |
a saved DM object (list) from runModel() with the objects "input","dat","result", and "tDat" |
dmObj.RData.file |
a saved DM object as a RData file with the objects "input","dat","result", and "tDat" |
output.file |
Filename to give saved report. |
rav.options |
list of the rav file options to use that do not come from the posteriors. |
output.format |
list of output formats. If input.type=="xRmd", any format that rmarkdown:::render allows is fine. If input.type=="xRnw", only "pdf_document" is allowed. |
output.dir |
Directory where to save the output files. Defaults to working directory. |
input.type |
xRmd for rmarkdown. xRnw for knitr. |
Nothing. The report is written to a PDF and tex file.
Write DM output (posteriors) to csv file and RData file.
writeResultsToFile(input, dat, tDat, mlEst, bdat, plist, filename)writeResultsToFile(input, dat, tDat, mlEst, bdat, plist, filename)
input |
a list with the other values needed for a DM run. The following are the defaults used, but can be passed in to specify something different: naturalMort = c(0.5,0.4,0.3,0.2,0.1,0), analysisType = "DM", SRfunction = "ricker", covariates = "no", includeMarineSurvival = "no", includeFlow = "no" |
dat |
data from the A & P file |
tDat |
a data frame of the posteriors and calculated process error and Smsy for each draw. From calculatePEandSmsy(). |
mlEst |
The maximum likelihood estimate for the four SR parameters |
bdat |
data for the bayesian specification of a DM run |
filename |
name to give the outputfiles file |
population |
Name of the population. |
nothing is returned but the csv and RData files are written.