scenarios

library(tidyverse)
library(knitr)
library(RLA)

Generating scenarios

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.

Loading the life-history parameters of the harbour porpoises

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 1

Simulations

We 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   25

Some 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’.