library(tidyverse)
library(knitr)
library(RLA)The RLA is a management procedure in which you need to use simulations to assess the objective of interest. The first step is to define your simulation scenarios. This vignette describes how to do that with the example of Harbour porpoise (Phocoena phocoena) in the North Sea. We’ll be simulating data with the function pellatomlinson_rla and our aims here to obtain a specified level of current depletion of the simulated populations (given the life-history parameters) in order to run the function forward_rla on these simulations to achieve tuning (see the vignette on ‘simulations’). To obtain simulations with the desired level of current depletion, we will use a sort of reject-accept procedure. This procedure is wasteful but needs only to be done once as results can be stored with the seeds you’re using.
library("tidyverse")
data("north_sea_hp")
str(north_sea_hp$life_history)
#> List of 6
#> $ K : num 5e+05
#> $ L : num 22
#> $ phi : num [1:23] 0.85 0.87 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.91 ...
#> $ maturity: num [1:23] 0 0.000335 0.002473 0.017986 0.119203 ...
#> $ eta : num [1:23] 2 2 1 1 1 1 1 1 1 1 ...
#> $ sexratio: num [1:2] 1 1We start by specifying a large number of simulations (well, not so large as this is just for illustrative purposes. We’ll use \(500\), which is too small.) out of which we will keep only a small fraction that conforms to our desiderata. The latter are: - Maximum Net Productivity of either \(2, ~4, ~6\)%; - Maximum Net Productivity Level of either \(50, ~60, ~70\)% of \(K\), the carrying capacity; and - current depletion level of the population of \(50, ~70, ~90\)% of \(K\).
n_sim <- 5e2
seeds <- sample.int(1e6, size = n_sim, replace = FALSE)
MNPL <- sample(c(0.5, 0.6, 0.7), n_sim, replace = TRUE)
cv_byc <- round(runif(n_sim, 0.05, 0.50), 3)
cv_env <- round(runif(n_sim, 0.00, 0.20), 3)
cv_obs <- round(runif(n_sim, 0.15, 0.30), 3)
MNP <- sample(c(1.02, 1.04, 1.06), size = n_sim, replace = TRUE)
r <- round(runif(n_sim, 0.001, 0.050), 3)The input \(r\) is used to simulate a rate of unmanaged bycatch with the function bycatch_rw. We set verbose to FALSE to avoid message printing.
set.seed(20210302)
rla <- lapply(1:n_sim,
function(i) {
pellatomlinson_rla(MNPL = MNPL[i],
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 = MNP[i],
# series of catches
catches = bycatch_rw(n = 51,
K = north_sea_hp$life_history$K,
rate = r[i],
cv = cv_byc[i],
seed_id = seeds[i]
),
CV = cv_obs[i],
CV_env = cv_env[i],
# scans survey
scans = c(29, 40, 51),
verbose = FALSE,
# for reproducibility purposes
seed_id = seeds[i]
)
}
)The next step is to identify the simulations with the desired current level of depletion:
depletion <- data.frame(current_depletion = do.call('c',
lapply(rla,
function(l) {
l$depletion[length(l$depletion)]
})
)
) %>%
mutate(MNP = MNP,
MNPL = MNPL,
level = cut(current_depletion, breaks = c(seq(0.4, 0.8, 0.2), 1.0)),
sim = 1:n() # indicator variable
) %>%
filter(!is.na(level)) # remove heavily depleted populations
with(depletion, table(MNPL, MNP, level))
#> , , level = (0.4,0.6]
#>
#> MNP
#> MNPL 1.02 1.04 1.06
#> 0.5 7 6 8
#> 0.6 6 2 3
#> 0.7 3 6 1
#>
#> , , level = (0.6,0.8]
#>
#> MNP
#> MNPL 1.02 1.04 1.06
#> 0.5 7 5 5
#> 0.6 10 2 3
#> 0.7 8 4 6
#>
#> , , level = (0.8,1]
#>
#> MNP
#> MNPL 1.02 1.04 1.06
#> 0.5 6 12 16
#> 0.6 7 14 12
#> 0.7 11 31 25Some scenarios are more difficult to obtain. For example having a population depleted to \(\approx 50\)% of \(\mathrm{K}\) with an \(\mathrm{MNP}\) of \(6\)% and the \(\mathrm{MNPL}\) occurring at \(\approx 70\)% of \(\mathrm{K}\). We can change that by increasing the value of inout r for the bycatch rate in the function bycatch_rw. We can save \(1\) of each of the \(3 \times 3 \times 3 = 18\) scenarios we’re interested in.
set.seed(19790428)
n_sim <- 1e2
keep <- depletion %>%
group_by(MNP, level, MNPL) %>%
slice_sample(n = n_sim) %>%
mutate(cv_byc = round(cv_byc[sim], 3),
cv_env = round(cv_env[sim], 3),
cv_obs = round(cv_obs[sim], 3),
rate = round(r[sim], 3),
seed = seeds[sim]
) %>%
dplyr::select(-sim) %>%
as.data.frame()
keep %>%
head(5) %>%
kable()| current_depletion | MNP | MNPL | level | cv_byc | cv_env | cv_obs | rate | seed |
|---|---|---|---|---|---|---|---|---|
| 0.5820464 | 1.02 | 0.5 | (0.4,0.6] | 0.204 | 0.057 | 0.215 | 0.012 | 826165 |
| 0.4865971 | 1.02 | 0.5 | (0.4,0.6] | 0.114 | 0.066 | 0.223 | 0.019 | 670250 |
| 0.5731640 | 1.02 | 0.5 | (0.4,0.6] | 0.173 | 0.122 | 0.160 | 0.013 | 425144 |
| 0.4209884 | 1.02 | 0.5 | (0.4,0.6] | 0.303 | 0.079 | 0.271 | 0.020 | 686324 |
| 0.5122780 | 1.02 | 0.5 | (0.4,0.6] | 0.254 | 0.004 | 0.279 | 0.018 | 104151 |
These scenarios can be used now for tuning the RLA. They are in fact stored as data in a table scenarios_north_sea_hp. See the vignette on ‘simulations’.