In TESA-workshops/TESAsamSim: Closed-loop salmon simulator for 2021 TESA workshop

This document illustrates how to simulate some stock-recruitment data using the TESAsamSim R package and then fit a linear state-space model (aka Dynamic Linear model) with potentially time-varying parameters to it. The model is fit via Maximum likelihood followed by Kalman filtering and smoothing using a wrapper for the dlm R package (see here). The wrapper is based in part on code generously provided by Cody Szuwalski (NOAA). For this illustrative example the simulations are based on the population characteristics of Interior Fraser River coho salmon.

First, load packages, read in data and set up scenarios

[Figure out how to do all this with purrr or furrr later]

## Check if necessary packages are available and install if necessary
listOfPackages <- c("here", "parallel", "doParallel", "foreach", 
                    "tidyverse", "tictoc", "TESAsamSim")
newPackages <- listOfPackages[!(listOfPackages %in% 
                                  installed.packages()[ , "Package"])]
if (length(newPackages)) {
  install.packages(newPackages)
}
lapply(listOfPackages, require, character.only = TRUE)

# load wrapper function to fit dynamic linear models
source(here("dlm-wrapper.R"))

simPar <- read.csv(here("data", "IFCohoPars",
                          "SimPar.csv"),
                        stringsAsFactors = F)

scenNames <- unique(simPar$scenario)

#if you need to get the sim files run this code (takes 15 mins)
#source("runTESAsamSim.R")

Put all scenarios in the same list then convert to a dataframe. Then create dataframes for output.

# read in list of files created in TESAsamSim/runTESAsamSim.R
simFiles <- list.files(here::here("outputs/simData/runTESA"))
nsc <- length(simFiles)

# Read the  sim scenarios into a list
simData <- list()

for(i in 1:nsc){
  filename <- here::here("outputs/simData/runTESA", simFiles[i])
  simData[[i]] <- readRDS(filename)$srDatout
}

nyr <- max(unique(simData[[i]]$year))


# filter for CU and rename cols
for(i in 1:nsc){
  simData[[i]] <- simData[[i]] %>% 
    dplyr::filter(CU == 1, year %in% (nyr-50+1):nyr) %>% 
    mutate() %>% 
    select(-CU) %>% 
    rename(byr=year, spwn=spawners, rec=recruits, alpha_true=alpha, beta_true=beta) %>% 
    mutate(scenario = scenNames[i], alpha=99., beta=99., alpha_se=99., beta_se=99.) %>% # cols for output
    select(scenario, everything())  #reorder cols
}

#, alpha=9999., beta=9999., alpha_se=9999., beta_se=9999.

# convert to df with all scenarios
dlm_Out <- simData[[1]]
for(i in 2:nsc) dlm_Out <- rbind(dlm_Out, simData[[i]])

# Now make copies to cover our alternative estimation models
estNames <- c("1_Stat", "2_Alpha_vary", "3_Beta_vary", "4_Alpha_Beta_vary")
dlm_Out_stat <- dlm_Out_alpha <- dlm_Out_beta <- dlm_Out_alphabeta <- dlm_Out

Now fit the models, borrowing from Brendan's code.

iter <- unique(simData[[1]]$iteration)
nsc <- length(scenNames)

#check to see if models have already been run (this file is in the Google drive)
dlm_filename <- here::here("outputs/diagnostics",   "model_estimates_all_combos.csv")

