In mcsiple/mmrefpoints: Project Marine Mammal Populations and Calculate Reference Points

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

Population model

Population projections are generated using a single-sex age-structured model. In this projection model, the number of calves or pups born each year is density dependent, with the extent of density dependence being a function of the number of mature adults $\tilde N$, the fecundity (pregnancy rate) at pre-exploitation equilibrium $f_0$, the maximum theoretical fecundity rate $f_{max}$, the degree of compensation $z$, and the abundance of individuals aged 1+ ($N_{y+1}^{1+}$) relative to carrying capacity $K^{1+}$.

The number of individuals age 1 and older is a function of calf/pup survival, $S_0$, and the survival rate of adults, $S_{1+}$, and removals due to bycatch mortality, $C_y$ (all symbols are defined in Table 1).

$$ \tag{Eq 1} N_{y+1,a} = \begin{cases} \tilde N_{y+1}\bigg(f_0 + (f_{max}-f_0)\bigg[1-\bigg(\frac{N_{y+1}^{1+}}{K^{1+}}\bigg)^z\bigg]\bigg) & a=0\ N_{y,0}S_0 & a=1 \ (N_{y,a-1}-C_{y,a-1})S_{1+} & 2\leq a<x \ (N_{y,x-1}-C_{y,x-1})S_{1+} + (N_{y,x}-C_{y,x})S_{1+} & a=x \ \end{cases} $$

where $N_{y,a}$ is the number of animals of age $a$ at the start of year $y$. The number of animals of age $a$ dying due to bycatch mortality during year $y$, $C_{y,a}$, is removed uniformly from the 1+ component of the population (calves and pups are assumed not to die due to bycatch), i.e.:

$$ \tag{Eq 2} C_{y,a}= \begin{cases} 0 & a=0 \ \frac{C_y N_{y,a}}{N_y^{1+}} & a>0 \ \end{cases} $$

Given this assumption, our model will not adequately characterize cases where bycatch mortality occurs predominantly among calves/pups.

If the user specifies a constant rate of bycatch mortality, individuals are removed from the population according to a bycatch mortality rate $E_y$ and vulnerability (which is assumed constant through time, and uniform on age 1+ individuals):

$$ \tag{Eq 3} C_{y,a}= \begin{cases} 0 & a=0 \ N_{y,a} E_{y} & a>0 \ \end{cases} $$ If bycatch is specified by the user as a number of individuals per year, annual bycatch rates $E_y$ are lognormally distributed with mean $log(C_{user})$ and standard deviation $\sigma_E$ of $\sqrt{log(CV_C^2 + 1)}$, where $C_{user}$ is the bycatch in numbers per year defined by the user and $CV_C$ is the bycatch CV defined by the user. If bycatch is specified as a constant rate, $E_y$ is beta-distributed with mean $E_{user}$ and standard deviation $E_{user} \cdot CV_E$.

Table 1. Symbols included in the projection model.

Symbol=c(#User-specified parameters
         "$S_0$",
         "$S_{1+}$",
         "$x$",
         "$\\lambda_{max}$",
          "$a_p$",
         "$E_y$",
         # Parameters derived from user-specified parameters
         "$f_0$",
         "$f_{max}$",
         "$K^{1+}$",
         "$z$",
         # Derived variables
         "$\\tilde N_{y,a}$", # need double dash for using in kableExtra
         "$N_{y,a}$",
         "$N_y^{1+}$",
         "$C_{y,a}$")

Description=c(#User-specified parameters
              "Pup or calf survival",
              "Survival of individuals aged 1+",
              "Plus-group age",
              "Maximum steady rate of increase (population growth rate)",
              "Age at first parturition",
              "Bycatch mortality rate in year $y$ (specified by the user or computed from the user-specified total bycatch mortality)", 
              # Parameters derived from user-specified parameters
              "Fecundity (pregnancy rate) at pre-exploitation equilibrium",
              "Maximum theoretical fecundity (pregnancy rate)",
              "Carrying capacity in terms of 1+ component of the population",
              "Degree of compensation",
              # Derived variables
              "Number of mature animals of age $a$ at the start of year $y$",
              "Number of animals of age $a$ at the start of year $y$",
              "Number of animals aged 1 and older at the start of year $y$",
              "Mortality due to bycatch of animals of age $a$ in year $y$")

Symbol_definitions <- data.frame(cbind(Symbol,
                                      Description))
knitr::kable(format="html", #html  - change to pandoc for ms word 
             Symbol_definitions, 
             escape = FALSE) %>% 
            kable_styling('bordered') %>%
            pack_rows('User-specified parameters', 1, 6, label_row_css = "font-style: italic;border-bottom: 1px solid black;") %>%
            pack_rows('Parameters derived from user-specified parameters', 7, 10, label_row_css = "font-style: italic;border-bottom: 1px solid black;") %>%
            pack_rows('Derived variables', 11, 14, label_row_css = "font-style: italic;border-bottom: 1px solid black;")

Performance measures

There are three primary performance measures included in the outputs:

