add fixed_percentage control rule (untested)

DESCRIPTION

@@ -35,6 +35,7 @@ Imports:
     gridExtra,
     magrittr,
     MCMCpack,
+    methods,
     parallel,
     patchwork,
     rstan,

NAMESPACE

@@ -13,6 +13,7 @@ export(cr_param)
 export(depensation)
 export(dtlnorm)
 export(dtnorm)
+export(fixed_percentage)
 export(forward_pbr)
 export(forward_rla)
 export(get_streaks)
@@ -61,6 +62,7 @@ importFrom(dplyr,summarize)
 importFrom(dplyr,ungroup)
 importFrom(gridExtra,tableGrob)
 importFrom(magrittr,"%>%")
+importFrom(methods,is)
 importFrom(parallel,detectCores)
 importFrom(rstan,Rhat)
 importFrom(rstan,extract)

R/fixed_percentage.R

+#' A function to carry out simulations according to the fixed percentage rule
+#'
+#' @param poplist a list; the output from one call to function \code{\link{pellatomlinson_pbr}} or \code{\link{pellatomlinson_rla}}
+#' @param bycatch_variation_cv a positive scalar; variability in bycatch
+#' @param distribution a character string; must be one of "truncnorm",
+#' "lognormal", or "gamma" to simulate catch numbers.
+#' @param frequency a count scalar; frequency of new surveys. Default to 6, i.e.
+#' a new survey every 6 years.
+#' @param horizon a count scalar; number of years to run the management procedure.
+#' Default to 100 years.
+#' @param p a positive scalar between 0 & 1; the percentage to evaluate. Default to 0.017.
+#' @param p_discount a positive scalar between 0 & 1; a targeted decrease in % in p to be achieved every cycle. Default to 0.0.
+#' @param bias_abund a positive scalar larger than 0; bias in abundance estimates.
+#' For example, bias_abund = 0.5 (2) means that estimate abundance is on average
+#' half (twice) the true abundance. Default to 1, i.e. no bias.
+#' @param Ktrend a positive scalar; change in carrying capacity over the horizon
+#' for example, 0.5, means a halving of carrying capacity. Default to 1, i.e. no change.
+#' @param bias_byc a positive scalar larger than 0; bias in removals. For example,
+#' bias_byc = 2 means that twice the pbr is actually removed. Default to 1, i.e. no bias.
+#' @param catastrophe a positive scalar between 0 & 1;
+#' catastrophic mortality event as a fraction of current abundance.
+#' Happens randomly once during the management procedure.
+#' default to 0 i.e. no catastrophic mortality event.
+#' @param random logical; inherit the seed from pbrlist or rlalist for reproducibility purposes.
+#' Default to TRUE.
+#' @param averaging logical: should a time- and precision-aware averaged estimate
+#' of abundance be used in the PBR control rule? Default to FALSE.
+#'
+#' @importFrom dplyr last mutate
+#' @importFrom methods is
+#' @importFrom assertthat assert_that is.scalar is.number is.count
+#' @importFrom stats rlnorm qlnorm rbinom
+#'
+#' @return
+#' a list with two dataframes "management" and "depletion"
+#' "management" recapitulates the management objective of the procedure, with
+#' each line corresponding to a management cycle over the horizon.
+#' N_hat is the simulated survey estimate and depletion is the true depletion
+#' "depletion" recapitulates the whole population dynamics over the simulation horizon
+#'
+#' @export
+#'
+#' @seealso \code{\link{pellatomlinson_pbr}}, \code{\link{pellatomlinson_rla}},
+#' \link{forward_pbr}}, \link{pbr_nouveau}}, \code{\link{bycatch_obs}}
+#'
+#' @examples
+#' ### first, one needs to generate a list from a call to pellatomlinson_pbr
+#' set.seed(20221201)
+#' removals <- floor(runif(30, 1e3, 5e3))
+#' ### Rmax is set as for a cetacean species
+#' hp <- pellatomlinson_pbr(burnin = 50,
+#'                          vital_param = vital_param_pbr(
+#'                            K = 100000,
+#'                            depletion0 = 0.5,
+#'                            MNPL = 0.50,
+#'                            Rmax = 0.04
+#'                            ),
+#'                            catches = removals,
+#'                            seed_id = 20221201
+#'                            )
+#' ### now run the management procedure and evaluate the 1.7% rule
+#' ### and the default recovery factor
+#' test <- fixed_percentage(poplist = hp,
+#'                          distribution = "truncnorm",
+#'                          frequency = 6,
+#'                          horizon = 100,
+#'                          random = FALSE
+#'                          )
+#' test2 <- fixed_percentage(poplist = hp,
+#'                           distribution = "truncnorm",
+#'                           frequency = 6,
+#'                           horizon = 100,
+#'                           p_discount = 0.1,
+#'                           random = FALSE
+#'                           )
+#' test3 <- fixed_percentage(poplist = hp,
+#'                           distribution = "truncnorm",
+#'                           frequency = 6,
+#'                           horizon = 100,
+#'                           p_discount = 0.1,
+#'                           bias_byc = 2,
+#'                           random = FALSE
+#'                           )
+#' ### compare. Note that arg. random is FALSE by default so the same seed is used
+#' ### across the different tests. This can be seen below.
+#' plot(test$depletion$time, test$depletion$depletion,
+#'      type = 'l', las = 1, bty = 'n', xlab = "#years",
+#'      ylab = "depletion", ylim = c(0, 1)
+#'      )
+#' lines(test2$depletion$time, test2$depletion$depletion, col = "tomato", lty = 1)
+#' lines(test3$depletion$time, test3$depletion$depletion, col = "seagreen", lty = 1)
+#'
+#' @references
+#' Brandon, J. R.; Punt, A. E.; Moreno, P. & Reeves, R. R. (2017) Toward a Tier
+#' System Approach for Calculating Limits on Human-Caused Mortality of Marine Mammals.
+#' ICES Journal of Marine Science, 2017, 74, 877-887
+#'
+#' Wade, P. R. (1998) Calculating Limits To the Total Allowable Human-Caused
+#' Mortality of Cetaceans and Pinnipeds. Marine Mammal Science, 14:1-37.
+#' https://doi.org/10.1111/j.1748-7692.1998.tb00688.x
+fixed_percentage <- function(poplist,
+                             bycatch_variation_cv = 0.0,
+                             distribution = c("truncnorm", "lognormal", "gamma"),
+                             frequency = 6,
+                             horizon = 100,
+                             bias_abund = 1,
+                             Ktrend = 1,
+                             bias_byc = 1,
+                             catastrophe = 0,
+                             p = 0.017,
+                             p_discount = 0.0,
+                             random = FALSE,
+                             averaging = FALSE
+                             ) {
+  assert_that(is.list(poplist))
+  assert_that(is.number(bycatch_variation_cv), (bycatch_variation_cv >= 0))
+  assert_that(is.character(distribution),
+              distribution %in% c("truncnorm", "lognormal", "gamma")
+              )
+  assert_that(is.count(frequency))
+  assert_that(is.count(horizon))
+  assert_that(is.number(bias_abund), (bias_abund > 0))
+  assert_that(is.number(Ktrend), (Ktrend > 0))
+  assert_that(is.number(bias_byc), (bias_byc > 0))
+  assert_that(is.number(catastrophe), (catastrophe <= 1), (catastrophe >= 0))
+  assert_that(is.number(p), (p < 1), (p > 0))
+  assert_that(is.number(p_discount), (p_discount < 1), (p_discount >= 0))
+  assert_that(is.logical(random))
+  assert_that(is.logical(averaging))
+
+  ### for reproducibility
+  if(random) {
+    set.seed(sample.int(1e6, size = 1))
+  }
+  else {
+    set.seed(poplist$seed_id)
+  }
+
+  ### store population trajectory and bycatch history
+  catches <- 0
+  pop <- last(poplist$depletion)
+  ### when do new survey take place?
+  when <- seq(frequency, horizon - 1, frequency)
+  ### when to run the PBR?
+  dd <- data.frame(start = c(1, when + 1),
+                   end = c(when, horizon)
+                   ) %>%
+    mutate(catch_limit = NA,
+           N_hat = NA,
+           CV = NA,
+           delta_t = NA,
+           time_w = NA,
+           w = NA,
+           depletion = NA,
+           K = ceiling(seq(poplist$K, poplist$K * Ktrend, length.out = n())),
+           catastrophe = 0,
+           percent = p,
+           seed = sample.int(1e6, size = n())
+           )
+  ### catastrophic mortality event
+  if(catastrophe != 0) {
+    dd$catastrophe[sample.int(nrow(dd) - 1, size = 1, replace = FALSE)] <- 1
+  }
+
+  pt <- poplist
+  current_p <- p
+  ### pella-tomlinson dynamics
+  if(is(poplist, "pellatomlinson_pbr")) {
+    for(t in 1:nrow(dd)) {
+      ### compute catch limit
+      dd$N_hat[t] <- last(pt$survey$mean) * bias_abund
+      dd$CV[t] <- last(pt$survey$cv)
+      dd$depletion[t] <- last(pt$depletion)
+
+      # weights
+      dd$delta_t[1:t] <- (t - dd$end[1:t])
+      dd$time_w[1:t] <- 0.9^dd$delta_t[1:t]
+      dd$w[1:t] <- dd$time_w[1:t] / (dd$CV[1:t] * dd$CV[1:t])
+
+      ### averaging estimates?
+      if(averaging) {
+        Nbest <- dd$N_hat[t]
+        Nbest_cv <- dd$CV[t]
+      }
+      # time- and precision-aware averaging
+      else {
+        Nbest <- exp(sum(dd$time_w[1:t] * log(dd$N_hat[1:t]) / (dd$CV[1:t] * dd$CV[1:t])) / sum(dd$w[1:t]))
+        Nbest_cv <- sqrt(sum(dd$time_w[1:t] * dd$time_w[1:t] / (dd$CV[1:t] * dd$CV[1:t]))) / sum(dd$w[1:t])
+      }
+      # PBR
+      dd$catch_limit[t] <- floor(current_p * Nbest)
+      nt <- dd$end[t] - dd$start[t] + 1
+      ### removals
+      byc <- rep(dd$catch_limit[t], nt)
+      ## add variability
+      byc <- bycatch_obs(mu = byc,
+                         cv = bycatch_variation_cv,
+                         distribution = distribution
+                         )
+      ### run dynamics
+      pt <- pellatomlinson_pbr(
+        vital_param = vital_param_pbr(
+          depletion0 = last(pt$depletion),
+          # Shape parameter of the Pella-Tomlinson process
+          z = pt$z,
+          # carrying capacity can change over the simulation
+          K = dd$K[t],
+          # maximum growth rate
+          Rmax = pt$Rmax
+        ),
+        # series of catches
+        catches = byc * bias_byc,
+        # expected coefficient of variation for the survey estimates
+        CV = pt$CV,
+        # cv of environmental stochasticity
+        CV_env = pt$CV_env,
+        # correlated random walk
+        rho = pt$rho,
+        # years when surveys took place
+        scans = nt,
+        verbose = FALSE,
+        # for reproducibility purposes
+        seed_id = dd$seed[t]
+      )
+      ### is there a catastrophic mortality event after the survey?
+      if(dd$catastrophe[t] == 1) {
+        pt$N[nt + 1] <- rbinom(1, size = last(pt$N), prob = (1 - catastrophe))
+        pt$depletion[nt + 1] <- last(pt$N) / dd$K[t]
+      }
+      ### depletion level
+      pop <- c(pop, pt$depletion[-1])
+      ### catches
+      catches = c(catches, byc * bias_byc)
+      ### environmental objective
+      current_p <- current_p * (1 - p_discount)
+    }
+  }
+
+  if(is(poplist, "pellatomlinson_rla")) {
+    for(t in 1:nrow(dd)) {
+      ### compute catch limit
+      dd$N_hat[t] <- last(pt$survey$mean) * bias_abund
+      dd$CV[t] <- last(pt$survey$cv)
+      dd$depletion[t] <- last(pt$depletion)
+
+      # weights
+      dd$delta_t[1:t] <- (t - dd$end[1:t])
+      dd$time_w[1:t] <- 0.9^dd$delta_t[1:t]
+      dd$w[1:t] <- dd$time_w[1:t] / (dd$CV[1:t] * dd$CV[1:t])
+
+      ### averaging estimates?
+      if(averaging) {
+        Nbest <- dd$N_hat[t]
+        Nbest_cv <- dd$CV[t]
+      }
+      # time- and precision-aware averaging
+      else {
+        Nbest <- exp(sum(dd$time_w[1:t] * log(dd$N_hat[1:t]) / (dd$CV[1:t] * dd$CV[1:t])) / sum(dd$w[1:t]))
+        Nbest_cv <- sqrt(sum(dd$time_w[1:t] * dd$time_w[1:t] / (dd$CV[1:t] * dd$CV[1:t]))) / sum(dd$w[1:t])
+      }
+      # PBR
+      dd$catch_limit[t] <- floor(current_p * Nbest)
+      nt = dd$end[t] - dd$start[t] + 1
+      ### removals
+      byc = rep(dd$catch_limit[t], nt)
+      ## add variability
+      byc <- bycatch_obs(mu = byc,
+                         cv = bycatch_variation_cv,
+                         distribution = distribution
+                         )
+      ### run dynamics
+      pt <- pellatomlinson_rla(
+        vital_param = vital_param_rla(
+          Nf = pt$Nf,
+          Nm = pt$Nm,
+          z = pt$z,
+          K = dd$K[t],
+          L = pt$L,
+          eta = pt$eta,
+          phi = pt$phi,
+          m = pt$m,
+          b_max = pt$b_max,
+          b_K = pt$b_K
+        ),
+        catches = byc * bias_byc,
+        CV = pt$CV,
+        CV_env = pt$CV_env,
+        rho = pt$rho,
+        Allee = pt$Allee,
+        scans = nt,
+        verbose = FALSE,
+        seed_id = dd$seed[t]
+      )
+      ### is there a catastrophic mortality event after the survey?
+      ### females and males are equally affected
+      if(dd$catastrophe[t] == 1) {
+        pt$Nf <- sapply(pt$Nf, function(x) {rbinom(1, size = x, prob = (1 - catastrophe))})
+        pt$Nm <- sapply(pt$Nm, function(x) {rbinom(1, size = x, prob = (1 - catastrophe))})
+      }
+      ### depletion level
+      pop <- c(pop, pt$depletion[-1])
+      ### catches
+      catches = c(catches, byc * bias_byc)
+      ### environmental objective
+      current_p <- current_p * (1 - p_discount)
+    }
+  }
+
+  ### wrap-up
+  out <- list(management = dd %>%
+                mutate(bias_abund = bias_abund,
+                       bias_bycatch = bias_byc,
+                       seed = poplist$seed_id,
+                       averaging = ifelse(averaging, "time- and precision-aware", "none")
+                       ) %>%
+                dplyr::select(catch_limit,
+                              averaging,
+                              N_hat,
+                              depletion,
+                              K,
+                              percent,
+                              bias_abund,
+                              bias_bycatch,
+                              catastrophe,
+                              seed
+                              ),
+              depletion = data.frame(depletion = pop,
+                                     removals = catches
+                                     ) %>%
+                mutate(time = 0:horizon,
+                       percent = p,
+                       environmental_obj = p_discount,
+                       seed = poplist$seed_id
+                       )
+  )
+
+  return(out)
+}

