In mairindeith/PVAInvasR_RPackage: PVA to assess invasive species control methods.

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(PVAInvasR)

For a simpler introduction to PVAInvadR, check out the quick start guide vignette.

Introduction

PVAInvadR is a population viability analysis (PVA) simulator designed to quickly and accurately reflect uncertainty in the fate of invasive populations. The software evaluates how candidate removal strategies that target different life stages can influence population persistence.

Moreover, given a suite of alternative removal strategies--differing by the types of gears (or combinations of gears) and effort regimes--PVAInvadR can run a PVA with each removal strategy and compare their efficacy with a multi-criteria decision table. The removal stragies can be compared according to the probability of eradicating the targetted species, the annual and cumulative cost of implementing the gears, and the probable range of target species abundance after the full simulation period.

The goal of the software is to provide rapid assessment of different strategies to assist workshop-type discussions of invasive species removal, thereby encouraging discussion among stakeholders involved in invasive species removal.

This user guide begins by describing the model's internal functions in Model description (including the parameters needed to run the model) and closes with how to run PVAInvasR functions

Model description {#modeldesc}

The PVA model can be used on any species with indeterminate growth (for example, most fish, reptiles, amphibians, and invertebrates). Populations are split into two life stages: pre-recruits and recruited juveniles/adults. Pre-recruits can be further split into stanzas or cohorts -- fish in each stanza are subject to density-dependent mortality. Users can define the strength of this dependence for each stanza. Recruited animals, on the other hand, are age structured and subjected to length-dependent mortality.

However, this approach can also be catered somewhat to specific case studies. For example, if density-dependence persists after recruitment, age at recruitment can be adjusted and the number of stanzas increased to overlap with juvenile life stages.

Biological and control parameters {#params}

The PVA model (called by the PVA() function) simulates population trajectories given a set of input biological and control parameters. It also depends on parameters that define the behaviour of the PVA simulation; for example, for how many iterations and over what period of time should the model be run?

Below are tables of the biological, control, and PVA parameters necessary to run a PVA simulation.

library(knitr)
library(dplyr)
library(kableExtra)

param_names <- c(
  'species',
  'A',
  'AR',
  'nS',
  'nT',
  'dt',
  'n_sim',
  'n_gear',
  't_start_R', #
  't_start_A', #
  # Population parameters
  'V0',
  'reck',
  'p_can',
  'K',
  'afec',
  'Wmat',
  't_spn',
  'Ms',
  'Bs',
  'V1',
  'bet',
  'cann_a',
  'sd_S',
  'UR', #
  'UA', #
  'samp_A',
  'E_R', #
  'E_A', #
  'C_f_R', #
  'C_f_A', #
  'C_E_R', #
  'C_E_A', #
  'r',
  'G',
  'v_a', #
  'v_b', #
  'v_c', #
  'v_d', #
  'init_NA'
)
notes <<- c(
  'Species abbreviation',
  'Maximum age (the age where only 1% of fish survive when unfished)',
  'Age at recruitment (years)',
  'Number of pre-recruit stanzas',
  'Length of simulation (number of years)',
  'Time-step length as a proportion of one year (<1 means more than 1 time step per \n year)',
  'Number of simulations',
  'Number of gears applied to recruited animals',
  'Time-step when removal starts for each pre-recruit stanza',
  'Time-step when removal starts for each gear used on recruited animals',
  'Range of unfished recruited abundance',
  'Recruitment compensation ratio (difference in juvenile survival between unfished \n and near-zero density)',
  'Proportion of mortality at equilibrium that is due to cannibalism',
  'von Bertalanffy growth parameter',
  'Fecundity multiplier on weight',
  'Minimum weight at maturity (as a proportion of maximum weight)',
  'Time-steps of the year when spawning occurs (indicated as a 0/1 switch)',
  'Stanza-specific maximum survival',
  'Stanza-specific maximum available habitat',
  'Initial number of vulnerable fish in the population',
  'Hyperstability parameter',
  'Age at which cannibalism begins',
  'Standard deviation in recruitment across years',
  'Proportion of pre-recruit animals removed by gear per unit effort',
  'Proportion of recruited animals removed by gear per unit effort',
  'Time step(s) in the year that each gear is fished',
  'Effort expended by each pre-recruit capture gear',
  'Effort expended by each recruit capture gear',
  'Fixed cost associated with each pre-recruit capture gear',
  'Fixed cost associated with each recruit capture gear',
  'Effort-based cost associated with each pre-recruit capture gear',
  'Effort-based cost associated with each recruit capture gear',
  'Discount rate (optional)',
  'Generation time (optional)',
  'Logistic ascending slope of removal gear for recruited animals (as a proportion \nof Linf)',
  'Ascending length at 50% selectivity of removal gear for recruited animals (as a \nproportion of Linf)',
  'Logistic descending slope of removal gear for recruited animals (as a proportion \nof Linf)',
  'Descending length at 50% selectivity of removal gear for recruited animals (as a \n proportion of Linf)',
  '(Optional) Initial age structure in the population'
)

file_out <- data.frame(
  Parameters = param_names,
  # values = "",
  Description = notes
)

# Rows for biological params
biol_rows <- c(1:4, 11:23, 39)
# Rows for control params
control_rows <- c(8:10, 24:38)
# Rows for PVA params
pva_rows <- c(5:7)

knit_table <- function(df, caption){
  if (is_html_output()) {
    df %>%
      kable("html", escape = F, caption = caption) %>%
      kable_styling()
  } else {
    df <- data.frame(lapply(df, function(x) {gsub("<br>", "\n", x)}), stringsAsFactors = F)

    df %>%  
      mutate_all(linebreak) %>%
      kable("latex", booktabs = T, escape = F, caption = caption)  
  }
}


knit_table(file_out[biol_rows,], caption = "Biological parameters")
knit_table(file_out[control_rows,], caption = "Control parameters")
knit_table(file_out[pva_rows,], caption = "PVA parmaters")

# pandoc.table(file_out[biol_rows,], wrap text in  align='l', row.names = NA, caption = "Biological parameters")
# pandoc.table(file_out[control_rows,], wrap text in  align='l', row.names = NA, caption = "Control parameters")
# pandoc.table(file_out[pva_rows,], wrap text in  align='l', row.names = NA, caption = "PVA parmaters")

When these user-input parameters are used to derive other parameters below, user-defined parameters are marked with an asterisk ($$). For example, calculations that use maximum age, A, is indicated with $A^{\text{}}$.

Derived parameters and population simulations

Length, weight, and fecundity at age

Using the user-input parameters (those listed in the above table), PVAInvadR simulates the population dynamics of the exploited population. First, PVAInvadR derives stock-recruit parameters from calculations of equilibrium population structure when the population is unexploited.

Length, weight, and fecundity at age $a$ ($l_a$, $w_a$, and $f_a$, equations 1-3) are calculated first, then instantaneous length-based mortality based on the Lorenzen model (such that mortality is inversely related to length; equation 4, Lorenzen 2000).

\begin{equation} \label{eq:1} l_a = 1 - e^{-K^{\text{\text{*}}} a} \end{equation}

\begin{equation} \label{eq:2} w_a = l_a^3 \end{equation}

\begin{equation} \label{eq:3} f_a = afec^{\text{}} (w_a - Wmat^{\text{}}) \end{equation}

Natural mortality

Adult natural mortality

Natural mortality rate is asymtomatic and calculated from $A^{\text{\text{}}}$, the maximum age observed in the population. This maximum age is assumed to have a lifetime survivorship of 1% (i.e. 1% of recruited animals live to age = $A^{\text{\text{}}}$). $M_{\infty}$, minimum instantaneous mortality, is back-calculated similar to the method of Hoenig (1983).

\begin{equation} \label{eq:4} M_{\infty} = \frac{ln(0.01)K^{\text{\text{}}}}{ln(l_{a=AR^{\text{\text{}}}}) - ln \left [ l_{a=AR^{\text{\text{}}}} + e^{\frac{K^{\text{\text{}}}(A^{\text{}} - AR^{\text{\text{}}})}{dt^{\text{\text{*}}}}}-1 \right ]} \end{equation}

From this, survival of each age group follows Lorenzen 2000.

\begin{equation} \label{eq:5} S_a = \left ( \frac{l_a}{l_a + e^{K^{\text{\text{}}}dt^{\text{}}}-1}\right )^\frac{M_{\infty}}{K^{\text{\text{}}}} \end{equation}

Then, survivorship to each age ($lx_a$) is calculated:

\begin{equation} \label{eq:6} lx_a = \begin{cases} 1 & \text{ if } a = AR^{\text{\text{}}}\ \prod_{a=AR^{\text{}}}^{a-1} S_a & \text{ if } a > AR^{\text{\text{*}}} \end{cases} \end{equation}

The product of survivorship and fecundity-at-age for females then provides equilibrium spawners per recruit based on the ages when spawning occurs (i.e. where $spn_a = 1$; non-spawning ages are indicated with $spn_a = 0$).

\begin{equation} \label{eq:7} \varphi_0 = \sum^{A^{\text{\text{*}}}}_{a=1} \frac{lx_a f_a spn_a}{2} \end{equation}

The age indicator is used to calculate Beverton-Holt recruitment parameters (Walters and Martell, 2004).

For pre-recruits, Beverton-Holt parameters are calculated for each stanza, $s$. The Beverton-Holt parameters are influenced by stanza-specific relative mortality, $M_s$, habitat capacity (which scales mortality by density of competitors within a cohort; (Walters and Korman, 1999; Pine et al., 2013), and cannibalism if it occurs.

This form of the Beverton-Holt model is from Walters and Korman's (1999) reformulation:

\begin{equation} \label{eq:8} R = \frac{N_0 e^{-M_0}}{1 + \frac{M_1}{M_0} \left ( 1 - e^{-M_0} \right ) } \end{equation}

Here, $e^{-M_0}$ is maximum survival at low spawning stock and $\frac{M_1}{M_0} \left ( 1 - e^{-M_0}\right )$ is the carrying capacity parameter.

When cannibalism occurs, $M_0$ can be modified by combining stanza-specific mortality without cannibals ($\rho_s$) with the product of likely cannibalistic age classes $C$ and the cannibalism dependent parameter, $\tau$.

\begin{equation} \label{eq:9} M_{0,s} = \rho_s + \tau C = \rho_s + \tau \left [ R_0 \sum_{a = \text{cann_a}^{\text{*}}} \left ( lx_a spn_a \right ) \right ] \end{equation}

\begin{equation} \label{eq:10} \tau_s = M_s^0 \frac{\text{p_cann}}{R_0 \sum_{a = \text{cann_a}^{\text{*}}} \left ( lx_a spn_a \right )} \end{equation}

Alternatively, mortality can be expressed by:

\begin{equation} \label{eq:11} M_{0,s} = -ln (\alpha_s) \end{equation}

Stanza-specific $\alpha$ and $\beta$ parameters are calculated as:

\begin{equation} \label{eq:9} \alpha_s = \frac{K}{\varphi_0} e^ {\left ( \frac{M_s^{\text{\text{}}}}{\sum M_s^{\text{\text{}}}} \right )} \end{equation}

\begin{equation} \label{eq:10} \beta_s = Bs^{\text{\text{}}} \frac{K^{\text{\text{}}}-1}{R0^{\text{\text{*}}}\varphi_0 \sum_{s'} \left ( \beta_{s'} \prod_{s''=0}^{s''=s'-1} \alpha_{s''}\right )} \end{equation}

This results in a Ricker-like recruitment (e.g. overcompensation at high spawner abundance).

Initializing populations

Initially high abundance populations (starting at >90% of carrying capacity)

There are two ways to determine the number of recruited animals in the first year, $t = 1$. If abundance is close to carrying capacity (i.e. $V_1 \geq 0.9 R_0 \sum lx_a$), abundance is randomly allocated among age classes (assuming a multinomial distribution):

\begin{equation} \label{eq:16a} N_{t=1,a} \sim Mult(V_1^*, lx_a, e^{\epsilon_a}); \epsilon_a \mathcal{N}(0, \sigma_R) \end{equation}

Random year class strength is given by $\epsilon_a$ above. This is a useful approximatino if the population is close to carrying capacity. However, it will under-represent early year classes if the population is still growing (owing to the steady state assumption). If the initial abundance is low, there is also a chance that all mature year-classes will be empty, causing a lag in egg production in early years; this will result in population growth being low in the first few years (depending on generation time).

Initially low abundance (starting at <90% of carrying capacity)

If initial abundance is less than 90% of carrying capacity, initial age structure is created by deterministically simulating the population from low abundance (i.e. with $10e^{-15}$) recruits until the population reaches the user-defined starting abundance $V_1^*$.

First, the abundance of each age class is calculated:

\begin{equation} \label{eq:16b} N_{t=1,a}^{i} = N_{t-1, a-1}^{i} S_{a-1} \end{equation}

where $N^{i}$ indicates the numbers in initialization. The product of numbers-at-age and fecundity-at-age (equation 3) are used to calculate egg production in the first year ($E_{t=1}$; using the stage-independent Beverton-Holt function).

\begin{equation} \label{eq:17} N_{t, a=1}^i = \frac{E^i_{t-1} e^{-(\rho + \tau C_t)}}{1 + \frac{e^{-(\rho + \tau C_t)}}{M_1} \left ( 1 - e^{-(\rho + \tau C_t)} \right )} \end{equation}

\begin{equation} \label{eq:18} M_1 = \frac {\left ( K^{\text{}} - 1 \right ) ln \left( \frac{K^{\text{}}}{\varphi_0} \right )}{R_0 \left( K^{\text{*}} - \varphi_0 \right ) } \end{equation}

\begin{equation} \label{eq:19} \rho = -ln \left ( \frac{K^{\text{*}}}{\varphi_0} \right ) \left ( 1 - p_{cann} \right ) \end{equation}

\begin{equation} \label{eq:20} \tau = -ln \left ( \frac{K^{\text{*}}}{\varphi_0} \right ) \frac{p_{cann}}{R_0 \sum_{a = cann_a} \left ( lx_a spn_a \right )} \end{equation}

Equations 18-20 are stage-independent analogues to equations 11, 14, and 15.

Subsequent egg production and recruitment

Regardless of how the population is initiated, egg production in the first year is given by:

\begin{equation} \label{eq:21} E_{t=1,a} = \sum_{i=1}^{N_t=1} f_a \end{equation}

Recruitment in the following years is based on survival through each pre-recruit stanza.

Mortality

Survival from one stanza to the next is given by the product of focused removals (i.e. capture) and the stanza-specific stock-recruitment function:

\begin{equation} \label{eq:22} N_{t, s+1} \sim Binomial \left ( N_{t, s}, \frac{e^{-F_s -M_{s,t}^0 + \psi_{s,t}}}{1 + \frac{M_s^0}{M^0_{s,t}} \left( 1-e^{-M^0_{s,t}} \right)} \right ); \psi_{s,t} = N(0, \sigma_r) \end{equation}

Maximum survival of each pre-recruit stanza is a function of $\rho$, mortality independent of cannibals, \begin{equation} \label{eq:23} M_{s,t}^0 = \rho_s + \tau_s \sum_{cann_a} (N_{t,a}) \end{equation}

and the effort used in each time step, $F_s$ given by:

\begin{equation} \label{eq:24} F_s = q_s E_s \end{equation}

Survival of each recruited age-class in each timestep is described by the Baranov equation (Ricker 1975):

\begin{equation} \label{eq:25} N_{t+1, a+1} \sim Binomial \left ( N_{t,a}, e^{-Z_{t,a}} \right) \end{equation}

Instantaneous mortality for each age and timestep, $Z_{t,a}$ is a function of natural mortality, $M_a$, and capture mortality, $F_{t,a}$ for each gear, $g$.

\begin{equation} \label{eq:26} Z_{t,a} = M_a + F_{t,a} \end{equation}

Removal mortality

Specifically, capture mortality is calculated as the sum-product of each gear's effort in time t, $E_{t,g}$ -- set to 0 except for where users indicate the gear is used -- and the selectivity of the gear to each age class, $v_{a,g}$.

\begin{equation} \label{eq:27} F_{t,a} = \sum_{g=1}^{n_g} q_{N,g} E_{t,g} v_{a,g} \end{equation}

Each control gear is described by a density-dependent constant of proportionality, $q_{N,g}$. This is often termed catchability in fisheries literature. Catch-per-unit effort is related to abundance, or in other terms, removal effort compared to instantaenous removal mortality (Arreguin-Sanchez 1996). There are many reasons why catchability may be density dependent, but a common mechanism in aquatic invasive species is due to change in individual home range size as the population grows (e.g. the population is more densely packed as it fills the spatial range that is occupied by the population).

The invaded area increases non-linearly with increased abundance according to

\begin{equation} \label{eq:28} A_N = \gamma N^{\beta} \end{equation}

$\gamma$ describes the (non-linear) density independent rate of increase in the area occupied by the population. $\beta$ is the rate at which animals reduce density with population abundance.

Catchability can then be calculated as the proportion of an area affected when swept by a removal gear.

\begin{equation} \label{eq:29} q_N = \frac{\delta}{A_N} \end{equation}

where $\delta$ is the area swept by that gear with one unit of effort.

Equations 28 and 29 can be combined into the standard hyperstability equation.

\begin{equation} \label{eq:30} q_N = qN^{-\beta} \end{equation}

where $q$ is the probability of capture when the population is at low abundance.

The two parameters $\phi$ and $\beta$ can be determined by considering how catch per unit effort would change as the population grows.

For example, imagine a system where a gear is fished with one unit of effort at low population abundances at

two population sizes

(e.g. $N_1 = 50$ and $N_2 = 200$) resulted in $CPUE_1 = 5$ and $CPUE_2 = 15$ you can calculate $\beta$ as:

\begin{equation} \label{eq:31} \beta = 1 - \frac{ln(5)-ln(15)}{ln(50)-ln(200)} = 0.207 \end{equation}

Note that $q_N = \frac{\text{CPUE}}{N}$, and Equation 31 can be substituted into Equation 30 to give $q = 0.225$.

PVAInvadR repeats these steps every timestep defined by the user for each iteration of the simulation.

Running the model and use for decision support {#modelrun}

When beginning a PVA session, the user can decide to consider one or multiple capture scenarios to control the invasive population following the calculations above.

Briefly, here are the basic steps to comparing multiple control scenarios using PVAInvadR with more details below:

1. Provide biological and control parameters to the model;
2. (Optional) Supply these parameters to PVA() as the "base model";
3. Define multiple control scenarios with decision_setup();
4. Compare these scenarios with decision();
5. Examine how well these scenarios perform when simulated under varying biological parameters with rank_uncertainty() (which assesses the impact of parameter uncertainty).

Upload parameters

PVAInvadR requires two types of parameters: biological and control parameters (described above with full list and definitions).

These parameters are uploaded in a .csv file saved on your hard drive. You can start by downloading a template .csv that includes parameter names, a blank column to enter parameter values, and descriptions of these parameters.

Download a template with pva_template("FILE/PATH"), where FILE/PATH is the location where your template should be saved.

When filling out the template, a few formatting rules must be followed:

• The number of stanzas influences how many pre-recruit parameters are needed. For example, if you are modelling a population with 5 stanzas, there must be 5 parameter values provided for each stanza-specific parameter (Ms, Bs, qR, ER, CFR, and CER). Put all parameter values in the same space, but separate them with a semicolon, ";"
• Similarly, for each control gear, gear-specific parameters must be provided in a single .csv cell and separated with a semicolon, ";".

[[IMAGE HERE]] - Any parameter where you would like to list multiple values (like cann_a, for example), these must again be separated with a ";". - Parameter names must stay the same as in the template - otherwise they will not be properly loaded by PVAInvadR.

All parameters in the template must be provided except for initial age structure. This can be left blank and will be initialized by PVAInvadR. Control options should be relatively easy to define, but biological parameters may be more difficult to estimate. We suggest a literature survey, researching parameters for similar species, or forming a good biological understanding of the species can help.

Once you have filled in the PVA parameters and saved the values to a .csv file, upload this to your R session (after replacing NEW/FILE/PATH with the path to your filled .csv):

pva_params <- load_pva_parameters("NEW/FILE/PATH")

Now, pva_params contains a named list of parameters from your .csv. You can check the value of any parameter using R's $-notation (for example, to check the value for the proportion of cannibals, use pva_params$p_can). You can also do this to modify values in the R session (but note that this will not modify the source .csv and changes to PVA parameters may be lost).

Evaluating a single removal strategy with PVA

Evaluating multiple removal strategies with decision

Evaluating multiple removal strategies under biological uncertainty with rank_uncertainty

References {#refs}

• Arreguin-Sanchez, F., 1996. Catchability: assessment a key parameter for fish stock. Rev. Fish Biol. Fish. 6, 221–242.
• Hoenig, J.M., 1983. Empirical use of longevity data to estimate mortality rates. Fish. Bull. 82, 898–903. doi:10.2307/1940001
• Lorenzen, K., 2000. Allometry of natural mortality as a basis for assessing optimal release size in fish-stocking programmes. Can. J. Fish. Aquat. Sci. 57, 2374–2381.
• Pine, W.E.I., Healy, B., Smith, E.O., Trammell, M., Speas, D., Valdez, R., Yard, M., Walters, C., Ahrens, R., Vanhaverbeke, R., Stone, D., Wilson, W., 2013. An individual-based model for population viability analysis of humpback chub in Grand Canyon. North Am. J. Fish. Manag. 33, 626–641.
• Ricker, W.E., 1975. Computation and interpretation of biological statistics of fish populations. Bull. Fish. Res. Board Canada 1913890, 382pp.
• Walters, C., Korman, J., 1999. Linking recruitment to trophic factors: revisiting the Beverton-Holt recruitment model from a life history and multispecies perspective. Rev. Fish Biol. Fish. 9, 187–202.
• Walters, C., Martell, S., 2004. Fisheries ecology and management. Princeton University Press, Princeton.

mairindeith/PVAInvasR_RPackage documentation built on Feb. 1, 2023, 3:17 p.m.