• the probability that the abundance of age 1+ animals exceeds the Maximum Net Productivity Level (MNPL) after 50 and 100 years;
• the abundance of age 1+ animals compared to 1+ carrying capacity after 10, 20, and 50 years; and
• the abundance of age 1+ animals compared to the abundance of age 1+ animals under a no-bycatch scenario after 10, 20, and 50 years.

These performance measures are all related to population recovery. Wade (1998) identifies a performance ‘standard’ of a 95% probability of recovery to MNPL after 100 years for a population initially at 30% of its carrying capacity, with an MNPL of 0.5K. MNPL is the lower bound for the U.S. Marine Mammal Protection Act Optimum Sustainable Population (OSP), which is defined as the “number of any animals which will result in the maximum productivity of the population or the species, keeping in mind the carrying capacity of the habitat and the health of the ecosystem of which they form a constituent element” (16 USCS § 1362 (9)). In the US management scheme, marine mammal stocks are considered depleted when they are below OSP (16 USCS § 1362 (1B)).

Single values reported in tables are the medians from the number of projections.

Solving for bycatch rate

The app can solve for the bycatch level that would result in recovery to a given proportion of K in a specific amount of time. The rebuilding probability is calculated as the proportion of simulations in which the population has recovered to the rebuilding goal by the rebuilding year. The “Solve for bycatch” tab uses root finding (the function uniroot() in our code) to find the bycatch mortality rate that minimizes the difference between the desired recovery probability and the recovery probability under bycatch mortality rate $E$.

Maximum theoretical fecundity

The maximum theoretical fecundity rate is calculated based on the population size when the population is growing at its maximum population growth rate $\lambda_{max}$. Parts of this derivation can be found in @butterworth_scientific_1992 and @breiwick_population_1984. Starting with the population model for the female component of the population:

$$ \tag{Eq 4} N_{t+1} = N_t(1-E)S_{1+}+N_{t-a_p}f q_f S_0 (S_{1+})^{a_p -2}(1-E)^{a_p-1-a_r} $$

where $a_p$ is the age at first parturition (assumed to be one year after the age at maturity), $f$ is the fecundity (sometimes expressed as the product of pregnancy rate and pup survival), $q_f$ is the proportion of calves/pups that are female and $a_r$ is the age at which the mammals first suffer any bycatch mortality (in our case, this is equal to 1 year). When the population is growing at its fastest rate, fishing is zero, and $N_{t+1}=N_t \lambda_{max}$. In our case $q_f = 1$ because we are modeling all adults, instead of just mature females. Thus, the above equation becomes:

$$ \tag{Eq 5} \lambda_{max}N_t=N_t S_{1+}+f_{max} N_{t-a_p-1} \lambda {max}^{a_p-1} S_0 S{1+}^{(a_p-2)} $$ Solving for $f_{max}$ gives the maximum theoretical fecundity as a function of $\lambda_{max}$, survival, and age at first parturition:

$$ \tag{Eq. 6} f_{max} = \frac {\lambda_{max}^{a_p -1} - \lambda_{max}^{{a_p -2}}} {S_0 {S_{1+}^{a_p-2}}} $$ This is referred to as $p$ in @butterworth_scientific_1992.

Maximum net productivity level

To calculate the maximum net productivity level ($MNPL$) given $z$, we first calculate the sustainable yield $C$ as a function of bycatch mortality rate $E$.

$$ \tag{Eq. 7} C = E \tilde B(E)\tilde P(E) $$ where $\tilde B (E)$ is the normalized (to numbers at unfished equilibrium) recruitment when the bycatch rate is fixed at $E$ and $\tilde P (E)$ is the equilibrium number of "recruited" (age 1+) animals per calf/pup when the bycatch mortality rate is fixed at $E$. The normalized recruitment $\tilde B(E)$ is calculated as follows:

$$ \tag{Eq. 8} \tilde B(E) = \bigg(1 - \frac{1-f_0 \tilde N(E)} {Af_0\tilde N (E)}\bigg)^{1/z} \bigg(\frac{\tilde B(0)\tilde P(0)}{\tilde P(E)}\bigg) $$ where $f_0 = \frac{1}{\tilde N(0)}$, $\tilde N(E)$ is the number of animals at the age of first parturition and older (i.e., reproducing animals) per recruit at no-bycatch-mortality equilibrium, and $A$ is the Pella-Tomlinson resilience parameter ($A=\frac{f_{max}-f_0}{f_0}$ ; @punt_a._e._annex_1999). $\tilde B(0)$ is assumed to be equal to 1, because all calculations are per-recruit.

The number of reproducing animals per recruit at exploitation rate $\boldsymbol{E}$ is the sum of adult animals per recruit $\tilde N_a$ from the age at first parturition $a_p$ to the plus group age $x$:

$$ \tag{Eq. 9} \tilde N(E) = \sum_{a=a_p}^{x} \tilde N_{a}(E) $$

We solve for the bycatch mortality rate at which $\tilde C$ is maximized (i.e., where $\frac{dC}{dE}$ is zero). This is $MSYR$, which is analogous to $F_{MSY}$ in fisheries.

The number of 1+ animals per recruit at bycatch mortality rate $E$, $\tilde P(E)$ is defined as: $$ \tag{Eq. 10} \tilde P(E)=\sum_{a=1}^{x} \tilde N_{a}(E) $$

