Commit 4958c54e authored by Mathieu Genu's avatar Mathieu Genu
Browse files

- modify `vital_param_args()` into vital_param_pbr() for **PBR** and add...

- modify `vital_param_args()` into vital_param_pbr() for **PBR** and add `vital_param_rla()` for **RLA** and **PBR_nouveau**

- all test modified and run correctly, examples of all function adapted
parent 5de92156
Pipeline #31350 failed with stage
......@@ -24,7 +24,8 @@ export(stableage)
export(standata)
export(summary_plot)
export(time2CO)
export(vital_param_args)
export(vital_param_pbr)
export(vital_param_rla)
export(zero_mortality)
import(dplyr)
import(ggplot2)
......
......@@ -54,14 +54,14 @@
#' @examples
#' ### first, one needs to generate a list from a call to pellatomlinson_pbr
#' ### Rmax is set as for a cetacean species
#' hp <- pellatomlinson_pbr(depletion0 = 0.3,
#' MNPL = 0.6,
#' Rmax = 0.04,
#' catches = 0,
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#'
#' hp <- pellatomlinson_pbr(
#' vital_param = vital_param_pbr(depletion0 = 0.3,
#' MNPL = 0.6,
#' Rmax = 0.04),
#' catches = 0,
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' ### now run the management procedure and evaluate the 20% quantile
#' ### and the default recovery factor
#' test <- forward_pbr(pbrlist = hp,
......@@ -70,8 +70,7 @@
#' horizon = 100,
#' q = 0.2,
#' random = FALSE
#' )
#'
#' )
#' ### includes a catastrophic event whereby 10% of the population is wiped out
#' test2 <- forward_pbr(pbrlist = hp,
#' distribution = "truncnorm",
......@@ -81,7 +80,7 @@
#' F_r = 0.5,
#' q = 0.2,
#' random = FALSE
#' )
#' )
#' ### includes a halving of the carrying capacity
#' test3 <- forward_pbr(pbrlist = hp,
#' distribution = "truncnorm",
......@@ -90,7 +89,7 @@
#' Ktrend = 0.5,
#' q = 0.2,
#' random = FALSE
#' )
#' )
#' ### includes a catastrophic event whereby 10% of the population is wiped out
#' ### and a halving of the carrying capacity
#' test4 <- forward_pbr(pbrlist = hp,
......@@ -101,17 +100,16 @@
#' catastrophe = 0.1,
#' q = 0.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 = "midnightblue", lty = 1)
#' lines(test4$depletion$time, test4$depletion$depletion, col = "seagreen", lty = 1)
#'
#' @references
#' Wade, P. R. (1998) Calculating Limits To the Total Allowable Human-Caused
#' Mortality of Cetaceans and Pinnipeds. Marine Mammal Science, 14:1-37.
......@@ -199,8 +197,8 @@ forward_pbr <- function(pbrlist,
distribution = distribution
)
### run dynamics
pt <- pellatomlinson_pbr(depletion0 = last(pt$depletion),
vital_param = vital_param_args(
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
......
......@@ -76,66 +76,68 @@
#'
#' @examples
#' \dontrun{
#' ## see also the vignette 'simulations'
#' library("tidyverse")
#' ### load data for the harbour porpoise in the North Sea
#' data("north_sea_hp")
#' data("scenarios_north_sea_hp")
#' ### use the first scenario
#' dyn <- pellatomlinson_rla(MNPL = scenarios_north_sea_hp$MNPL[1],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[1],
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[1],
#' cv = scenarios_north_sea_hp$cv_byc[1],
#' seed_id = scenarios_north_sea_hp$seed[1]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[1],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[1],
#' # scans survey
#' scans = c(29, 40, 51),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[1]
#' )
#' library(rstan)
#' options(mc.cores = parallel::detectCores())
#' rstan_options(auto_write = TRUE)
#'
#' ## load Stan models
#' data(rlastan_models)
#' # use uniform priors: have a look at the Stan code
#' cat(rlastan_models$uniform)
#' # compile model
#' rlastan <- rstan::stan_model(model_code = rlastan_models$uniform,
#' model_name = "RLA"
#' )
#' ## call the rla management procedure
#' rla <- forward_rla(rlalist = dyn,
#' rlastan = rlastan,
#' q = 0.5,
#' distribution = "truncnorm",
#' # upper bound for the prior on r, the population growth rate
#' upper = 0.1
#' )
#' ## have a look at the output
#' # management outcome
#' rla$management
#' # population trajectory
#' rla$depletion %>%
#' ggplot() +
#' geom_line(aes(x = 551 + time, y = depletion), color = "midnightblue") +
#' theme_bw()
#' # data passed to Stan in the final call
#' rla$standata
#' ## see also the vignette 'simulations'
#' library("tidyverse")
#' ### load data for the harbour porpoise in the North Sea
#' data("north_sea_hp")
#' data("scenarios_north_sea_hp")
#' ### use the first scenario
#' dyn <- pellatomlinson_rla(
#' vital_param = vital_param_rla(
#' MNPL = scenarios_north_sea_hp$MNPL[1],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[1]
#' ),
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[1],
#' cv = scenarios_north_sea_hp$cv_byc[1],
#' seed_id = scenarios_north_sea_hp$seed[1]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[1],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[1],
#' # scans survey
#' scans = c(29, 40, 51),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[1]
#' )
#' library(rstan)
#' options(mc.cores = parallel::detectCores())
#' rstan_options(auto_write = TRUE)
#' ## load Stan models
#' data(rlastan_models)
#' # use uniform priors: have a look at the Stan code
#' cat(rlastan_models$uniform)
#' # compile model
#' rlastan <- rstan::stan_model(model_code = rlastan_models$uniform,
#' model_name = "RLA"
#' )
#' ## call the rla management procedure
#' rla <- forward_rla(rlalist = dyn,
#' rlastan = rlastan,
#' q = 0.5,
#' distribution = "truncnorm",
#' # upper bound for the prior on r, the population growth rate
#' upper = 0.1
#' )
#' ## have a look at the output
#' # management outcome
#' rla$management
#' # population trajectory
#' rla$depletion %>%
#' ggplot() +
#' geom_line(aes(x = 551 + time, y = depletion), color = "midnightblue") +
#' theme_bw()
#' # data passed to Stan in the final call
#' rla$standata
#' }
forward_rla <- function(rlalist,
rlastan,
......@@ -264,27 +266,30 @@ forward_rla <- function(rlalist,
distribution = distribution
)
### run dynamics
pt <- pellatomlinson_rla(Nf = pt$Nf,
Nm = pt$Nm,
z = pt$z,
K = dd$K[t],
L = pt$L,
eta = pt$eta,
sexbias = pt$sexbias,
phi = pt$phi,
m = pt$m,
sexratio = pt$sexratio,
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]
)
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,
sexbias = pt$sexbias,
phi = pt$phi,
m = pt$m,
sexratio = pt$sexratio,
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) {
......
......@@ -51,96 +51,98 @@
#'
#' @examples
#' \dontrun{
#' ## see also the vignette 'simulations'
#' ### load data for the harbour porpoise in the North Sea
#' data("north_sea_hp")
#' data("scenarios_north_sea_hp")
#' sim <- 1
#' dyn <- pellatomlinson_rla(MNPL = scenarios_north_sea_hp$MNPL[sim],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[sim],
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[sim],
#' cv = scenarios_north_sea_hp$cv_byc[sim],
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[sim],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[sim],
#' # scans survey
#' scans = c(29, 40, 51),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' )
#'
#' ## call the pbr_nouveau management procedure
#' res <- pbr_nouveau(rlalist = dyn,
#' distribution = "truncnorm",
#' q = 0.5,
#' F_r = 0.5
#' )
#' ## have a look at the output
#' # management outcome
#' res$management
#' # population trajectory
#' plot(res$depletion$time,
#' res$depletion$depletion,
#' las = 1, type = 'l', xlab = "# years", ylab = "depletion",
#' bty = 'n'
#' )
#'
#' ### compare with depensation: assume that depensation occurs at 50% of male depletion
#' dyn <- pellatomlinson_rla(MNPL = scenarios_north_sea_hp$MNPL[sim],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[sim],
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[sim],
#' cv = scenarios_north_sea_hp$cv_byc[sim],
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[sim],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[sim],
#' # scans survey
#' scans = c(29, 40, 51),
#' # depensation: only 0.95 of productivity when male depletion is 30%
#' Allee = list(meetup = 0.95, male_depletion = 0.3),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' )
#'
#' ## call the pbr_nouveau management procedure
#' res2 <- pbr_nouveau(rlalist = dyn,
#' distribution = "truncnorm",
#' q = 0.5,
#' F_r = 0.5
#' )
#' # population trajectory: compare with no depensation
#' plot(res$depletion$time,
#' res$depletion$depletion,,
#' bty = 'n', las = 1, type = 'l', lty = 2,
#' xlab = "# years", ylab = "depletion", ylim = c(0.5, 0.65)
#' )
#' lines(res2$depletion$time,
#' res2$depletion$depletion,
#' color = "midnightblue"
#' )
#' ## see also the vignette 'simulations'
#' ### load data for the harbour porpoise in the North Sea
#' data("north_sea_hp")
#' data("scenarios_north_sea_hp")
#' sim <- 1
#' dyn <- pellatomlinson_rla(
#' vital_param = vital_param_rla(MNPL = scenarios_north_sea_hp$MNPL[sim],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[sim]
#' ),
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[sim],
#' cv = scenarios_north_sea_hp$cv_byc[sim],
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[sim],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[sim],
#' # scans survey
#' scans = c(29, 40, 51),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' )
#' ## call the pbr_nouveau management procedure
#' res <- pbr_nouveau(rlalist = dyn,
#' distribution = "truncnorm",
#' q = 0.5,
#' F_r = 0.5
#' )
#' ## have a look at the output
#' # management outcome
#' res$management
#' # population trajectory
#' plot(res$depletion$time,
#' res$depletion$depletion,
#' las = 1, type = 'l', xlab = "# years", ylab = "depletion",
#' bty = 'n'
#' )
#' ### compare with depensation: assume that depensation occurs at 50% of male depletion
#' dyn <- pellatomlinson_rla(
#' vital_param = vital_param_rla(
#' MNPL = scenarios_north_sea_hp$MNPL[sim],
#' K = north_sea_hp$life_history$K,
#' L = north_sea_hp$life_history$L,
#' eta = north_sea_hp$life_history$eta,
#' phi = north_sea_hp$life_history$phi,
#' m = north_sea_hp$life_history$maturity,
#' MNP = scenarios_north_sea_hp$MNP[sim]
#' ),
#' # series of catches
#' catches = bycatch_rw(n = 51,
#' K = north_sea_hp$life_history$K,
#' rate = scenarios_north_sea_hp$rate[sim],
#' cv = scenarios_north_sea_hp$cv_byc[sim],
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' ),
#' # expected coefficient of variation for the survey estimates
#' CV = scenarios_north_sea_hp$cv_obs[sim],
#' # CV of environmental stochasticity on birth rate
#' CV_env = scenarios_north_sea_hp$cv_env[sim],
#' # scans survey
#' scans = c(29, 40, 51),
#' # depensation: only 0.95 of productivity when male depletion is 30%
#' Allee = list(meetup = 0.95, male_depletion = 0.3),
#' verbose = TRUE,
#' everything = TRUE,
#' seed_id = scenarios_north_sea_hp$seed[sim]
#' )
#' ## call the pbr_nouveau management procedure
#' res2 <- pbr_nouveau(rlalist = dyn,
#' distribution = "truncnorm",
#' q = 0.5,
#' F_r = 0.5
#' )
#' # population trajectory: compare with no depensation
#' plot(res$depletion$time,
#' res$depletion$depletion,,
#' bty = 'n', las = 1, type = 'l', lty = 2,
#' xlab = "# years", ylab = "depletion", ylim = c(0.5, 0.65)
#' )
#' lines(res2$depletion$time,
#' res2$depletion$depletion,
#' col = "midnightblue"
#' )
#' }
pbr_nouveau <- function(rlalist,
bycatch_variation_cv = 0.0,
......@@ -224,25 +226,28 @@ pbr_nouveau <- function(rlalist,
distribution = distribution
)
### run dynamics
pt <- pellatomlinson_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]
)
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) {
......
......@@ -3,15 +3,16 @@
#'
#' @param burnin a count scalar. The number of time steps to run the Pella-Tomlison
#' process without anthropogenic removals. Default to 0.
#' @param depletion0 positive scale; depletion at the beginning of the
#' population trajectory.
#' @param vital_param
#' argument have to be included in the function `vital_param_pbr()` (ex : `vital_param = vital_param_pbr(MNPL = 0.6)`)
#' - **MNPL** a proportion; the maximum net productivity level as fraction of
#' carrying capacity. Default to 0.6.
#' - **z** a positive scalar (internal). Default to NULL.
#' - **K** a count scalar; carrying capacity. Default to 2.5E5.
#' - **Rmax** a positive scalar; maximum growth rate,
#' set traditionally to 0.04 for cetaceans and 0.12 for seals.
#' - **depletion0** positive scale; depletion at the beginning of the
#' population trajectory.
#' @param catches a count vector; time-series of anthropogenic catches of animals.
#' Default to 0.
#' @param CV a positive scalar; expected coefficient of variation for the survey estimates.
......@@ -61,48 +62,43 @@
#' set.seed(123)
#' ## series of catches
#' hp_catches <- floor(runif(30, 1e3, 5e3))
#'
#' system.time(
#' hp <- pellatomlinson_pbr(burnin = 5e2,
#' depletion0 = 0.3,
#' vital_param = list(Rmax = 0.04),
#' vital_param = vital_param_pbr(
#' depletion0 = 0.3,
#' Rmax = 0.04),
#' catches = hp_catches,
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' )
#' )
#'
#' plot(hp$N, type = 'l', las = 1, bty = 'n', xlab = "#years", ylab = "Abundance")
#'
#' ### an example to simulate directly (no burnin)
#' set.seed(123)
#' ## series of catches
#' hp_catches <- floor(runif(30, 1e3, 5e3))
#'
#' system.time(
#' hp <- pellatomlinson_pbr(depletion0 = 0.3,
#' vital_param = list(Rmax = 0.04),
#' catches = hp_catches,
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' hp <- pellatomlinson_pbr(
#' vital_param = vital_param_pbr(depletion0 = 0.3,
#' Rmax = 0.04),
#' catches = hp_catches,
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' )
#'
#' plot(hp$N, type = 'l', las = 1, bty = 'n', xlab = "#years", ylab = "Abundance")
#'
#' ### an example with no simulation (for use with forward_rla)
#' system.time(
#' hp <- pellatomlinson_pbr(depletion0 = 0.3,
#' vital_param = list(Rmax = 0.04),
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' hp <- pellatomlinson_pbr(
#' vital_param = vital_param_pbr(depletion0 = 0.3,
#' Rmax = 0.04),
#' CV_env = 0.01,
#' seed_id = 20210219
#' )
#' )
#'
pellatomlinson_pbr <- function(
vital_param = vital_param_args(),
vital_param = vital_param_pbr(),
burnin = 0,
depletion0 = NULL,
catches = NULL,
CV = 0.3,
CV_env = 0,
......@@ -114,10 +110,11 @@ pellatomlinson_pbr <- function(
) {
# define args = list$args
MNPL <- vital_param$MNPL
z <- vital_param$z
K <- vital_param$K
Rmax <- vital_param$Rmax
depletion0 <- vital_param$depletion0
MNPL <- vital_param$MNPL
z <- vital_param$z
K <- vital_param$K
Rmax <- vital_param$Rmax
### sanity checks
assert_that(is.count(burnin) || burnin == 0)
......@@ -282,15 +279,28 @@ pellatomlinson_pbr <- function(
#' Default vital_param
#' @noRd
vital_param_args <- function(
MNPL = 0.6,z = NULL,K = 2.5E5,Rmax = NULL
#' Default vital_param for PBR
#'
#' @param MNPL Maximum Net Productivity Level
#' @param z shape parameter of pella-tomlinson growth model
#' @param K Carrying capacity
#' @param Rmax Productivity rate
#'
#' @seealso [`RLA::pellatomlinson_pbr()`]
#'
#' @export
vital_param_pbr <- function(
MNPL = 0.6,
z = NULL,
K = 2.5E5,
Rmax = NULL,