Change params in dynamics() into a list (now matches projections())

Author: mcsiple

R/20_dynamics.R

@@ -5,31 +5,41 @@
 #' @details
 #' The population model is a single-sex age-structured model in which the number of calves or pups born each year is density dependent, with the extent of density dependence a function of the number of mature adults \eqn{\tildeN}, the fecundity (pregnancy rate) at pre-exploitation equilibrium \eqn{f_0}, the maximum theoretical fecundity rate fmax, the degree of compensation \eqn{z}, and the abundance of individuals aged 1+ \eqn{N_{y+1}^{1+}} relative to carrying capacity \eqn{K^{1+}}. This function can be used alone but is intended to be used with \code{Projections()} to generate multiple simulations. NOTE: Either \code{ConstantCatch} or \code{ConstantF} can be specified, but not both.
 #'
-#' @param S0 Calf/pup survival, a numeric value between 0 and 1
-#' @param S1plus Survival for animals age 1 year and older, a numeric value between 0 and 1
-#' @param K1plus The pre-exploitation population size of individuals aged 1 and older.  If this value is unavailable, it can be approximated by using the initial depletion and the estimate of current abundance
-#' @param AgeMat Age at maturity in years (assumed to be age at first parturition - 1).
+#' @param lh.params A list containing life history parameters: \cr\cr
+#' \code{S0} Calf/pup survival, a numeric value between 0 and 1 \cr\cr
+#' \code{S1plus} Survival for animals age 1 year and older, a numeric value between 0 and 1 \cr\cr
+#' \code{K1plus} The pre-exploitation population size of individuals aged 1 and older.  If this value is unavailable, it can be approximated by using the initial depletion and the estimate of current abundance \cr\cr
+#' \code{AgeMat} Age at maturity in years (assumed to be age at first parturition - 1) \cr\cr
+#' \code{nages} "Maximum" age, treated as the plus group age. The plus group age can be set equal to the age at maturity +2 years without losing accuracy. Must be greater than \code{AgeMat}.\cr\cr
+#' \code{z} The degree of compensation.  The default value is \code{z = 2.39}.\cr\cr
+#' \code{lambdaMax} Maximum steady rate of increase (population growth rate)
 #' @param InitDepl Starting depletion level
 #' @param ConstantCatch Total bycatch each year, expressed as a vector of length \code{nyears}
 #' @param ConstantF vector (length = \code{nyears}) rate of bycatch each year
-#' @param z The degree of compensation.  The default value is \code{z = 2.39}.
 #' @param nyears Number of years to project
-#' @param nages "Maximum" age, treated as the plus group age. The plus group age can be set equal to the age at maturity +2 years without losing accuracy. Must be greater than AgeMat.
-#' @param lambdaMax Maximum steady rate of increase (population growth rate)
 #'
 #' @return A list containing a matrix \code{N} of numbers at age (dimensions \code{nyears} (rows) x \code{nages} (columns)) and one vector \code{TotalPop} (a vector of length \code{nyears}), containing the number of age 1+ individuals in the population.
-
-# Note, nages = Plus Group Age, and Plus Group Age can = AgeMat+2 without losing accuracy (per AEP 11/30/18)
 #'
 #' @examples
 #' # Generate a time series of abundance for a bowhead whale
-#' dynamics(S0 = 0.944, S1plus = 0.99, K1plus = 9000, AgeMat = 17,
+#' dynamics(lh.params = list(S0 = 0.944, S1plus = 0.99, K1plus = 9000, AgeMat = 17,nages = 25,z = 2.39, lambdaMax = 1.04),
 #'  InitDepl = 0.6, ConstantCatch = NA, ConstantF = rep(0.01, times = 100), 
-#'  z = 2.39, nyears = 100, nages = 25, lambdaMax = 1.04)
+#'  nyears = 100)
 #' @export
