Objective

Sensitivity 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

Scaling Methodology

  1. 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.

  2. 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

  3. 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.

Running Sensitivity Analysis

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.

Package Dependencies and Parallel Mode Set Up

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 library

Next, 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'))

Sensitivility Analysis Process

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.

Run Scenario

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)))
}

Run Scenarios with Scaled Parameters

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]))
}

Runnning Sensitivity

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)