danstich/shadia

shadia: dam passage performance standard models for R

shadia

Package containing functions and data for American shad population modeling (click here for project website).

Note

This is the development version of the shadia package for R (shadia v2.0). This package is under continual development at this time, as it is part of ongoing research and management. For that reason, it is important to correspond with package developers before using the dam passage performance standard models in decision making or research.

Frozen (stable) versions of the software are provided for consistency during dam relisencing studies and reproducibility of research through the shadia website. The site now also provides up-to-date information about the status of all models.

Please check regularly for updates, changes, and bug patches. Please submit issues directly to this repository using the links above.

Installation

The development version of shadia can be installed with devtools in R using the repository url:

devtools::install_github("danStich/shadia")

To install shadia, you will need to have devtools installed ahead of time in R, but that requires some special tools. To install on Windows, you will need to download and install the appropriate version of Rtools. To install on Mac, you will need to have the XCode command-line tools installed. And, if running from Linux, you will need to install the developer version of R (r-base-dev) if you have not already.

Use

The purpose of this package is to distribute code used to run dam passage performance standard models for anadromous alosines (shad and herring). Currently, the model is implemented for the Connecticut, Kennebec, Merrimack, Mohawk and Hudson, Penobscot, Saco, and Susquehanna rivers, USA, but we are actively adding new rivers.

The main package functions can be run without any arguments to estimate population abundance in various reaches or in whole rivers under 'no dam' passage scenarios (see examples). Alternatively, the user can pass one or more values for upstream and downstream fish passage at a given dam which can then be applied throughout the watershed, or separately at each dam. Outputs include population abundance of spawners in the watersheds, within specific production units of each river, and the proportion of repeat spawners in each age class.

The models take several (10-30) seconds to run for one iteration on most standard workstations.

Warning

Management decisions should not be based on a single model run. The models rely on stochastic inputs for parameterization, as detailed in Stich et al. (2019). As such, any two model runs might result in substantially different predictions, even under the same passage scenario. We recommend at least 100 model runs per scenario to provide a minimal characterization of stochasticity, and a cursory understanding of variability in the response(s) of interest. In these cases, we strongly recommend running the model using the snowfall package as demonstrated in the help file for each model, which can be accessed by typing ?...RiverModel (where '...' is the name of each river in lowercase) in the console and pressing < Enter >.

Directories

data/ Contains built-in data sets for the package

man/ help files and documentation

R/ R functions in scripts

src/ C++ source files written with Rcpp

shadia's Issues

Add @param for initial spawning population size

Would like to add this as "number of adult spawners". John S suggested reverse-engineering stable age distribution from a Leslie matrix. Can be done pretty easy by importing from demogR package as we do in anadrofish.

populationSize doesn't match sum of PUS

populationSize is calculated as the sum of males and females across production units before application of pre-spawn mortality, fishing mortality, and post-spawn mortality.

Production-unit population sizes are currently calculated after application of those rates.

Should re-code R files used in penobscotRiverModel() and merrimackRiverModel() so that the sum of production-unit spawners equals populationSize.

Saco: dam passage being filled incorrectly

Upstream and downstream passage efficiencies are not being filled correctly during data write. May be occurring in fillOutputVectors.R or in writeData functions. The passage efficiencies are being implemented in the individual-based migration model (moveC) correctly as evidenced by movement of fish throughout the river.

Timing for watershed-scale models

Related to issue #24.

It appears that timing values for the first dam are not recycled at remaining dams as intended when watershed = TRUE.

Temporary workaround: can pre-sample the value for timing and then pass to function using object name, like this:

# Sample timing
times <- sample(c(1,2,3),1)  