R/rla-global.R

@@ -4,6 +4,7 @@ utils::globalVariables(
     "sigma", "mu", "Nmin", "x", "y", "l", "convergence",
     "time", "Abundance", "q10_abundance", "q90_abundance", "proj_min",
     "proj_max", "ymin", "ymax", "Birthrate", "Productivity", "survey_estimate",
-    "hp", "time", ".", "Depletion", "Removals", "removals", "recovery_factor"
+    "hp", "time", ".", "Depletion", "Removals", "removals", "recovery_factor",
+    "percent"
   )
 )

man/fixed_percentage.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fixed_percentage.R
+\name{fixed_percentage}
+\alias{fixed_percentage}
+\title{A function to carry out simulations according to the fixed percentage rule}
+\usage{
+fixed_percentage(
+  poplist,
+  bycatch_variation_cv = 0,
+  distribution = c("truncnorm", "lognormal", "gamma"),
+  frequency = 6,
+  horizon = 100,
+  bias_abund = 1,
+  Ktrend = 1,
+  bias_byc = 1,
+  catastrophe = 0,
+  p = 0.017,
+  p_discount = 0,
+  random = FALSE,
+  averaging = FALSE
+)
+}
+\arguments{
+\item{poplist}{a list; the output from one call to function \code{\link{pellatomlinson_pbr}} or \code{\link{pellatomlinson_rla}}}
+
+\item{bycatch_variation_cv}{a positive scalar; variability in bycatch}
+
+\item{distribution}{a character string; must be one of "truncnorm",
+"lognormal", or "gamma" to simulate catch numbers.}
+
+\item{frequency}{a count scalar; frequency of new surveys. Default to 6, i.e.
+a new survey every 6 years.}
+
+\item{horizon}{a count scalar; number of years to run the management procedure.
+Default to 100 years.}
+
+\item{bias_abund}{a positive scalar larger than 0; bias in abundance estimates.
+For example, bias_abund = 0.5 (2) means that estimate abundance is on average
+half (twice) the true abundance. Default to 1, i.e. no bias.}
+
+\item{Ktrend}{a positive scalar; change in carrying capacity over the horizon
+for example, 0.5, means a halving of carrying capacity. Default to 1, i.e. no change.}
+
+\item{bias_byc}{a positive scalar larger than 0; bias in removals. For example,
+bias_byc = 2 means that twice the pbr is actually removed. Default to 1, i.e. no bias.}
+
+\item{catastrophe}{a positive scalar between 0 & 1;
+catastrophic mortality event as a fraction of current abundance.
+Happens randomly once during the management procedure.
+default to 0 i.e. no catastrophic mortality event.}
+
+\item{p}{a positive scalar between 0 & 1; the percentage to evaluate. Default to 0.017.}
+
+\item{p_discount}{a positive scalar between 0 & 1; a targeted decrease in \% in p to be achieved every cycle. Default to 0.0.}
+
+\item{random}{logical; inherit the seed from pbrlist or rlalist for reproducibility purposes.
+Default to TRUE.}
+
+\item{averaging}{logical: should a time- and precision-aware averaged estimate
+of abundance be used in the PBR control rule? Default to FALSE.}
+}
+\value{
+a list with two dataframes "management" and "depletion"
+"management" recapitulates the management objective of the procedure, with
+each line corresponding to a management cycle over the horizon.
+N_hat is the simulated survey estimate and depletion is the true depletion
+"depletion" recapitulates the whole population dynamics over the simulation horizon
+}
+\description{
+A function to carry out simulations according to the fixed percentage rule
+}
+\examples{
+### first, one needs to generate a list from a call to pellatomlinson_pbr
+set.seed(20221201)
+removals <- floor(runif(30, 1e3, 5e3))
+### Rmax is set as for a cetacean species
+hp <- pellatomlinson_pbr(burnin = 50,
+                         vital_param = vital_param_pbr(
+                           K = 100000,
+                           depletion0 = 0.5,
+                           MNPL = 0.50,
+                           Rmax = 0.04
+                           ),
+                           catches = removals,
+                           seed_id = 20221201
+                           )
+### now run the management procedure and evaluate the 1.7\% rule
+### and the default recovery factor
+test <- fixed_percentage(poplist = hp,
+                         distribution = "truncnorm",
+                         frequency = 6,
+                         horizon = 100,
+                         random = FALSE
+                         )
+test2 <- fixed_percentage(poplist = hp,
+                          distribution = "truncnorm",
+                          frequency = 6,
+                          horizon = 100,
+                          p_discount = 0.1,
+                          random = FALSE
+                          )
+test3 <- fixed_percentage(poplist = hp,
+                          distribution = "truncnorm",
+                          frequency = 6,
+                          horizon = 100,
+                          p_discount = 0.1,
+                          bias_byc = 2,
+                          random = FALSE
+                          )
+### compare. Note that arg. random is FALSE by default so the same seed is used
+### across the different tests. This can be seen below.
+plot(test$depletion$time, test$depletion$depletion,
+     type = 'l', las = 1, bty = 'n', xlab = "#years",
+     ylab = "depletion", ylim = c(0, 1)
+     )
+lines(test2$depletion$time, test2$depletion$depletion, col = "tomato", lty = 1)
+lines(test3$depletion$time, test3$depletion$depletion, col = "seagreen", lty = 1)
+
+}
+\references{
+Brandon, J. R.; Punt, A. E.; Moreno, P. & Reeves, R. R. (2017) Toward a Tier
+System Approach for Calculating Limits on Human-Caused Mortality of Marine Mammals.
+ICES Journal of Marine Science, 2017, 74, 877-887
+
+Wade, P. R. (1998) Calculating Limits To the Total Allowable Human-Caused
+Mortality of Cetaceans and Pinnipeds. Marine Mammal Science, 14:1-37.
+https://doi.org/10.1111/j.1748-7692.1998.tb00688.x
+}
+\seealso{
+\code{\link{pellatomlinson_pbr}}, \code{\link{pellatomlinson_rla}},
+\link{forward_pbr}}, \link{pbr_nouveau}}, \code{\link{bycatch_obs}}
+}

man/rla_nouveau.Rd

@@ -197,6 +197,13 @@ Uses a likelihood from a Stochastic Surplus Production Model (Ouzoulias 2022)
 Ouzoulias, F. (2022) Bayesian exploration of Surplus Production Models
 for cetaceans by-catch management, using strandings and abundance data.
 MSc Thesis. L'Institut Agro Rennes-Angers.
+
+Bordet, C. and Rivest, L.-P. (2014) A stochastic Pella Tomlinson model and its
+maximum sustainable yield. Journal of Theoretical Biology, 360:46–53.
+
+Bousquet, N., Duchesne, T., and Rivest, L. P. (2008) Redefining the maximum
+sustainable yield for the Schaefer population model including multiplicative
+environmental noise. Journal of theoretical biology, 254 1:65–75.
 }
 \seealso{
 \code{\link{pellatomlinson_rla}}, \code{\link{bycatch_obs}},

tests/testthat/test-fixed_percentage.R

+test_that("multiplication works", {
+  expect_equal(2 * 2, 4)
+})