sensitivity-analysis.RmdSensitivity analysis is performed to understand how changes in model parameters or inputs affect the model output. A One-way sensitivity analysis varies each model parameter or input individually across a known parameter space while holding all other inputs and parameters constant.
The model parameters varied for sensitivity analysis were all the parameters documented in the springRunDSM::params
Most of the model inputs are scaled by 0.5 to 1.5 in increments of 0.1. For example, if an input has an original value of 100, then for sensitivity analysis the model is run eleven times with the input set to the values: 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, and 150.
When the scaled input value need to be between 0 and 1 (ex: proportions or survival rates), the input is scaled as follows: \[Inv.Logit\left(log\left(\frac{x}{1-x}\right) \times scaler\right)\] where \(x\) is the model input being varied for sensitivity analysis, and \(scaler\) is the same scaling from .5 to 1.5 described above. In this example, if \(x = .3\) then the scaled inputs used in the sensitivity analysis would be 0.396, 0.376, 0.356, 0.337, 0.318, 0.3, 0.283, 0.266, 0.249, 0.234, and 0.219
The following model inputs are discrete values and handled as
special cases. For cc_gates_days_closed, first the values
of cc_gates_prop_days_closed (proportion of days within
a month the cross channel gates are closed) are scaled using the
inverse logit method described above and then multiplying this
proportion by the number of days in the month to produce the number of
days closed within a month. For weeks_flooded, the original
input is modified by adding \(\pm 2\)
weeks or \(\pm 1\) week, with a minimum
of 1 week and a maximum of 4 weeks.
The functions developed below automate as much of the sensitivity analysis as possible. It captures the scaling described above and allows the analysis to be run in parallel.
In order to run models in parallel, both parallel and
doParallel packages need to be installed.
library(springRunDSM)
# remotes::install_github("cvpia-osc/DSMscenario")
library(DSMscenario)
# install.packages("parallel")
# install.packages("doParallel")
library(parallel)
library(doParallel)
# install.packages("purrr")
library(purrr) # functional programming libraryNext, register the CPU cores that will be utilized for the parallel
process. Also specify any functions or data that will be used during the
parallel process. Here the run_scenario,
spring_run_model, and scenarios elements are
registered. These elements are defined further down in this
document.
no_cores <- detectCores(logical = TRUE) # total cores available
cl <- makeCluster(no_cores - 1) # create a cluster from cores
registerDoParallel(cl) # register cluster with package
# register the functions for use in parallel mode
clusterExport(cl, list('run_scenario', 'spring_run_model', 'scenarios'))The sensitivity analysis process is made up of a series of functions each doing a specific job. Below we describe each one of these and explain how they work together to create the sensitivity analysis.
The run_scenarios function calls the
spring_run_model in both the seeding and the full 20-year
simulation modes with the user provided scenario object and modified
model parameters. The function returns the valley-wide averages of total
number of spawners for each year in the simulation. This metric can be
modified to be whatever the user desires to compare between model runs
to gauge model sensitivity. The run_scenarios function is
called within the run_scenarios_scaled_param function
described in the next section.
run_scenario <- function(scenario, sensi_params) {
seeds <- spring_run_model(mode = "seed", ..params = sensi_params, stochastic = FALSE)
run <- spring_run_model(scenario = scenario,
mode = "simulate", seeds = seeds,
..params = sensi_params, stochastic = FALSE)
return(mean(colSums(run$spawners * run$proportion_natural, na.rm = TRUE)))
}The run_scenarios_scaled_param runs
run_scenarios in parallel with a scaled version of model
parameters (see section Scaling Methodology above for details).
The function selects the appropriate scaling methodology based on the
param name. This function returns a data frame with results
for each of the 13 SIT Candidate Restoration Strategies.
# create a list of the scenario for convenience
scenarios <- list(DSMscenario::scenarios$NO_ACTION, DSMscenario::scenarios$ONE,
DSMscenario::scenarios$TWO, DSMscenario::scenarios$THREE,
DSMscenario::scenarios$FOUR, DSMscenario::scenarios$FIVE,
DSMscenario::scenarios$SIX, DSMscenario::scenarios$SEVEN,
DSMscenario::scenarios$EIGHT, DSMscenario::scenarios$NINE,
DSMscenario::scenarios$TEN, DSMscenario::scenarios$ELEVEN,
DSMscenario::scenarios$TWELVE, DSMscenario::scenarios$THIRTEEN)
run_scenarios_scaled_param <- function(param, scalar, index = NULL) {
sensi_params <- springRunDSM::params
if (param %in% c("cc_gates_prop_days_closed", "cc_gates_days_closed")) {
# scale prop days closed using the 0-1 restriction
sensi_params["cc_gates_prop_days_closed"][[1]] <- boot::inv.logit(log((sensi_params["cc_gates_prop_days_closed"][[1]] + 1e-7) /
((1 - sensi_params["cc_gates_prop_days_closed"][[1]]) + 1e-7)) * scalar)
sensi_params["cc_gates_days_closed"][[1]] <- floor(lubridate::days_in_month(1:12) * sensi_params["cc_gates_prop_days_closed"][[1]])
} else {
sensi_params[param][[1]] <-
if (param %in% c("cross_channel_stray_rate",
"delta_prop_high_predation", "delta_proportion_diverted",
"growth_rates", "growth_rates_floodplain",
"hatchery_allocation", "mean_egg_temp_effect",
"migratory_temperature_proportion_over_20", "min_survival_rate",
"month_return_proportions", "natural_adult_removal_rate",
"prob_nest_scoured", "prob_strand_early", "prob_strand_late",
"prop_flow_natal",
"prop_high_predation", "prop_pulse_flows", "proportion_diverted",
"proportion_flow_bypass", "proportion_hatchery", "rear_decay_rate",
"spawn_decay_rate",
"spawn_success_sex_ratio", "stray_rate")) {
boot::inv.logit(log((sensi_params[param][[1]] + 1e-7) / ((1 - sensi_params[param][[1]]) + 1e-7)) * scalar)
} else if (param %in% c("weeks_flooded")) {
scalar
} else {
sensi_params[param][[1]] * scalar
}
}
scenario_results_list <- parLapply(cl, scenarios,
fun = function(scenario) {
run_scenario(scenario, sensi_params)
})
scenario_results <- unlist(scenario_results_list)
if (param == "weeks_flooded") {
wf <- c("+2", "-2", "1", "-1")
scalar <- wf[index]
} else {
scalar <- paste(round(scalar, 2), collapse = ",")
}
return(data.frame(param, scalar, base = scenario_results[1],
scenario_1 = scenario_results[2], scenario_2 = scenario_results[3],
scenario_3 = scenario_results[4], scenario_4 = scenario_results[5],
scenario_5 = scenario_results[6], scenario_6 = scenario_results[7],
scenario_7 = scenario_results[8], scenario_8 = scenario_results[9],
scenario_9 = scenario_results[10], scenario_10 = scenario_results[11],
scenario_11 = scenario_results[12], scenario_12 = scenario_results[13],
scenario_13 = scenario_results[14]))
}The param_sensitivity function brings everything
together, it takes a model parameter or input as an argument and
iterates through different scaling values each time running all 13 SIT
Candidate Restoration Strategies by calling
run_scenarios_scaled_param.
param_sensitivity <- function(param) {
scalars <- if (param == "weeks_flooded") {
(function() {
original <- springRunDSM::params$weeks_flooded
fp_filter <- springRunDSM::params$floodplain_habitat > 0
weeks_flooded_scaled <- list()
weeks_flooded_scaled[[1]] <- pmin(original + 2, 4) * fp_filter
weeks_flooded_scaled[[2]] <- pmax(original - 2, 1) * fp_filter
weeks_flooded_scaled[[3]] <- pmin(original + 1, 4) * fp_filter
weeks_flooded_scaled[[4]] <- pmax(original - 1, 1) * fp_filter
return(weeks_flooded_scaled)
})()
} else {
# default steps
seq(.5, 1.5, by = .1)
}
purrr::imap_dfr(scalars, ~run_scenarios_scaled_param(param, .x, .y))
}Below is an example of running sensitivity for one model input.
hatchery_allocation_results <- param_sensitivity("hatchery_allocation")
View(hatchery_allocation_results)