# Run the model
res1 <- connecticutRiverModel(
  nRuns = 1,
  nYears = 2,
  timing = list(times, times, times, times, times),
  upstream = list(
    holyoke = c(0.50, 0.60, 0.70),
    cabot = 1,
    spillway = 1,
    gatehouse = 1,
    vernon = 1
  ),
  downstream = list(
    holyoke = c(0.80, 0.90, 0.95),
    cabot = 1,
    gatehouse= 1,
    vernon = 1
  ),
  northfield = list(
    turnersA = 0.95,
    turnersJ = 0.95,
    vernonA = 0.95,
    vernonJ = 0.95
  ),
  pSpillway = 1,
  inRiverF = 0,
  commercialF = 0,
  bycatchF = 0,
  indirect = 1,
  latent = 1,
  watershed = TRUE
  )

Parameter check

Need to remove unused parameters following code updates

Merrimack stock crashing after adding second passage route

Added a second passage route for Merrimack River. Following the addition, the population appears to blink out in < 20 years. Need to look at downstream survival calculation in downstreamMigration.R.

R code cleanup

Remove deprecated or unused R code from release version

User-specified arguments

Right now, changes to dam passage performance standards need to be specified in outer-loop-sampling.R

We need a way that these can be specified in the call to penobscotRiverModel() once the package is loaded. Ideally, this would be able to incorporate:

1. Dam-specific upstream and downstream passage efficiencies for individual dams (or all dams set to be the same).

2. Some standardization so that it can be passed from river to river as we develop more models.

Remote parallel execution

Need to develop a direction/examples for parallel execution on remote clusters. Right now, it is up to the user to know how to do this.

Bug: invalid length argument with no traceback

Looks like a bug was introduced in one of the variable storage routines during transition to a package. I think this may be related to garbage collection in R/innerLoopSampling.R. Need to look into it, because it causes cores to crash and execution is halted.

v-2.0-updates

Need to implement significant updates to the package prior to submission of new models to peer review. These will include issues
#18, #22, #24, #38, #46, #47

These and other changes will include:
[x] Change watershed default value to FALSE because this will be used more commonly moving forward
[ ] Implementation of the downstream delay parameter
[ ] Implementation of truck & transport and user-facing arguments
[ ] Implementation of juvenile stocking at arbitrary locations in the watershed and user-facing arguments
[x] Implementation of seasonal passage performance standard specification or derivation
[x] Separation of adult and juvenile downstream dam survival and user-facing argument for downstream_juv
[x] Memory management and speed improvements (reduction of redundant object returns)
[ ] Ability to read upstream and downstream passage efficiencies from .csv files
[x] Ability to specify starting population size as number of spawners instead of age-1 fish
[ ] Implementation of Salmon Falls River model
[ ] Add watershed output toggle to models other than mohawkHudsonRiverModel()

Generalize application of moveC and delayC

I want to generalize the application of the individual-based upstream migration model (moveC) and delay calculations (delayC) based on the number of routes (nRoutes). Right now, we use logical queries based on river name, but several of the routines are identical based on number of routes.

`popStart` not output correctly to `sens`

Starting population size for sensitivity analysis (popStart) is not being output to sens correctly for sensitivity analysis.

setUpData slowing down

As we add more data sets, the data set up is slowing down. A lot of the fat can be trimmed from some of the data sets to speed load because we don't use a lot of the data but manipulation is done in place for transparency. Can do ahead of time, comment out code, and leave in place unused for reproducibility.

Fish not moving into NY

Fish in the Susquehanna River model are not moving into NY waters in the 'North Branch' movement route (upstream_path == 4). Thought that it was due to temperatures or long distance from mouth of river, but I cranked up movement speed parameters and closed off other routes, and fish are still not making it. Might be fine, but need to take a closer look.

Daily vs seasonal passage

Need: functionality for users to supply seasonal passage efficiencies that we can convert to time-based rates.

Justification: we traditionally think about upstream and downstream passage in terms of seasonal efficiency and most projection models make assumptions about bulk efficiency to simulate effects of dams (habitat access) on populations. The models in shadia use individual-based migration models instead of (state-based) projection matrices. Upstream passage through dams is parameterized as a probability of passage per unit time (e.g. daily rates). This was done to align with how dam passage standards are prescribed and assessed (e.g. telemetry studies).

Dam passage efficiencies not aligning