-dynamics <- function(S0, S1plus, K1plus, AgeMat, InitDepl, 
-                     ConstantCatch = NA, ConstantF = NA, z, 
-                     nyears, nages, lambdaMax) {
+dynamics <- function(lh.params, InitDepl, 
+                     ConstantCatch = NA, ConstantF = NA, 
+                     nyears) {
+  
+  #S0, S1plus, K1plus, AgeMat, nages, z, lambdaMax
+  # Life history params
+  S0 <- lh.params$S0
+  S1plus <- lh.params$S1plus
+  K1plus <- lh.params$K1plus
+  AgeMat <- lh.params$AgeMat
+  nages <- lh.params$nages
+  z <- lh.params$z
+  lambdaMax <- lh.params$lambdaMax
+  
   # Checks
   if (length(ConstantCatch) > 1 & length(ConstantF) > 1) {
     stop("Cannot have both constant F and constant catch- choose one and set the other to NA!")

R/21_projections.R

@@ -2,14 +2,20 @@
 #'
 #' @description Generates several projections, stochasticity is in the number of catches from year to year
 #' @importFrom stats rlnorm
-#' @param NOut  number of simulations
+#' @param NOut  Number of simulations
 #' @param ConstantBycatch Mean and CV of number of animals killed as bycatch per year (assumed lognormal)
 #' @param ConstantRateBycatch Mean and CV of bycatch rate (assumed normal)
-#' @param InitDepl initial depletion. If obs_CV>0, this is the mean depletion.
-#' @param lh.params - life history parameters as a list. The list must include S0, S1plus, K1plus, AgeMat, nages, z, and lambdaMax
-#' @param nyears number of years to do projections
-#' @param obs_CV observation CV. Default to 1 for simple projections
-#' @return list of outputs from simulations: \code{params} contains parameter values for each trajectory as a matrix; \code{trajectories} contains simulation outputs as a matrix; \code{fishing.rates} contain the bycatch rates for each year in each simulation as a matrix; \code{InitDepl} returns the initial depletion for the projections; \code{ConstantBycatch} provides Catch (total individuals killed in bycatch events per year) and CV of Catch (if the user has specified bycatch as a constant number); \code{ConstantRateBycatch} contains Bycatch Rate (additional mortality from bycatch each year) and CV of ByCatch rate. Other parameters are the same as in the \code{dynamics()} function.
+#' @param InitDepl Initial depletion. If obs_CV>0, this is the mean depletion.
+#' @param lh.params Life history parameters as a list. The list must include S0, S1plus, K1plus, AgeMat, nages, z, and lambdaMax.
+#' @param nyears Number of years to do projections
+#' @param obs_CV Observation CV. Default to 1 for simple projections
+#' @return A list of outputs from simulations: \cr\cr
+#' * \code{params} contains parameter values for each trajectory as a matrix; \cr
+#' * \code{trajectories} contains simulation outputs as a matrix; \cr
+#' * \code{fishing.rates} contain the bycatch rates for each year in each simulation as a matrix; \cr
+#' * \code{InitDepl} returns the initial depletion for the projections; \cr
+#' * \code{ConstantBycatch} provides Catch (total individuals killed in bycatch events per year) and CV of Catch (if the user has specified bycatch as a constant number); \cr
+#' * \code{ConstantRateBycatch} contains Bycatch Rate (additional mortality from bycatch each year) and CV of ByCatch rate. Other parameters are the same as in the \code{dynamics()} function.
 #'
 #' @examples
 #' projections(
@@ -33,7 +39,7 @@ projections <- function(NOut,
                         obs_CV = 0) {
   trajectories <- rep(0, times = nyears)
   fishing.rates <- rep(0, times = nyears)
-  params <- 0 # S0, S1plus, fmax, K1plus, z # also need: nyears = 50, nages = 15
+  params <- 0 
 
   # Life history params
   S0 <- lh.params$S0
@@ -45,23 +51,23 @@ projections <- function(NOut,
   lambdaMax <- lh.params$lambdaMax
 
   # Checks for parameter values
-  if (S0 >= 1) {
+  if (lh.params$S0 >= 1) {
     warning("Calf/pup survival must be between 0 and 1")
   }
-  if (S1plus >= 1) {
+  if (lh.params$S1plus >= 1) {
     warning("Adult survival must be between 0 and 1")
   }
-  if (AgeMat > nages) {
+  if (lh.params$AgeMat > nages) {
     warning("Age at maturity must be less than plus group age")
   }
+  if (lh.params$lambdaMax < 1) {
+    warning("LambdaMax should be greater than 1.
+                         Typical values are lambdaMax = 1.04 for cetaceans and lambdaMax = 1.12 for pinnipeds.")
+  }
   if (!is.na(ConstantBycatch$Catch) & !is.na(ConstantRateBycatch$Rate)) {
     warning("You cannot provide both bycatch as a whole number and as a rate.
          Please specify either ConstantBycatch or ConstantRateBycatch.")
   }
-  if (lambdaMax < 1) {
-    warning("LambdaMax should be greater than 1.
-                         Typical values are lambdaMax = 1.04 for cetaceans and lambdaMax = 1.12 for pinnipeds.")
-  }
 
   set.seed(123)
 
@@ -79,9 +85,9 @@ projections <- function(NOut,
       sdlog.catches <- sqrt(log(ConstantBycatch$CV^2 + 1)) # should be ~cv when sd is low (<0.5)
       catch.vec <- stats::rlnorm(n = nyears, meanlog = log(ConstantBycatch$Catch), sdlog = sdlog.catches)
       traj <- dynamics(
-        S0 = S0, S1plus = S1plus, K1plus = K1plus, AgeMat = AgeMat, InitDepl = InitVal / K1plus,
+        lh.params = lh.params,  InitDepl = InitVal / lh.params$K1plus,
         ConstantCatch = catch.vec,
-        z = z, nyears = nyears, nages = nages, lambdaMax = lambdaMax
+        nyears = nyears
       )$TotalPop # TotalPop is 1+ component
       param.vec <- unlist(lh.params)
       if (NTries < NOut) {
@@ -107,10 +113,10 @@ projections <- function(NOut,
       }
       fishing.rates <- rbind(fishing.rates, f.vec)
 
-      traj <- dynamics(
-        S0 = S0, S1plus = S1plus, K1plus = K1plus, AgeMat = AgeMat, InitDepl = InitVal / K1plus,
+      traj <- dynamics(lh.params = lh.params,
+        InitDepl = InitVal / K1plus,
         ConstantF = f.vec,
-        z = z, nyears = nyears, nages = nages, lambdaMax = lambdaMax
+        nyears = nyears
       )$TotalPop # 1+ component of population
 
       param.vec <- unlist(lh.params)

R/28_extract_df.R

@@ -1,6 +1,6 @@
-#' Extract results from `dynamics()` function as a dataframe
+#' Extract results from `projections()` function as a dataframe
 #'
-#' @description #unnest results from dynamics() and turn them into a dataframe that you can bind together
+#' @description Unnest results from projections() and turn them into a dataframe that you can bind together
 #' @param x Outputs from a call to projections(), which are a named list
 #'
 #' @return a dataframe with population trajectories, initial depletion, and bycatch levels that can be printed as a table.