if(file.exists(dlm_filename)){
  dlm_out_all_combo <- readr::read_csv(dlm_filename)
} else{
  for(j in 1:nsc){
    for(i in seq_along(iter)){

      dat <- dlm_Out %>% 
        dplyr::filter(scenario == scenNames[j], iteration==i) %>% 
        select(-c(alpha,beta,alpha_se,beta_se)) #need to remove these for fitting

      # alpha and beta fixed in estimation model
      dlm_model_stat <- fitDLM(data = dat,
                          alpha_vary = FALSE,
                          beta_vary = FALSE)

      dlm_Out_stat[which(dlm_Out_stat$scenario==scenNames[j] & dlm_Out_stat$iteration==i),10:13] <- dlm_model_stat$results[,10:13]

      # alpha varies in estimation model
      dlm_model_alpha <- fitDLM(data = dat,
                          alpha_vary = TRUE,
                          beta_vary = FALSE)

      dlm_Out_alpha[which(dlm_Out_alpha$scenario==scenNames[j] & dlm_Out_alpha$iteration==i),10:13] <- dlm_model_alpha$results[,10:13]

      # beta varies in estimation model
      dlm_model_beta <- fitDLM(data = dat,
                          alpha_vary = FALSE,
                          beta_vary = TRUE)

      dlm_Out_beta[which(dlm_Out_beta$scenario==scenNames[j] & dlm_Out_beta$iteration==i),10:13] <- dlm_model_beta$results[,10:13]

      # alpha and beta vary in estimation model
      dlm_model_alphabeta <- fitDLM(data = dat,
                          alpha_vary = TRUE,
                          beta_vary = TRUE)

      dlm_Out_alphabeta[which(dlm_Out_alphabeta$scenario==scenNames[j] & dlm_Out_alphabeta$iteration==i),10:13] <- dlm_model_alphabeta$results[,10:13]

    }#end j
  }# end i

  # Now append the estimation model name to each dataframe
dlm_Out_stat <- dlm_Out_stat %>% 
  mutate(estModel = estNames[1])
dlm_Out_alpha <- dlm_Out_alpha %>% 
  mutate(estModel = estNames[2])
dlm_Out_beta <- dlm_Out_beta %>% 
  mutate(estModel = estNames[3])
dlm_Out_alphabeta <- dlm_Out_alphabeta %>% 
  mutate(estModel = estNames[4])

# Now make one gigantic dataframe (not sure if we want this)
dlm_out_all_combo <- rbind(dlm_Out_stat, dlm_Out_alpha, dlm_Out_beta,dlm_Out_alphabeta)

readr::write_csv(dlm_out_all_combo, here::here("outputs/diagnostics",   "model_estimates_all_combos.csv"))

}# end if

So now all the outputs are in one giant dataframe

Get the bias for all combos

# Now get bias ... adapted from Brendan's code

# change sign of beta_true
dlm_out_all_combo$beta_true <- -dlm_out_all_combo$beta_true


allbias <- dlm_out_all_combo %>%
  group_by(scenario,estModel, iteration) %>%
  dplyr::summarize(
    alpha_mpb=mean((alpha_true-alpha)/alpha_true)*100,
    beta_mpb=mean((beta_true-beta)/beta_true)*100) %>%
  pivot_longer(alpha_mpb:beta_mpb,names_to="parameter",values_to="mpb")

# View(allbias)

Plot bias.

source(here::here("plotOutputs_newfuncs.R"))


# Subset the scenarios and make box plots
plotBias <- allbias %>% 
  filter(scenario %in% scenNames[1:5])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc1-5.png"),  width = 960, height = 720)
  plot_bias_boxplots(plotBias)
dev.off()

plotBias <- allbias %>% 
  filter(scenario %in% scenNames[6:10])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc6-10.png"), width = 960, height = 720 )
  plot_bias_boxplots(plotBias)
dev.off()

plotBias <- allbias %>% 
  filter(scenario %in% scenNames[11:15])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc11-15.png"), width = 960, height = 720 )
  plot_bias_boxplots(plotBias)
dev.off()

plotBias <- allbias %>% 
  filter(scenario %in% scenNames[16:20])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc16-20.png"), width = 960, height = 720 )
  plot_bias_boxplots(plotBias)
dev.off()

plotBias <- allbias %>% 
  filter(scenario %in% scenNames[21:25])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc21-25.png"), width = 960, height = 720 )
  plot_bias_boxplots(plotBias)
dev.off()

plotBias <- allbias %>% 
  filter(scenario %in% scenNames[26:31])
png(here::here("outputs/diagnostics","Bias_boxplots_Sc26-31.png"), width = 960, height = 720 )
  plot_bias_boxplots(plotBias)
dev.off()