During first watershed application for a re-licensing, I noticed that the passage efficiencies for each dam were being randomly drawn from the values of the first dam, instead of all being set to the same randomly drawn value (intention).

I have added patches to all of the main ...RiverModel functions to resolve the issue

Big population in first PU of Kennebec River

High numbers of spawners predicted for first PU. Possibly related to #40

Multiple values per dam for timing

Cannot currently pass multiple values per dam to the timing argument in ...RiverModel due to an error in data write.

Temporary workaround: can pass a vector, to be sampled randomly, which achieves the same result as passing a vector (these are sampled randomly by internal functions to evaluate all input combinations:

res1 <- connecticutRiverModel(
  nRuns = 1,
  nYears = 10,
  timing = list(sample(c(1,2,3),1),1,1,1,1),
  upstream = list(
    holyoke = c(0.50, 0.60, 0.70),
    cabot = 1,
    spillway = 1,
    gatehouse = 1,
    vernon = 1
  ),
  downstream = list(
    holyoke = c(0.80, 0.90, 0.95),
    cabot = 1,
    gatehouse= 1,
    vernon = 1
  ),
  northfield = list(
    turnersA = 0.95,
    turnersJ = 0.95,
    vernonA = 0.95,
    vernonJ = 0.95
  ),
  pSpillway = 1,
  inRiverF = 0,
  commercialF = 0,
  bycatchF = 0,
  indirect = 1,
  latent = 1,
  watershed = TRUE
  )

Make passage time dam-specific in all models

Would like the ability to assign different passage times to each dam in a system so that individual performance standards for passage time can be evaluated.

Juvenile downstream survival

Need to add ability to specify downstream survival through dams separately for juveniles and adults. Not sure what the best way is to do this. Could use a multiplier, but that would not allow for dam-specific juvenile survival evaluation. Can add as user-specified arguments, but the functions for each of the river models are getting really cumbersome to write out with all of the optional arguments...

Check scalars for CT River model

Base population numbers appear to have increased in CT River model following package updates.

Transient dynamics resulting from simStartingPop()

There are some odd transient dynamics occurring as a result of the new routine used to simulate the starting population.

Merrimack River needs second migration route

Following conversations with Michael Bailey, US Fish and Wildlife Service, Nashua, NH:

We need to add a second migration route for the Merrimack River. The second route needs to come in at the Pawtucket Dam, where there are two options for upstream and downstream migration. Should be two paths up, and two paths down. The two routes split about 1 km apart (upstream to downstream). Habitat is the same in both routes-- just need to sum within PUs across routes.

Comment updates

Need to update and double-check comments in .R files to make sure they are accurate.

Testing for kennebecRiverModel and sacoRiverModel

New models need testing to make sure things aren't crazy

Susquehanna River needs truck and transport + stocking

Need to add capacity for truck and transport and stocking scenarios to the Susquehanna River model.

Would like for this to be generalized so it is applicable to all river systems future and current.

CTR Fish not passing

Fish in the Connecticut River ABM are not passing upstream in the manner in which they should. Numbers in PUs 3-5 are shockingly low, and they do not match previous results from old code. This is probably a result of retroactively fitting this code into the PNR and MMR code. Need to have a look at this. Otherwise, the model is looking good.

Huge numbers in first PU for `mohawkHudsonRiverModel()`

The number of spawners being returned in the first downstream production unit seems unrealistically high in this model. Possibly related to similar bug/data issue #42 in kennebecRiverModel()

Documentation missing

Need to go through and integrate roxygen2 comments to compile documentation

Results cleanup

Need to clean up results so that:

1. only relevant outputs (excluding vars for sensitivity analysis) are returned in results list
2. names match input/output parameter notation from manuscript

Results check

Double check to make sure that results across temporal and spatial replicates are equal

Standardize naming conventions

This one is a long way off and fairly low priority right now. But, would be nice to standardize naming conventions between inputs and outputs of all river models.

Saco: temperature data QC

Need to run QC on Saco River temperature data (tempD_saco.rda).

Parallel implementation

Create wrapper for parallel execution on local work station