where $\tilde N_{0,a}(E)$ is the numbers per recruit at each age $a$ given a stable age structure:

$$ \tag{Eq. 11} \tilde N_{0,a}(E) = \begin{cases} 1 & a=0 \ S_0[S_{1+}(1-E)]^{a-1} & 1\leq a<x \ \frac{S_0[S_{1+}(1-E)]^{x-1}}{1-[S_{1+}(1-E)]} & a=x \ \end{cases} $$

$MSYR$ is the value of $E$ for which the derivative of $C$ with respect to $E$ is zero, which we determined through numerical differentiation:

$$ \tag{Eq. 12} \frac{dC}{dE} = \frac {C(E+0.001) - C(E-0.001)} {0.002} $$

Then the relative abundance that corresponds to $MSYR$, $MNPL$, is determined by calculating the total 1+ population size at $MSYR$ relative to the equilibrium 1+ population size with no bycatch mortality:

$$ \tag{Eq. 13} MNPL = \frac{\tilde P(E=MSYR)\tilde B(E=MSYR)}{\tilde P(0)\tilde B(0)} = \frac{\tilde P(E=MSYR)\tilde B(E=MSYR)}{\tilde P(0)} $$

Parameterization

We assume that the population starts with a stable age structure in year 1 of the projection period (Eq. 11). The numbers at age at the start of projections correspond to a constant bycatch mortality rate $E$, which is calculated by solving the following equation for $E$:

$$ \tag{Eq. 14} \frac{\tilde B(E)\tilde P(E)}{\tilde P(0)}= \frac{N_0^{1+}}{K^{1+}} $$

The initial depletion $\frac{N_0^{1+}}{K^{1+}}$ is based on the history of human-caused mortality for the population, which is provided by the user.

For cases where observation error is given for the initial population size, the starting abundance is drawn from a lognormal distribution with a standard deviation proportional to the observation CV.

Life history types

Each marine mammal life history type presented as a default option in this app corresponds to a unique value of calf/pup survival $S_0$, adult survival $S_{1+}$, age at first parturition, $a_p$, and intrinsic rate of population growth $\lambda_{max}$. These values are presented in Table 2. For computation purposes, we assumed that the plus group age $x$ is two years past the age at maturity ($x=a_p+1$).

Compensation

We solve for the value of the degree of compensation $z$ that corresponds to the value of MNPL provided by the user. This involves solving the equation $\tilde P(E^) \tilde B(E^) = MSYL * \tilde P(0)$ for $z$ where $E^*$ depends on $z$ as outlined above.

Table 2.

lh.sources <- data.frame(stringsAsFactors=FALSE,
                           Type = c("Bowhead whale", "Bottlenose dolphin", "Humpback whale",
                                    "Phocid seal", "Fur seal", "Sea lion",
                                    "Porpoise", "Minke whale",
                                    "False killer whale/killer whale", "Pilot whale", "Right whale"),
                 Representative = c("Balaena mysticetus", "Tursiops truncatus",
                                    "Megaptera novaeangliae", "Phoca vitulina",
                                    "Arctocephalus pusillus pusillus",
                                    "Zalophus californianus", "Phocoena phocoena",
                                    "Balaenoptera bonaerensis", "Orcinus orca",
                                    "Globicephala macrorhynchus", "Eubalaena glacialis"),
                             S0 = c(0.944, 0.865, 0.9, 0.802, 0.77, 0.83, 0.8096, 0.84216,
                                    0.84744, 0.85008, 0.85536),
                         S1plus = c(0.99, 0.951, 0.95, 0.92, 0.88, 0.95, 0.92, 0.957, 0.963,
                                    0.966, 0.972),
                         AgeMat = c(17, 6, 10, 6, 3, 4, 3, 7, 9, 9, 8),
                   #PlusGroupAge = c(25, 10, 15, 8, 10, 5, 7, 9, 11, 11, 10),
                         Source = c("@punt_conserving_2018 and references therein",
                                    "@punt_conserving_2018 and references therein,
                                    except juvenile survival which is set to 0.88($S_{1+}$)", "@punt_conserving_2018,
                                    @arso_civil_variations_2019, @speakman_mark-recapture_2010",
                                    "@punt_conserving_2018; @hastings_sex-_2012",
                                    "@punt_conserving_2018; @butterworth_effects_1995",
                                    "@punt_conserving_2018; @delong_age-_2017",
                                    "@moore_unpublished_2019; G. Vikingsson pers. comm.; olafsdottir_growth_2003",
                                    "@moore_unpublished_2019",
                                    "@moore_unpublished_2019",
                                    "@moore_unpublished_2019",
                                    "@moore_unpublished_2019")
              )


x <- knitr::kable(format="html",
             col.names = c('Type','Representative','$S_0$','$S_{1+}$','Age at maturity','Source'),
             lh.sources, 
             escape = FALSE) %>% 
            kable_styling('bordered')
x <- column_spec(x, column = 2,italic = TRUE)
x

mcsiple/mmrefpoints documentation built on June 17, 2022, 8:41 p.m.