library(ggplot2)
library(dplyr)
#>
#> Attachement du package : 'dplyr'
#> Les objets suivants sont masqués depuis 'package:stats':
#>
#> filter, lag
#> Les objets suivants sont masqués depuis 'package:base':
#>
#> intersect, setdiff, setequal, union
library(RLA)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.
hp <- pellatomlinson_pbr(depletion0 = 0.3,
K = 5E5,
Rmax = 0.04,
catches = floor(runif(30, 1e3, 5e3)),
CV_env = 0.01,
seed_id = 20210306
)
#> No burn-in initial depletion level: 30 % of carrying capacity
#> Correlated random walk for environmental stochasticityWe’ve simulated a population of \(\frac{30}{100} \times \mathrm{K} = \frac{30}{100} \times 500000 = 150000\) animals, of which between \(1000\) and \(5000\) are removed annually because of human activities. The maximum population growth rate is \(\mathrm{R_{max}} = 4\)% and the \(\mathrm{R_{max}} = 60\)% of the carrying capacity \(\mathrm{K}\).
data.frame(time = rep(0:length(hp$removals), 2),
y = c(hp$depletion, 0, hp$removals/hp$K),
response = rep(c("depletion", "removal_rate"), each = length(hp$depletion))
) %>%
ggplot(aes(x = time, y = y, color = response)) +
geom_line() +
facet_wrap(~response, scales = "free_y") +
theme_bw()In this simulation, the population is starting to recover. However, we can change \(\mathrm{K}\) to be only \(250000\) and see what this changes:
hp <- pellatomlinson_pbr(depletion0 = 0.3,
K = 2.5E5,
Rmax = 0.04,
catches = floor(runif(30, 1e3, 5e3)),
CV_env = 0.01,
seed_id = 20210306
)
#> No burn-in initial depletion level: 30 % of carrying capacity
#> Correlated random walk for environmental stochasticity
data.frame(time = rep(0:length(hp$removals), 2),
y = c(hp$depletion, 0, hp$removals/hp$K),
response = rep(c("depletion", "removal_rate"), each = length(hp$depletion))
) %>%
ggplot(aes(x = time, y = y, color = response)) +
geom_line() +
facet_wrap(~response, scales = "free_y") +
theme_bw()The removals are now representing a larger fraction of the total population, and the population is declining.
In order to tune PBR to a conservation objective, you need to assess many scenarios, each with several simulations. This is illustrated below. We first need to generate a several simulations, here \(100\).
n_sim <- 100
set.seed(20210306)
scenario1 <- lapply(1:n_sim,
function(i) {
pellatomlinson_pbr(depletion = 0.3, Rmax = 0.04, CV_env = 0.01, verbose = FALSE)
}
)These can now be used with the function forward_pbr. Here we’ll start by evaluating the \(20\)% quantile of the best available abundance estimate (simulated as part of the function pellatomlinson_pbr):
q20 <- lapply(scenario1,
forward_pbr,
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:
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%, Quantile = 20%") +
guides(color = "none") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.The point of PBR is to test several quantiles to tune the procedure
quantile2evaluate <- seq(0.1, 0.9, 0.1)
tuning <- lapply(quantile2evaluate,
function(q) {
lapply(scenario1,
forward_pbr,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
F_r = 0.5,
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 = 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%") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.For example, in this case, only the smallest 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 now a catastrophic event that will happen once in those \(100\) years and that will kill immediately \(20\)% of the population.
tuning <- lapply(quantile2evaluate,
function(q) {
lapply(scenario1,
forward_pbr,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
catastrophe = 0.2,
F_r = 0.5,
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 = 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%, Catastrophic event wiping out 20% of animals") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.Here the catastrophe is happening always at the same moment across the different choice of quantiles for a given simulation in scenario1 because we’ve left the arg. random = FALSE in forward_pbr. This allows to control the seed for the simulations across the different choices of a quantile, and thus it is like playing the exact same simulations with only the choice of quantile that is different. This can be turned off by setting random = TRUE in forward_pbr.
tuning <- lapply(quantile2evaluate,
function(q) {
lapply(scenario1,
forward_pbr,
bycatch_variation_cv = 0.3,
distribution = "truncnorm",
catastrophe = 0.2,
F_r = 0.5,
q = q,
random = TRUE
)
}
)
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 = 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%, Catastrophic event wiping out 20% of animals") +
theme_bw()
#> `summarise()` has grouped output by 'time'. You can override using the `.groups` argument.The average curves do not look very different but the individual trajectories are:
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 = 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%, Catastrophic event wiping out 20% 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.This also shows that the forward_* functions in the package are systematically simulating a catastrophic event at the end of a MSFD cycle (the jagged pattern every six years). This is because it was convenient to code and easy to do so, and tweaking it would not make a big difference.
New simulations can thus easily be carried out, and tuning to new management objectives is straightforward.