library("tidyverse")
#> -- Attaching packages --------------------------------------- tidyverse 1.3.0 --
#> v ggplot2 3.3.3 v purrr 0.3.4
#> v tibble 3.1.0 v dplyr 1.0.5
#> v tidyr 1.1.3 v stringr 1.4.0
#> v readr 1.4.0 v forcats 0.5.1
#> -- Conflicts ------------------------------------------ tidyverse_conflicts() --
#> x dplyr::filter() masks stats::filter()
#> x dplyr::lag() masks stats::lag()
library("knitr")
library("RLA")This vignette is to showcase some experiments with the RLA package, using the data on the harbour porpoises (Phocoena phocoena) in the North Sea. The data is included in the package.
Potential Biological Removal (Wade 1998) is a management procedure to set limits to anthropogenic removals. The operation model behind PBR is a deterministic Pella-Tomlinson process, which is emulated in function pellatomlinson_pbr. However, if there are data on survival and fecundity, but no data on removals, it is still possible to use the RLA operating model and tune it to a conservation objective using a PBR-like approach. This what the function pbr_nouveau does.
To use it, we need to load scenarios to pass to pella_tomlinson_rla.
data("north_sea_hp")
data("scenarios_north_sea_hp")
scenarios_north_sea_hp %>%
head(5) %>%
kable()| current_depletion | MNP | MNPL | level | cv_obs | cv_env | cv_byc | rate | seed |
|---|---|---|---|---|---|---|---|---|
| 0.5417200 | 1.02 | 0.5 | (0.4,0.6] | 0.231 | 0.052 | 0.195 | 0.014 | 271317 |
| 0.4005977 | 1.02 | 0.5 | (0.4,0.6] | 0.239 | 0.065 | 0.388 | 0.021 | 921891 |
| 0.4947520 | 1.02 | 0.5 | (0.4,0.6] | 0.193 | 0.110 | 0.430 | 0.012 | 803073 |
| 0.5601648 | 1.02 | 0.5 | (0.4,0.6] | 0.298 | 0.155 | 0.428 | 0.008 | 737284 |
| 0.4948804 | 1.02 | 0.5 | (0.4,0.6] | 0.226 | 0.070 | 0.179 | 0.014 | 191156 |
We illustrate the function the first \(100\) scenarios of the table scenario_north_sea_hp.
scenario1 <- lapply(1:100, function(i) {
pellatomlinson_rla(MNPL = scenarios_north_sea_hp$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 = scenarios_north_sea_hp$MNP[i],
# series of catches
catches = bycatch_rw(n = 51,
K = north_sea_hp$life_history$K,
rate = scenarios_north_sea_hp$rate[i],
cv = scenarios_north_sea_hp$cv_byc[i],
seed_id = scenarios_north_sea_hp$seed[i]
),
# expected coefficient of variation for the survey estimates
CV = scenarios_north_sea_hp$cv_obs[i],
# CV of environmental stochasticity on birth rate
CV_env = scenarios_north_sea_hp$cv_env[i],
# scans survey
scans = c(29, 40, 51),
verbose = FALSE,
everything = TRUE,
seed_id = scenarios_north_sea_hp$seed[i]
)
})These scenarios can be visualized:
trajectory <- do.call('rbind',
lapply(1:length(scenario1), function(i) {
data.frame(time = 0:(length(scenario1[[i]]$depletion) - 1),
depletion = scenario1[[i]]$depletion,
seed = rep(scenario1[[i]]$seed_id, length(scenario1[[i]]$depletion))
) %>%
left_join(scenarios_north_sea_hp,
by = "seed"
)
})
)
trajectory %>%
ggplot(aes(x = time, y = depletion, group = seed)) +
geom_line(color = "midnightblue", alpha = 0.1) +
scale_x_continuous(name = "# years",
breaks = seq(0, 600, 50)
) +
scale_y_continuous(name = "Depletion (as % of K)",
breaks = seq(0.0, 2.0, 0.2),
labels = 100 * seq(0.0, 2.0, 0.2)
) +
facet_grid(MNP ~ MNPL) +
theme_bw()Each of these simulation can be fed to pbr_nouveau to tune the procedure. This is different from the traditional PBR as information on survival and reproduction that are specific to a population of interest. In this was, even if information on bycatch may be missing, any information on local vital rates may be used.
We’ll start by evaluating the \(20\)% quantile of the best available abundance estimate (simulated as part of the function pellatomlinson_rla):
q20 <- lapply(scenario1,
pbr_nouveau,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
F_r = 0.5,
q = 0.2
)The different population trajectories associated with this choice of a quantile can easily be visualized. Note that in pbr_nouveau, the \(\mathrm{MNP}\) input of the operating model is recycled as the \(0.5 \times \mathrm{R_{max}}\) input if \(\mathrm{R_{max}}\) is not specified by the user. See Wade (1998) for the PBR formula for the removal limit, and the relationship between \(\mathrm{MNP}\) and \(\mathrm{R_{max}}\).
trajectory <- do.call('rbind',
lapply(1:length(q20), function(i) { q20[[i]]$depletion })
)
ggplot() +
geom_hline(yintercept = 0.8, linetype = "dotted", color = "tomato", size = 1) +
geom_line(data = trajectory,
aes(x = time, y = depletion, group = factor(seed * quantile), color = quantile), alpha = 0.1
) +
geom_line(data = trajectory %>%
group_by(time, quantile) %>%
summarize(mean = mean(depletion),
lower = as.numeric(quantile(depletion, probs = 0.025)),
upper = as.numeric(quantile(depletion, probs = 0.975))
),
aes(x = time, y = mean, group = quantile, color = quantile), size = 2
) +
scale_y_continuous(name = "Depletion (as % of K)",
breaks = seq(0.0, 1.0, 0.1),
labels = 100 * seq(0.0, 1.0, 0.1)
) +
ggtitle("Rmax = 4%, MNPL = 50%\nQuantile = 20%") +
guides(color = "none") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.The argument \(\mathrm{R_{max}}\) can be specified directly. The point is to test several quantiles to tune the procedure:
quantile2evaluate <- seq(0.1, 0.9, 0.1)
tuning <- lapply(quantile2evaluate,
function(q) {
lapply(scenario1,
pbr_nouveau,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
F_r = 0.5,
Rmax = 0.04,
q = q
)
}
)
trajectory <- do.call('rbind',
lapply(1:length(tuning), function(l) {
do.call('rbind',
lapply(1:length(tuning[[l]]), function(i) { tuning[[l]][[i]]$depletion })
)
})
)
ggplot() +
geom_hline(yintercept = 0.8, linetype = "dotted", color = "tomato", size = 1) +
geom_line(data = trajectory %>%
group_by(time, quantile) %>%
summarize(mean = mean(depletion),
lower = as.numeric(quantile(depletion, probs = 0.025)),
upper = as.numeric(quantile(depletion, probs = 0.975))
),
aes(x = time, y = mean, group = quantile, color = quantile), size = 1
) +
scale_y_continuous(name = "Depletion (as % of K)",
breaks = seq(0.0, 1.0, 0.1),
labels = 100 * seq(0.0, 1.0, 0.1)
) +
ggtitle("Rmax = 4%, MNPL = 50%") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.In this case, none of the quantile allows the population to recover to at least \(80\)% of \(\mathrm{K}\) over the next \(100\) years. This is a simple scenario, which can be elaborated upon: we’ll simulate again a catastrophic event that will happen once in those \(100\) years and that will kill immediately \(10\)% of the population.
tuning <- lapply(quantile2evaluate,
function(q) {
lapply(scenario1,
pbr_nouveau,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
catastrophe = 0.1,
F_r = 0.5,
Rmax = 0.04,
q = q
)
}
)
trajectory <- do.call('rbind',
lapply(1:length(tuning), function(l) {
do.call('rbind',
lapply(1:length(tuning[[l]]), function(i) { tuning[[l]][[i]]$depletion })
)
})
)
ggplot() +
geom_hline(yintercept = 0.8, linetype = "dotted", color = "tomato", size = 1) +
geom_line(data = trajectory %>%
group_by(time, quantile) %>%
summarize(mean = mean(depletion),
lower = as.numeric(quantile(depletion, probs = 0.025)),
upper = as.numeric(quantile(depletion, probs = 0.975))
),
aes(x = time, y = mean, group = quantile, color = quantile), size = 1
) +
scale_y_continuous(name = "Depletion (as % of K)",
breaks = seq(0.0, 1.0, 0.1),
labels = 100 * seq(0.0, 1.0, 0.1)
) +
ggtitle("Rmax = 4%, MNPL = 50%\nCatastrophic event wiping\nout 10% of animals") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.And finally, looking at the individual trajectories:
ggplot() +
geom_hline(yintercept = 0.8, linetype = "dotted", color = "tomato", size = 1) +
geom_line(data = trajectory,
aes(x = time, y = depletion, group = factor(seed * quantile), color = quantile), alpha = 0.05
) +
geom_line(data = trajectory %>%
group_by(time, quantile) %>%
summarize(mean = mean(depletion),
lower = as.numeric(quantile(depletion, probs = 0.025)),
upper = as.numeric(quantile(depletion, probs = 0.975))
),
aes(x = time, y = mean, group = quantile, color = quantile), size = 1
) +
scale_y_continuous(name = "Depletion (as % of K)",
breaks = seq(0.0, 1.0, 0.1),
labels = 100 * seq(0.0, 1.0, 0.1)
) +
ggtitle("Rmax = 4%, MNPL = 50%\nCatastrophic event wiping\nout 10% of animals") +
guides(color = "none") +
facet_wrap(~ quantile, ncol = 3) +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.