In andrew-edwards/EDMsimulate: Simulating Fish Populations for Empirical Dynamic Modeling
knitr::opts_chunk$set(echo = TRUE)
Plot realisations
This code runs the sim_and_fit_realisations over numerous trials and extracts the forecasts
library(EDMsimulate)
library(dplyr)
library(ggplot2)
library(larkin)
library(parallel)
library(tidyr)
library(RColorBrewer)
# Run using stan-optim for Larkin model to speed up model runs
# See larkin-compareMLE-Bayes.Rmd for a comparison
# larkin_args() and ricker_args() have default values, but if one element
# changes, the entire list needs to be included here
# # List of elements for running function separately
# salmon_sim_args <- list(sigma_nu = 1, omega = 0.4)
# target_hr <- 0.2
# sigma_ou <- 0.05
# edm_fit <- TRUE
# larkin_fit <- TRUE
# ricker_fit <- TRUE
# mvs_fit <- FALSE
# R_switch <- "R_t"#R_prime_t"#"R_t"#
# pbsEDM_args <- list(lags = list(R_switch = 0,
# S_t = 0:3),
# first_difference = TRUE,
# centre_and_scale = TRUE)
# mve_args <- list(lags = list(R_t = 0:4,
# S_t = 0:8))
# larkin_args = list( run_stan = FALSE,
# prior_mean_alpha = 2,
# prior_mean_beta = -rep(1,4),
# prior_mean_sigma = 0.5,
# prior_sd_alpha = 0.5,
# prior_sd_beta = rep(0.25,4),
# prior_sd_sigma = 0.25)
# ricker_args <- list(run_stan = FALSE,
# prior_mean_alpha = 2,
# prior_mean_beta = -1,
# prior_mean_sigma = 0.5,
# prior_sd_alpha = 0.5,
# prior_sd_beta = 0.25,
# prior_sd_sigma = 0.25)
# M <- 2
# do_parallel <- TRUE#FALSE
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 1,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
larkin_fit = TRUE,
ricker_fit = TRUE,
R_switch = "R_t",#"R_prime_t",#"R_t",#
mvs_fit = FALSE,
M = 2)
Then plot outputs.
# Extract forecasts and forecast errors
results_bc <- forecast_errors(out, label="base_case")
if (file.exists(here::here("report/out")) == FALSE){
dir.create(here::here("report/out"))
}
write.csv(results_bc, here::here("report/out/results_bc.csv"))
# Plot forecasts against true (simulated) values
p1 <- results_bc$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_bc$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_bc$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density")
p4 <- results_bc$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
The above example assumes returns are the variable being predicted. Another option is recruitment (by brood year), which is the more direct predictor of Larkin model. The simulation is repeated with recruitment, and EDM forecasts are (tentatively) less biased (why?)
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 1,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
larkin_fit = TRUE,
ricker_fit = TRUE,
R_switch = "R_prime_t",#"R_t",#
M = 100)
# Extract forecasts and forecast errors
results_rec <- forecast_errors(out, label="recruits")
# Plot forecasts against true (simulated) values
p1 <- results_rec$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_rec$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_rec$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density") + xlim(-1000,1000)
p4 <- results_rec$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
The Larkin and Ricker models tend to have absolute errors in forecast that are closer to zero, compared with EDM.
The simulation is then repeated, but with altered beta parameters where beta0 is dominant over the remaining beta parameters (more similar to a Ricker model).
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 1,
omega = 0.4,
beta = c(0.45, 0.25, 0.15, 0.15)),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_ric <- forecast_errors(out, label="ricker")
p1 <- results_ric$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_ric$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_ric$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density") + xlim(-1000,1000)
p4 <- results_ric$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
Simulations were then run with low process variation (sigma_nu=0.1)
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 0.1,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
R_switch = "R_t",#"R_t",#
M = 1000)
# Extract forecasts and forecast errors
results_sigma_nu_0.1 <- forecast_errors(out, label="sigma_nu_0.1")
write.csv(results_sigma_nu_0.1, here::here("report/out/results_sigma_nu_0.1.csv"))
# Plot forecasts against true (simulated) values
p1 <- results_sigma_nu_0.1$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density") + xlim(-100,100)
p4 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2)
Simulations were then run including temporal variability in productivity (Larkin alpha parameter), with a trends similar to those observed historically for Fraser River Sockeye salmon (Pestal et al. in review; Huang et al. 2020 [Fraser Sockeye RPA]). Three patterns in temporal variability in productivity are evident in the historical time-series: (1) stable productivity throughout most of the time-series with the exception of declines in the most recent generation (e.g., Nadina, Seymour, Gates, Upper Pitt River, Scotch, Chilko, Late Shuswap) , (2) gradual declines starting in the 1970s (e.g.,Upper Barriere, Portage, Early Stuart, Bowron), and (3) relatively stable productivity in early time-series until 1990 followed by steeper declines (e.g., Weaver Creek, Raft, Birkenhead, Stellako, Quesnel, Late Stuart, Weaver Creek).
# Fraser Sockeye time-series run 1950-2016 BY, so 1970 is T minus 46 (or year 34 in 80 year time-series) and 1990 is T minus 26 (or year 54 in 80 year time-series)
alpha1 <- c(rep(exp(2), 100), rep(exp(2), 76), seq(exp(2), exp(1.75), length.out=5))
alpha2 <- c(rep(exp(2), 100), rep(exp(2), 34), seq(exp(2), exp(0.5), length.out=47))
alpha3 <- c(rep(exp(2), 100), rep(exp(2), 54), seq(exp(2), exp(0.5), length.out=27))
alpha4 <- c(rep(exp(2), 100), rep(exp(2), 34), rep(exp(0.5),47))
df <- data.frame(year = rep(1:181,4),
productivity = c(alpha1, alpha2, alpha3, alpha4),
trend=c(rep(1,181), rep(2, 181), rep(3,181), rep(4,181)))
p1 <- df |> ggplot2::ggplot(aes(x=year, y=productivity, colour=as.factor(trend))) + geom_line()
p1
These trends were then included in the simulation runs
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(alpha = alpha2,
sigma_nu = 1,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 500)
# Extract forecasts and forecast errors
results_alpha_decline500 <- forecast_errors(out, label = "alpha_decline")
write.csv(results_alpha_decline500, here::here("report/out/results_alpha_decline500.csv"))
p1 <- results_alpha_decline$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_alpha_decline$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_alpha_decline$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density") + xlim(-1000,1000)
p4 <- results_alpha_decline$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
#And repeated with low process variation and alpha decline.
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(alpha = alpha2,
sigma_nu = 0.1,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 500)
# Extract forecasts and forecast errors
results_alpha_decline_sigma_nu_0.1_500 <- forecast_errors(out, label = "alpha_decline_sigma_nu_0.1")
write.csv(results_alpha_decline_sigma_nu_0.1_500, here::here("report/out/results_alpha_decline_sigma_nu_0.1_500.csv"))
p1 <- results_alpha_decline_sigma_nu_0.1$forecasts %>%
ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) +
geom_point() +
geom_smooth(method='lm', formula= y~x) +
xlab("Simulated 'true' recruitment") +
ylab("Forecasted recruitment") +
geom_abline(slope=1, intercept = 0, col="dark grey")
# Plot distribution of absolute errors in forecast
p2 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "Err") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Errors in forecast") +
ylab("Density")
p3 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Percent errors in forecast") +
ylab("Density") + xlim(-100,100)
p4 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr") %>%
ggplot(aes(x=value, colour=forMethod)) +
geom_density() +
xlab("Absolute errors in forecast") +
ylab("Density")
ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
Or, results can be presented from sensitivity analyses together
#Aggregate the df's of discrete sensitivity analyses, Percent Error
discrete_sens_pe <- results_bc$errors %>% dplyr::filter(error == "perErr") %>%
add_row(results_rec$errors %>% dplyr::filter(error == "perErr")) %>%
add_row(results_ric$errors %>% dplyr::filter(error == "perErr")) %>%
add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>%
add_row(results_alpha_decline$errors %>% dplyr::filter(error == "perErr")) %>%
add_row(results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("base_case", "recruits", "ricker",
"sigma_nu_0.1", "alpha_decline",
"alpha_decline_sigma_nu_0.1")))
discrete_sens_ae <- results_bc$errors %>% dplyr::filter(error == "absErr") %>%
add_row(results_rec$errors %>% dplyr::filter(error == "absErr")) %>%
add_row(results_ric$errors %>% dplyr::filter(error == "absErr")) %>%
add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>%
add_row(results_alpha_decline$errors %>% dplyr::filter(error == "absErr")) %>%
add_row(results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("base_case", "recruits", "ricker",
"sigma_nu_0.1", "alpha_decline",
"alpha_decline_sigma_nu_0.1")))
dodge_pe <- position_dodge(width = 0.6)
dodge_ae <- position_dodge(width = 0.7)
p1 <- discrete_sens_pe %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+
geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) +
ylim(-100,100) + guides(alpha = FALSE) + ylab("Percent error") +
ggtitle("Discrete sensitivity analyses")
p2 <- discrete_sens_ae %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+
geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) +
ylim(0,1) + guides(alpha = FALSE) + ylab("Absolute error") +
ggtitle("Discrete sensitivity analyses")
p1
p2
Run sensitivity analyses over a gradient in sigma_nu
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 0.3,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_sigma_nu_0.3 <- forecast_errors(out, label = "sigma_nu_0.3")
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 0.5,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_sigma_nu_0.5 <- forecast_errors(out, label = "sigma_nu_0.5")
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(sigma_nu = 0.8,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_sigma_nu_0.8 <- forecast_errors(out, label = "sigma_nu_0.8")
#### Get saved results:
# results_bc <- read.csv(here::here('report/out/results_bc.csv'),
# stringsAsFactors = F,check.names = F)
#
# results_bc$errors <- results_bc %>% select(grep('errors',names(results_bc)))
#
# names(results_bc$errors) <- gsub(x = names(results_bc$errors),
# pattern = "errors.", replacement = "")
# results_bc$forecasts <- results_bc %>%
# select(grep('forecasts',names(results_bc)))
#
# names(results_bc$forecasts) <- gsub(x = names(results_bc$forecasts),
# pattern = "forecasts.", replacement = "")
#
# results_sigma_nu_0.1 <-
# read.csv(here::here('report/out/results_sigma_nu_0.1.csv'),
# stringsAsFactors = F,check.names = F)
#
# results_sigma_nu_0.1$errors <-
# results_sigma_nu_0.1 %>% select(grep('errors',names(results_sigma_nu_0.1)))
#
# names(results_sigma_nu_0.1$errors) <-
# gsub(x = names(results_sigma_nu_0.1$errors), pattern = "errors.",
# replacement = "")
#
# results_sigma_nu_0.1$forecasts <-
# results_sigma_nu_0.1 %>%
# select(grep('forecasts',names(results_sigma_nu_0.1)))
#
# names(results_sigma_nu_0.1$forecasts) <-
# gsub(x = names(results_sigma_nu_0.1$forecasts), pattern = "forecasts.",
# replacement = "")
# Compile results from various model runs
sigma_nu_sens_pe <- results_bc$errors %>% dplyr::filter(error == "perErr") %>%
add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>%
# add_row(results_sigma_nu_0.3$errors %>% dplyr::filter(error == "perErr")) %>%
# add_row(results_sigma_nu_0.5$errors %>% dplyr::filter(error == "perErr")) %>%
# add_row(results_sigma_nu_0.8$errors %>% dplyr::filter(error == "perErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("sigma_nu_0.1",
# "sigma_nu_0.3", "sigma_nu_0.5",
# "sigma_nu_0.8",
"base_case")))
sigma_nu_sens_pe <- sigma_nu_sens_pe %>% dplyr::filter(forMethod == "mve"|
forMethod == "lar"|
forMethod =="ric")
sigma_nu_sens_pe$forMethod <- factor(sigma_nu_sens_pe$forMethod, levels = c("ric", "lar", "mve"))
sigma_nu_sens_ae <- results_bc$errors %>% dplyr::filter(error == "absErr") %>%
add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>%
# add_row(results_sigma_nu_0.3$errors %>% dplyr::filter(error == "absErr")) %>%
# add_row(results_sigma_nu_0.5$errors %>% dplyr::filter(error == "absErr")) %>%
# add_row(results_sigma_nu_0.8$errors %>% dplyr::filter(error == "absErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("sigma_nu_0.1",
# "sigma_nu_0.3", "sigma_nu_0.5",
# "sigma_nu_0.8",
"base_case")))
sigma_nu_sens_ae <- sigma_nu_sens_ae %>% dplyr::filter(forMethod == "mve"|
forMethod=="lar"|
forMethod =="ric")
sigma_nu_sens_ae$forMethod <- factor(sigma_nu_sens_ae$forMethod, levels = c("ric", "lar", "mve"))
# dodge_pe <- position_dodge(width = 0.6)
# dodge_ae <- position_dodge(width = 0.7)
dodge_pe <- position_dodge(width = 0.8)
dodge_ae <- position_dodge(width = 0.8)
p1 <- sigma_nu_sens_pe %>%
ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+
geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) +
ylim(-100,150) + guides(alpha = FALSE) + ylab("Percent error") +
labs(title = "Effects of increasing sigma_nu (recruitment residual error)",
subtitle = "") +
scale_x_discrete(labels = c("sigma_nu_0.1" = "0.1",
"sigma_nu_0.3" = "0.3",
"sigma_nu_0.5" = "0.5",
"sigma_nu_0.8" = "0.8",
"base_case" = "1.0")) +
labs(x="SD of recruitment residuals") +
scale_fill_discrete(name = element_blank(),
labels = c("Ricker", "Larkin", "MVE")) +
geom_hline(yintercept = 0, col=grey(0.4)) +
theme_classic()
p2 <- sigma_nu_sens_ae %>%
ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+
geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) +
ylim(0,1.5) + guides(alpha = FALSE) + ylab("Absolute error") +
labs(title = "Effects of increasing sigma_nu (recruitment residual error)",
subtitle = "")+
scale_x_discrete(labels = c("sigma_nu_0.1" = "0.1",
"sigma_nu_0.3" = "0.3",
"sigma_nu_0.5" = "0.5",
"sigma_nu_0.8" = "0.8",
"base_case" = "1.0")) +
labs(x="SD of recruitment residuals") +
scale_fill_discrete(name = element_blank(),
labels = c("Ricker", "Larkin", "MVE")) + theme_classic()
p1
p2
# ggsave(filename = here::here("report/out/pe.png"), p1,
# width = 6, height = 5, units ="in", dpi=350)
# ggsave(filename = here::here("report/out/ae.png"), p2,
# width = 6, height = 5, units ="in", dpi=350)
And declines in productivity and a gradient in sigma_nu
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(alpha = alpha2,
sigma_nu = 0.3,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_alpha_decline_sigma_nu_0.3 <- forecast_errors(out, label = "sigma_nu_0.3")
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(alpha = alpha2,
sigma_nu = 0.5,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_alpha_decline_sigma_nu_0.5 <- forecast_errors(out, label = "sigma_nu_0.5")
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args =
list(alpha = alpha2,
sigma_nu = 0.8,
omega = 0.4),
target_hr = 0.2,
sigma_ou = 0.05,
R_switch = "R_t",#"R_prime_t",#"R_t",#
larkin_fit = TRUE,
ricker_fit = TRUE,
mvs_fit = FALSE,
M = 100)
# Extract forecasts and forecast errors
results_alpha_decline_sigma_nu_0.8 <- forecast_errors(out, label = "sigma_nu_0.8")
alpha_decline_sigma_nu_sens_pe <- results_alpha_decline500$errors %>%
dplyr::filter(error == "perErr") %>%
add_row(results_alpha_decline_sigma_nu_0.1_500$errors %>%
dplyr::filter(error == "perErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.3$errors %>% dplyr::filter(error == "perErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.5$errors %>% dplyr::filter(error == "perErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.8$errors %>% dplyr::filter(error == "perErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("alpha_decline_sigma_nu_0.1",
# "sigma_nu_0.3", "sigma_nu_0.5",
# "sigma_nu_0.8",
"alpha_decline")))
alpha_decline_sigma_nu_sens_pe <- alpha_decline_sigma_nu_sens_pe %>%
dplyr::filter(forMethod == "mve"|
forMethod=="lar"|
forMethod=="ric")
alpha_decline_sigma_nu_sens_pe$forMethod <-
factor( alpha_decline_sigma_nu_sens_pe$forMethod,
levels = c("ric", "lar", "mve") )
alpha_decline_sigma_nu_sens_ae <- results_alpha_decline500$errors %>%
dplyr::filter(error == "absErr") %>%
add_row(results_alpha_decline_sigma_nu_0.1_500$errors %>%
dplyr::filter(error == "absErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.3$errors %>% dplyr::filter(error == "absErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.5$errors %>% dplyr::filter(error == "absErr")) %>%
# add_row(results_alpha_decline_sigma_nu_0.8$errors %>% dplyr::filter(error == "absErr")) %>%
dplyr::mutate(model_label=factor(model_label,
levels=c("alpha_decline_sigma_nu_0.1",
# "sigma_nu_0.3", "sigma_nu_0.5",
# "sigma_nu_0.8",
"alpha_decline")))
alpha_decline_sigma_nu_sens_ae <- alpha_decline_sigma_nu_sens_ae %>%
dplyr::filter(forMethod == "mve"|
forMethod=="lar"|
forMethod=="ric")
alpha_decline_sigma_nu_sens_ae$forMethod <-
factor( alpha_decline_sigma_nu_sens_ae$forMethod,
levels = c("ric", "lar", "mve") )
dodge_pe <- position_dodge(width = 0.8)
dodge_ae <- position_dodge(width = 0.8)
p1 <- alpha_decline_sigma_nu_sens_pe %>%
ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+
geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) +
ylim(-100,200) + guides(alpha = FALSE) + ylab("Percent error") +
labs(title = "Effects of increasing sigma_nu (recruitment residual error)",
subtitle = "under scenario of declining alpha over time") +
scale_x_discrete(labels = c("alpha_decline_sigma_nu_0.1" = "0.1",
"sigma_nu_0.3" = "0.3",
"sigma_nu_0.5" = "0.5",
"sigma_nu_0.8" = "0.8",
"alpha_decline" = "1.0")) +
labs(x="SD of recruitment residuals") +
scale_fill_discrete(name = element_blank(),
labels = c("Ricker", "Larkin", "MVE")) +
geom_hline(yintercept = 0, col=grey(0.4)) +
theme_classic()
p2 <- alpha_decline_sigma_nu_sens_ae %>%
ggplot(aes(x = model_label, y = value, fill = forMethod)) +
geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+
geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) +
ylim(0,1.5) + guides(alpha = FALSE) + ylab("Absolute error") +
labs(title = "Effects of increasing sigma_nu (recruitment residual error)",
subtitle = "under scenario of declining alpha over time")+
scale_x_discrete(labels = c("alpha_decline_sigma_nu_0.1" = "0.1",
"sigma_nu_0.3" = "0.3",
"sigma_nu_0.5" = "0.5",
"sigma_nu_0.8" = "0.8",
"alpha_decline" = "1.0")) +
labs(x="SD of recruitment residuals") +
scale_fill_discrete(name = element_blank(),
labels = c("Ricker", "Larkin", "MVE")) + theme_classic()
p1
p2
# ggsave(filename = here::here("report/out/pe_alpha_decline.png"), p1,
# width = 6, height = 5, units ="in", dpi=350)
# ggsave(filename = here::here("report/out/ae_alpha_decline.png"), p2,
# width = 6, height = 5, units ="in", dpi=350)
andrew-edwards/EDMsimulate documentation built on Oct. 25, 2023, 2:43 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
Plot realisations
This code runs the sim_and_fit_realisations over numerous trials and extracts the forecasts
library(EDMsimulate) library(dplyr) library(ggplot2) library(larkin) library(parallel) library(tidyr) library(RColorBrewer) # Run using stan-optim for Larkin model to speed up model runs # See larkin-compareMLE-Bayes.Rmd for a comparison # larkin_args() and ricker_args() have default values, but if one element # changes, the entire list needs to be included here # # List of elements for running function separately # salmon_sim_args <- list(sigma_nu = 1, omega = 0.4) # target_hr <- 0.2 # sigma_ou <- 0.05 # edm_fit <- TRUE # larkin_fit <- TRUE # ricker_fit <- TRUE # mvs_fit <- FALSE # R_switch <- "R_t"#R_prime_t"#"R_t"# # pbsEDM_args <- list(lags = list(R_switch = 0, # S_t = 0:3), # first_difference = TRUE, # centre_and_scale = TRUE) # mve_args <- list(lags = list(R_t = 0:4, # S_t = 0:8)) # larkin_args = list( run_stan = FALSE, # prior_mean_alpha = 2, # prior_mean_beta = -rep(1,4), # prior_mean_sigma = 0.5, # prior_sd_alpha = 0.5, # prior_sd_beta = rep(0.25,4), # prior_sd_sigma = 0.25) # ricker_args <- list(run_stan = FALSE, # prior_mean_alpha = 2, # prior_mean_beta = -1, # prior_mean_sigma = 0.5, # prior_sd_alpha = 0.5, # prior_sd_beta = 0.25, # prior_sd_sigma = 0.25) # M <- 2 # do_parallel <- TRUE#FALSE out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 1, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, larkin_fit = TRUE, ricker_fit = TRUE, R_switch = "R_t",#"R_prime_t",#"R_t",# mvs_fit = FALSE, M = 2)
Then plot outputs.
# Extract forecasts and forecast errors results_bc <- forecast_errors(out, label="base_case") if (file.exists(here::here("report/out")) == FALSE){ dir.create(here::here("report/out")) } write.csv(results_bc, here::here("report/out/results_bc.csv")) # Plot forecasts against true (simulated) values p1 <- results_bc$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_bc$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_bc$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") p4 <- results_bc$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
The above example assumes returns are the variable being predicted. Another option is recruitment (by brood year), which is the more direct predictor of Larkin model. The simulation is repeated with recruitment, and EDM forecasts are (tentatively) less biased (why?)
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 1, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, larkin_fit = TRUE, ricker_fit = TRUE, R_switch = "R_prime_t",#"R_t",# M = 100) # Extract forecasts and forecast errors results_rec <- forecast_errors(out, label="recruits") # Plot forecasts against true (simulated) values p1 <- results_rec$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_rec$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_rec$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") + xlim(-1000,1000) p4 <- results_rec$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
The Larkin and Ricker models tend to have absolute errors in forecast that are closer to zero, compared with EDM.
The simulation is then repeated, but with altered beta parameters where beta0 is dominant over the remaining beta parameters (more similar to a Ricker model).
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 1, omega = 0.4, beta = c(0.45, 0.25, 0.15, 0.15)), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_ric <- forecast_errors(out, label="ricker") p1 <- results_ric$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_ric$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_ric$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") + xlim(-1000,1000) p4 <- results_ric$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
Simulations were then run with low process variation (sigma_nu=0.1)
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 0.1, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, R_switch = "R_t",#"R_t",# M = 1000) # Extract forecasts and forecast errors results_sigma_nu_0.1 <- forecast_errors(out, label="sigma_nu_0.1") write.csv(results_sigma_nu_0.1, here::here("report/out/results_sigma_nu_0.1.csv")) # Plot forecasts against true (simulated) values p1 <- results_sigma_nu_0.1$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") + xlim(-100,100) p4 <- results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2)
Simulations were then run including temporal variability in productivity (Larkin alpha parameter), with a trends similar to those observed historically for Fraser River Sockeye salmon (Pestal et al. in review; Huang et al. 2020 [Fraser Sockeye RPA]). Three patterns in temporal variability in productivity are evident in the historical time-series: (1) stable productivity throughout most of the time-series with the exception of declines in the most recent generation (e.g., Nadina, Seymour, Gates, Upper Pitt River, Scotch, Chilko, Late Shuswap) , (2) gradual declines starting in the 1970s (e.g.,Upper Barriere, Portage, Early Stuart, Bowron), and (3) relatively stable productivity in early time-series until 1990 followed by steeper declines (e.g., Weaver Creek, Raft, Birkenhead, Stellako, Quesnel, Late Stuart, Weaver Creek).
# Fraser Sockeye time-series run 1950-2016 BY, so 1970 is T minus 46 (or year 34 in 80 year time-series) and 1990 is T minus 26 (or year 54 in 80 year time-series) alpha1 <- c(rep(exp(2), 100), rep(exp(2), 76), seq(exp(2), exp(1.75), length.out=5)) alpha2 <- c(rep(exp(2), 100), rep(exp(2), 34), seq(exp(2), exp(0.5), length.out=47)) alpha3 <- c(rep(exp(2), 100), rep(exp(2), 54), seq(exp(2), exp(0.5), length.out=27)) alpha4 <- c(rep(exp(2), 100), rep(exp(2), 34), rep(exp(0.5),47)) df <- data.frame(year = rep(1:181,4), productivity = c(alpha1, alpha2, alpha3, alpha4), trend=c(rep(1,181), rep(2, 181), rep(3,181), rep(4,181))) p1 <- df |> ggplot2::ggplot(aes(x=year, y=productivity, colour=as.factor(trend))) + geom_line() p1
These trends were then included in the simulation runs
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(alpha = alpha2, sigma_nu = 1, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 500) # Extract forecasts and forecast errors results_alpha_decline500 <- forecast_errors(out, label = "alpha_decline") write.csv(results_alpha_decline500, here::here("report/out/results_alpha_decline500.csv")) p1 <- results_alpha_decline$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_alpha_decline$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_alpha_decline$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") + xlim(-1000,1000) p4 <- results_alpha_decline$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 ) #And repeated with low process variation and alpha decline. out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(alpha = alpha2, sigma_nu = 0.1, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 500) # Extract forecasts and forecast errors results_alpha_decline_sigma_nu_0.1_500 <- forecast_errors(out, label = "alpha_decline_sigma_nu_0.1") write.csv(results_alpha_decline_sigma_nu_0.1_500, here::here("report/out/results_alpha_decline_sigma_nu_0.1_500.csv")) p1 <- results_alpha_decline_sigma_nu_0.1$forecasts %>% ggplot(aes(x=R_sim, y=R_fit, colour=forMethod)) + geom_point() + geom_smooth(method='lm', formula= y~x) + xlab("Simulated 'true' recruitment") + ylab("Forecasted recruitment") + geom_abline(slope=1, intercept = 0, col="dark grey") # Plot distribution of absolute errors in forecast p2 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "Err") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Errors in forecast") + ylab("Density") p3 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Percent errors in forecast") + ylab("Density") + xlim(-100,100) p4 <- results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr") %>% ggplot(aes(x=value, colour=forMethod)) + geom_density() + xlab("Absolute errors in forecast") + ylab("Density") ggpubr::ggarrange(p1, p2, p3, p4, ncol=2, nrow=2 )
Or, results can be presented from sensitivity analyses together
#Aggregate the df's of discrete sensitivity analyses, Percent Error discrete_sens_pe <- results_bc$errors %>% dplyr::filter(error == "perErr") %>% add_row(results_rec$errors %>% dplyr::filter(error == "perErr")) %>% add_row(results_ric$errors %>% dplyr::filter(error == "perErr")) %>% add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>% add_row(results_alpha_decline$errors %>% dplyr::filter(error == "perErr")) %>% add_row(results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("base_case", "recruits", "ricker", "sigma_nu_0.1", "alpha_decline", "alpha_decline_sigma_nu_0.1"))) discrete_sens_ae <- results_bc$errors %>% dplyr::filter(error == "absErr") %>% add_row(results_rec$errors %>% dplyr::filter(error == "absErr")) %>% add_row(results_ric$errors %>% dplyr::filter(error == "absErr")) %>% add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>% add_row(results_alpha_decline$errors %>% dplyr::filter(error == "absErr")) %>% add_row(results_alpha_decline_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("base_case", "recruits", "ricker", "sigma_nu_0.1", "alpha_decline", "alpha_decline_sigma_nu_0.1"))) dodge_pe <- position_dodge(width = 0.6) dodge_ae <- position_dodge(width = 0.7) p1 <- discrete_sens_pe %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+ geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) + ylim(-100,100) + guides(alpha = FALSE) + ylab("Percent error") + ggtitle("Discrete sensitivity analyses") p2 <- discrete_sens_ae %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+ geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) + ylim(0,1) + guides(alpha = FALSE) + ylab("Absolute error") + ggtitle("Discrete sensitivity analyses") p1 p2
Run sensitivity analyses over a gradient in sigma_nu
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 0.3, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_sigma_nu_0.3 <- forecast_errors(out, label = "sigma_nu_0.3") out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 0.5, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_sigma_nu_0.5 <- forecast_errors(out, label = "sigma_nu_0.5") out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(sigma_nu = 0.8, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_sigma_nu_0.8 <- forecast_errors(out, label = "sigma_nu_0.8") #### Get saved results: # results_bc <- read.csv(here::here('report/out/results_bc.csv'), # stringsAsFactors = F,check.names = F) # # results_bc$errors <- results_bc %>% select(grep('errors',names(results_bc))) # # names(results_bc$errors) <- gsub(x = names(results_bc$errors), # pattern = "errors.", replacement = "") # results_bc$forecasts <- results_bc %>% # select(grep('forecasts',names(results_bc))) # # names(results_bc$forecasts) <- gsub(x = names(results_bc$forecasts), # pattern = "forecasts.", replacement = "") # # results_sigma_nu_0.1 <- # read.csv(here::here('report/out/results_sigma_nu_0.1.csv'), # stringsAsFactors = F,check.names = F) # # results_sigma_nu_0.1$errors <- # results_sigma_nu_0.1 %>% select(grep('errors',names(results_sigma_nu_0.1))) # # names(results_sigma_nu_0.1$errors) <- # gsub(x = names(results_sigma_nu_0.1$errors), pattern = "errors.", # replacement = "") # # results_sigma_nu_0.1$forecasts <- # results_sigma_nu_0.1 %>% # select(grep('forecasts',names(results_sigma_nu_0.1))) # # names(results_sigma_nu_0.1$forecasts) <- # gsub(x = names(results_sigma_nu_0.1$forecasts), pattern = "forecasts.", # replacement = "") # Compile results from various model runs sigma_nu_sens_pe <- results_bc$errors %>% dplyr::filter(error == "perErr") %>% add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_sigma_nu_0.3$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_sigma_nu_0.5$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_sigma_nu_0.8$errors %>% dplyr::filter(error == "perErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("sigma_nu_0.1", # "sigma_nu_0.3", "sigma_nu_0.5", # "sigma_nu_0.8", "base_case"))) sigma_nu_sens_pe <- sigma_nu_sens_pe %>% dplyr::filter(forMethod == "mve"| forMethod == "lar"| forMethod =="ric") sigma_nu_sens_pe$forMethod <- factor(sigma_nu_sens_pe$forMethod, levels = c("ric", "lar", "mve")) sigma_nu_sens_ae <- results_bc$errors %>% dplyr::filter(error == "absErr") %>% add_row(results_sigma_nu_0.1$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_sigma_nu_0.3$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_sigma_nu_0.5$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_sigma_nu_0.8$errors %>% dplyr::filter(error == "absErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("sigma_nu_0.1", # "sigma_nu_0.3", "sigma_nu_0.5", # "sigma_nu_0.8", "base_case"))) sigma_nu_sens_ae <- sigma_nu_sens_ae %>% dplyr::filter(forMethod == "mve"| forMethod=="lar"| forMethod =="ric") sigma_nu_sens_ae$forMethod <- factor(sigma_nu_sens_ae$forMethod, levels = c("ric", "lar", "mve")) # dodge_pe <- position_dodge(width = 0.6) # dodge_ae <- position_dodge(width = 0.7) dodge_pe <- position_dodge(width = 0.8) dodge_ae <- position_dodge(width = 0.8) p1 <- sigma_nu_sens_pe %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+ geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) + ylim(-100,150) + guides(alpha = FALSE) + ylab("Percent error") + labs(title = "Effects of increasing sigma_nu (recruitment residual error)", subtitle = "") + scale_x_discrete(labels = c("sigma_nu_0.1" = "0.1", "sigma_nu_0.3" = "0.3", "sigma_nu_0.5" = "0.5", "sigma_nu_0.8" = "0.8", "base_case" = "1.0")) + labs(x="SD of recruitment residuals") + scale_fill_discrete(name = element_blank(), labels = c("Ricker", "Larkin", "MVE")) + geom_hline(yintercept = 0, col=grey(0.4)) + theme_classic() p2 <- sigma_nu_sens_ae %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+ geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) + ylim(0,1.5) + guides(alpha = FALSE) + ylab("Absolute error") + labs(title = "Effects of increasing sigma_nu (recruitment residual error)", subtitle = "")+ scale_x_discrete(labels = c("sigma_nu_0.1" = "0.1", "sigma_nu_0.3" = "0.3", "sigma_nu_0.5" = "0.5", "sigma_nu_0.8" = "0.8", "base_case" = "1.0")) + labs(x="SD of recruitment residuals") + scale_fill_discrete(name = element_blank(), labels = c("Ricker", "Larkin", "MVE")) + theme_classic() p1 p2 # ggsave(filename = here::here("report/out/pe.png"), p1, # width = 6, height = 5, units ="in", dpi=350) # ggsave(filename = here::here("report/out/ae.png"), p2, # width = 6, height = 5, units ="in", dpi=350)
And declines in productivity and a gradient in sigma_nu
out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(alpha = alpha2, sigma_nu = 0.3, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_alpha_decline_sigma_nu_0.3 <- forecast_errors(out, label = "sigma_nu_0.3") out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(alpha = alpha2, sigma_nu = 0.5, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_alpha_decline_sigma_nu_0.5 <- forecast_errors(out, label = "sigma_nu_0.5") out <- EDMsimulate::sim_and_fit_realisations(salmon_sim_args = list(alpha = alpha2, sigma_nu = 0.8, omega = 0.4), target_hr = 0.2, sigma_ou = 0.05, R_switch = "R_t",#"R_prime_t",#"R_t",# larkin_fit = TRUE, ricker_fit = TRUE, mvs_fit = FALSE, M = 100) # Extract forecasts and forecast errors results_alpha_decline_sigma_nu_0.8 <- forecast_errors(out, label = "sigma_nu_0.8") alpha_decline_sigma_nu_sens_pe <- results_alpha_decline500$errors %>% dplyr::filter(error == "perErr") %>% add_row(results_alpha_decline_sigma_nu_0.1_500$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.3$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.5$errors %>% dplyr::filter(error == "perErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.8$errors %>% dplyr::filter(error == "perErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("alpha_decline_sigma_nu_0.1", # "sigma_nu_0.3", "sigma_nu_0.5", # "sigma_nu_0.8", "alpha_decline"))) alpha_decline_sigma_nu_sens_pe <- alpha_decline_sigma_nu_sens_pe %>% dplyr::filter(forMethod == "mve"| forMethod=="lar"| forMethod=="ric") alpha_decline_sigma_nu_sens_pe$forMethod <- factor( alpha_decline_sigma_nu_sens_pe$forMethod, levels = c("ric", "lar", "mve") ) alpha_decline_sigma_nu_sens_ae <- results_alpha_decline500$errors %>% dplyr::filter(error == "absErr") %>% add_row(results_alpha_decline_sigma_nu_0.1_500$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.3$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.5$errors %>% dplyr::filter(error == "absErr")) %>% # add_row(results_alpha_decline_sigma_nu_0.8$errors %>% dplyr::filter(error == "absErr")) %>% dplyr::mutate(model_label=factor(model_label, levels=c("alpha_decline_sigma_nu_0.1", # "sigma_nu_0.3", "sigma_nu_0.5", # "sigma_nu_0.8", "alpha_decline"))) alpha_decline_sigma_nu_sens_ae <- alpha_decline_sigma_nu_sens_ae %>% dplyr::filter(forMethod == "mve"| forMethod=="lar"| forMethod=="ric") alpha_decline_sigma_nu_sens_ae$forMethod <- factor( alpha_decline_sigma_nu_sens_ae$forMethod, levels = c("ric", "lar", "mve") ) dodge_pe <- position_dodge(width = 0.8) dodge_ae <- position_dodge(width = 0.8) p1 <- alpha_decline_sigma_nu_sens_pe %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_pe, width=1.6, aes(alpha = 0.5))+ geom_boxplot(width = .15, outlier.colour = NA, position = dodge_pe) + ylim(-100,200) + guides(alpha = FALSE) + ylab("Percent error") + labs(title = "Effects of increasing sigma_nu (recruitment residual error)", subtitle = "under scenario of declining alpha over time") + scale_x_discrete(labels = c("alpha_decline_sigma_nu_0.1" = "0.1", "sigma_nu_0.3" = "0.3", "sigma_nu_0.5" = "0.5", "sigma_nu_0.8" = "0.8", "alpha_decline" = "1.0")) + labs(x="SD of recruitment residuals") + scale_fill_discrete(name = element_blank(), labels = c("Ricker", "Larkin", "MVE")) + geom_hline(yintercept = 0, col=grey(0.4)) + theme_classic() p2 <- alpha_decline_sigma_nu_sens_ae %>% ggplot(aes(x = model_label, y = value, fill = forMethod)) + geom_violin(position = dodge_ae, width=1.5, aes(alpha = 0.5))+ geom_boxplot(width = .2, outlier.colour = NA, position = dodge_ae) + ylim(0,1.5) + guides(alpha = FALSE) + ylab("Absolute error") + labs(title = "Effects of increasing sigma_nu (recruitment residual error)", subtitle = "under scenario of declining alpha over time")+ scale_x_discrete(labels = c("alpha_decline_sigma_nu_0.1" = "0.1", "sigma_nu_0.3" = "0.3", "sigma_nu_0.5" = "0.5", "sigma_nu_0.8" = "0.8", "alpha_decline" = "1.0")) + labs(x="SD of recruitment residuals") + scale_fill_discrete(name = element_blank(), labels = c("Ricker", "Larkin", "MVE")) + theme_classic() p1 p2 # ggsave(filename = here::here("report/out/pe_alpha_decline.png"), p1, # width = 6, height = 5, units ="in", dpi=350) # ggsave(filename = here::here("report/out/ae_alpha_decline.png"), p2, # width = 6, height = 5, units ="in", dpi=350)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.