From 6d4cec95eff127705496992deca04820cfd6c6c1 Mon Sep 17 00:00:00 2001 From: Sarah Gaichas Date: Wed, 13 Nov 2024 18:06:45 -0500 Subject: [PATCH] first draft zoop and temp covariates WP --- WP_RecCovariates_Gaichas.Rmd | 595 ++++++- docs/WP_RecCovariates_Gaichas.html | 2308 ++++++++++++++++++++++++++++ docs/index.Rmd | 3 + docs/index.html | 2 + 4 files changed, 2901 insertions(+), 7 deletions(-) create mode 100644 docs/WP_RecCovariates_Gaichas.html diff --git a/WP_RecCovariates_Gaichas.Rmd b/WP_RecCovariates_Gaichas.Rmd index 99a50c3..8df9274 100644 --- a/WP_RecCovariates_Gaichas.Rmd +++ b/WP_RecCovariates_Gaichas.Rmd @@ -47,13 +47,13 @@ We are using the `devel` version of WHAM: https://github.com/timjmiller/wham/tre Model [mm192](https://drive.google.com/drive/folders/1sQdDsfdnVbiiY4X7Rgr-fvegwT7Fa1Az?usp=drive_link) is our starting point. -Haddock egg predation, Zooplankton, and temperature indices were explored as covariates on herring recruitment. +Haddock egg predation, Zooplankton, and temperature indices were explored as covariates on herring recruitment. Methods for developing each indicator are in separate working papers. Recruitment is modeled as deviations from the "recruitment scaling parameter", leaving one option for modeling effects of covariates on recruitment: "controlling". A "controlling" recruitment covariate results in a time-varying recruitment scaling parameter. -[Jon's text on the haddock egg predation testing here... including trying to fit a stock recruitment function that didnt work] +[merge with Jon's ecisting text on the haddock egg predation testing here... ] We explored indices with different zooplankton groups, seasons, and regions according to herring life history and results from the boosted regression tree: @@ -62,28 +62,609 @@ We explored indices with different zooplankton groups, seasons, and regions acco * Sep-Feb small copepods in herring larval area with lag-1 to represent food for larvae more specifically * Combinations of large and small copepod covariates above +In addition, we explored an index of optimal temperature duration (days) during the fall larval season, September-December + We evaluated * Options for covariate input (millions of cells vs. log(cells), VAST estimated SE vs. WHAM estimated SE) * Options for covariate observation model ("rw" vs. "ar1") -* Options for recruitment link ("none" vs. "controlling-linear" with lag-0 for large copepods and lag-1 for small) +* Options for recruitment link ("none" vs. "controlling-linear" with lag-0 for large copepods, lag-1 for small copepods, and lag-1 for larval temperature duration) +* No attempts were made to fit polynomial effects although this is possible in WHAM + +The main diagnostics we used to determine if the model was improved by covariates were: + +* model converged and Hessian matrix invertable +* dAIC lowest +* recruitment sigma reduced (how much?) +* estimated covariate effect CI does not include 0 +* direction of covariate effect is sensible + +Sensitivity runs + +Turn off NAA RE and then try fitting with the most robust recruitment covariates. + + +# Results Short story: * Models with covariates input on the log scale generally converged * Models with WHAM estimated covariate SE ("est_1") generally converged -* Under the above conditions, most models with and without the recruitment link converged for all covariates +* Under the above conditions, most models with and without the recruitment link converged * Models with the Jan-Jun (Spring) large copepods covariate also converged with input as millions of cells and VAST estimated SE -* I'm still figuring out where to find all the diagnostics in WHAM, so "converged" may not be "a good model" - [Way too much detail including false starts](https://noaa-edab.github.io/zooplanktonindex/WHAMcovariate_tests.html) +Summary table with each model, recruitment sigma compared with base, covariate beta with CI -# Results + + +We show diagnostics and results from the best fit models for each covariate in sections below. + +Diagnostics include: + +* WHAM's fit to the covariate time series +* One step ahead residuals for WHAM's fit +* Estimated time varying recruitment scaling parameter with N age 1 + +## Spring large copepods + +*Models with no covariates had slightly better AIC than comparable models with recruitment links* + +```{r, results='hide'} + +config <- "lgcopeSpring2" + +# name the model directory +name <- paste0("mm192_", config) + +write.dir <- here::here(sprintf("WHAMfits/%s/", name)) + +# larger set of ecov setups to compare +df.mods <- data.frame(Recruitment = c(rep(2, 16)), + ecov_process = c(rep("rw",8),rep("ar1",8)), + ecov_how = c(rep("none",4), + rep("controlling-lag-0-linear",4)), + ecovdat = rep(c("logmean-logsig", + "logmean-est_1", + "meanmil-logsigmil", + "meanmil-est_1"),4), + stringsAsFactors=FALSE) +n.mods <- dim(df.mods)[1] +df.mods$Model <- paste0("m",1:n.mods) +df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col + + +mod.list <- paste0(write.dir, df.mods$Model,".rds") +mods <- lapply(mod.list, readRDS) + +n.mods <- length(mod.list) + +vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x))) +for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n') + + +opt_conv = 1-sapply(mods, function(x) x$opt$convergence) +ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0) +df.mods$conv <- as.logical(opt_conv) +df.mods$pdHess <- as.logical(ok_sdrep) +df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3)) + +not_conv <- !df.mods$conv | !df.mods$pdHess +mods2 <- mods +mods2[not_conv] <- NULL +df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab) +df.aic <- df.aic.tmp[FALSE,] +ct = 1 +for(i in 1:n.mods){ + if(not_conv[i]){ + df.aic[i,] <- rep(NA,5) + } else { + df.aic[i,] <- df.aic.tmp[ct,] + ct <- ct + 1 + } +} +df.mods <- cbind(df.mods, df.aic) +df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),] +df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---" +rownames(df.mods) <- NULL + +top4 <- df.mods |> + dplyr::filter(ecovdat %in% ("logmean-est_1")) |> + dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC) + +``` + +```{r} + +flextable::flextable(top4) |> + flextable::set_header_labels(ecov_process = "Process", + ecov_how = "Rec. link") |> + flextable::set_table_properties(layout = "autofit") +``` + + +### Spring large copepods covariate: logscale ar1 diagnostics + +The WHAM ar1 fit looks strange + +![lgCopeSp2 m10 fit](WHAMfits/mm192_lgcopeSpring2/m10/plots_png/diagnostics/lgCopeSpring2_diagnostic.png) + + +![lgCopeSp2 m10 osa](WHAMfits/mm192_lgcopeSpring2/m10/plots_png/diagnostics/OSA_resid_ecov_4panel_lgCopeSpring2.png) + + + +### Spring large copepods covariate: logscale ar1 recruitment + + +![lgCopeSp2 m10 rec](WHAMfits/mm192_lgcopeSpring2/m10/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +![lgCopeSp2 m14 rec](WHAMfits/mm192_lgcopeSpring2/m14/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + +```{r} +m10par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m10/res_tables/parameter_estimates_table.RDS")) + +m14par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m14/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m10par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m14par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; lgCopeSpring2 beta_1 is `r round(m14par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",1],3)`, CI `r round(m14par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",3],3)`, `r round(m14par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",4],3)` + +### Spring large copepods covariate: logscale rw diagnostics + + +![lgCopeSp2 m2 fit](WHAMfits/mm192_lgcopeSpring2/m2/plots_png/diagnostics/lgCopeSpring2_diagnostic.png) + + +![lgCopeSp2 m2 osa](WHAMfits/mm192_lgcopeSpring2/m2/plots_png/diagnostics/OSA_resid_ecov_4panel_lgCopeSpring2.png) + + + +### Spring large copepods covariate: logscale rw recruitment + + +![lgCopeSp2 m10 rec](WHAMfits/mm192_lgcopeSpring2/m2/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + + +![lgCopeSp2 m14 rec](WHAMfits/mm192_lgcopeSpring2/m6/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +```{r} +m2par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m2/res_tables/parameter_estimates_table.RDS")) + +m6par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m6/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m2par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m6par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; lgCopeSpring2 beta_1 is `r round(m6par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",1],3)`, CI `r round(m6par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",3],3)`, `r round(m6par["stock 1 Recruitment Ecov: lgCopeSpring2 $\\beta_1$",4],3)` + + + + +## Fall small copepods + +*Models with no covariates had slightly better AIC than the ar1 model with recruitment links* + +```{r, results='hide'} + +config <- "smcopeFall2" + +# name the model directory +name <- paste0("mm192_", config) + +write.dir <- here::here(sprintf("WHAMfits/%s/", name)) + +# larger set of ecov setups to compare +df.mods <- data.frame(Recruitment = c(rep(2, 16)), + ecov_process = c(rep("rw",8),rep("ar1",8)), + ecov_how = c(rep("none",4), + rep("controlling-lag-0-linear",4)), + ecovdat = rep(c("logmean-logsig", + "logmean-est_1", + "meanmil-logsigmil", + "meanmil-est_1"),4), + stringsAsFactors=FALSE) +n.mods <- dim(df.mods)[1] +df.mods$Model <- paste0("m",1:n.mods) +df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col + + +mod.list <- paste0(write.dir, df.mods$Model,".rds") +mods <- lapply(mod.list, readRDS) + +n.mods <- length(mod.list) + +vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x))) +for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n') + + +opt_conv = 1-sapply(mods, function(x) x$opt$convergence) +ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0) +df.mods$conv <- as.logical(opt_conv) +df.mods$pdHess <- as.logical(ok_sdrep) +df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3)) + +not_conv <- !df.mods$conv | !df.mods$pdHess +mods2 <- mods +mods2[not_conv] <- NULL +df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab) +df.aic <- df.aic.tmp[FALSE,] +ct = 1 +for(i in 1:n.mods){ + if(not_conv[i]){ + df.aic[i,] <- rep(NA,5) + } else { + df.aic[i,] <- df.aic.tmp[ct,] + ct <- ct + 1 + } +} +df.mods <- cbind(df.mods, df.aic) +df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),] +df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---" +rownames(df.mods) <- NULL + +top4 <- df.mods |> + dplyr::filter(ecovdat %in% ("logmean-est_1")) |> + dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC) + +``` + +```{r} + +flextable::flextable(top4) |> + flextable::set_header_labels(ecov_process = "Process", + ecov_how = "Rec. link") |> + flextable::set_table_properties(layout = "autofit") +``` + +### Fall small copepods covariate: logscale ar1 diagnostics + +![smCopeFall2 m10 fit](WHAMfits/mm192_smcopeFall2/m10/plots_png/diagnostics/smCopeFall2_diagnostic.png) + + +![smCopeFall2 m10 osa](WHAMfits/mm192_smcopeFall2/m10/plots_png/diagnostics/OSA_resid_ecov_4panel_smCopeFall2.png) + + +### Fall small copepods covariate: logscale ar1 recruitment + +![smCopeFall2 m10 rec](WHAMfits/mm192_smcopeFall2/m10/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +![smCopeFall2 m14 rec](WHAMfits/mm192_smcopeFall2/m14/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +```{r} +m10par <- readRDS(here::here("WHAMfits/mm192_smcopeFall2/m10/res_tables/parameter_estimates_table.RDS")) + +m14par <- readRDS(here::here("WHAMfits/mm192_smcopeFall2/m14/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m10par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m14par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; smcopeFall2 beta_1 is `r round(m14par["stock 1 Recruitment Ecov: smCopeFall2 $\\beta_1$",1],3)`, CI `r round(m14par["stock 1 Recruitment Ecov: smCopeFall2 $\\beta_1$",3],3)`, `r round(m14par["stock 1 Recruitment Ecov: smCopeFall2 $\\beta_1$",4],3)` + + +## September - February small copepods in the larval herring area + +*Model with ar1 covariate had slightly better AIC than models without recruitment links* + +```{r, results='hide'} + +config <- "smcopeSepFeb2" + +# name the model directory +name <- paste0("mm192_", config) + +write.dir <- here::here(sprintf("WHAMfits/%s/", name)) + +# larger set of ecov setups to compare +df.mods <- data.frame(Recruitment = c(rep(2, 16)), + ecov_process = c(rep("rw",8),rep("ar1",8)), + ecov_how = c(rep("none",4), + rep("controlling-lag-1-linear",4)), + ecovdat = rep(c("logmean-logsig", + "logmean-est_1", + "meanmil-logsigmil", + "meanmil-est_1"),4), + stringsAsFactors=FALSE) +n.mods <- dim(df.mods)[1] +df.mods$Model <- paste0("m",1:n.mods) +df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col + + +mod.list <- paste0(write.dir, df.mods$Model,".rds") +mod.list <- mod.list[-16] # mod 16 didnt run +mods <- lapply(mod.list, readRDS) + +df.mods <- df.mods[-16,] # mod 16 didnt run + +n.mods <- length(mod.list) + +vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x))) +for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n') + + +opt_conv = 1-sapply(mods, function(x) x$opt$convergence) +ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0) +df.mods$conv <- as.logical(opt_conv) +df.mods$pdHess <- as.logical(ok_sdrep) +df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3)) + +not_conv <- !df.mods$conv | !df.mods$pdHess +mods2 <- mods +mods2[not_conv] <- NULL +df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab) +df.aic <- df.aic.tmp[FALSE,] +ct = 1 +for(i in 1:n.mods){ + if(not_conv[i]){ + df.aic[i,] <- rep(NA,5) + } else { + df.aic[i,] <- df.aic.tmp[ct,] + ct <- ct + 1 + } +} +df.mods <- cbind(df.mods, df.aic) +df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),] +df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---" +rownames(df.mods) <- NULL + +top4 <- df.mods |> + dplyr::filter(ecovdat %in% ("logmean-est_1")) |> + dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC) + +``` + +```{r} + +flextable::flextable(top4) |> + flextable::set_header_labels(ecov_process = "Process", + ecov_how = "Rec. link") |> + flextable::set_table_properties(layout = "autofit") +``` + +### Sep-Feb small copepods in herring larval area covariate: logscale ar1 diagnostics + + +![smCopeSepFeb2 m10 fit](WHAMfits/mm192_smcopeSepFeb2/m10/plots_png/diagnostics/smCopeSepFeb2_diagnostic.png) + + +![smCopeSepFeb2 m10 osa](WHAMfits/mm192_smcopeSepFeb2/m10/plots_png/diagnostics/OSA_resid_ecov_4panel_smCopeSepFeb2.png) + + +### Sep-Feb small copepods in herring larval area covariate: logscale ar1 recruitment + + +![smCopeSepFeb2 m10 rec](WHAMfits/mm192_smcopeSepFeb2/m10/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +![smCopeSepFeb2 m14 rec](WHAMfits/mm192_smcopeSepFeb2/m14/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +```{r} +m10par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2/m10/res_tables/parameter_estimates_table.RDS")) + +m14par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2/m14/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m10par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m14par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; smcopeSepFeb2 beta_1 is `r round(m14par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",1],3)`, CI `r round(m14par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",3],3)`, `r round(m14par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",4],3)` + +## Alternate: Sep-Feb small copepods in fall herring survey strata covariate + +*Both rw and ar1 models worked with covariates, rw better fit* + + +```{r, results='hide'} + +config <- "smcopeSepFeb2_fallstrata" + +# name the model directory +name <- paste0("mm192_", config) + +write.dir <- here::here(sprintf("WHAMfits/%s/", name)) + +# larger set of ecov setups to compare +df.mods <- data.frame(Recruitment = c(rep(2, 16)), + ecov_process = c(rep("rw",8),rep("ar1",8)), + ecov_how = c(rep("none",4), + rep("controlling-lag-1-linear",4)), + ecovdat = rep(c("logmean-logsig", + "logmean-est_1", + "meanmil-logsigmil", + "meanmil-est_1"),4), + stringsAsFactors=FALSE) +n.mods <- dim(df.mods)[1] +df.mods$Model <- paste0("m",1:n.mods) +df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col + + +mod.list <- paste0(write.dir, df.mods$Model,".rds") +mod.list <- mod.list[-16] # mod 16 didnt run +mods <- lapply(mod.list, readRDS) + +df.mods <- df.mods[-16,] # mod 16 didnt run + +n.mods <- length(mod.list) + +vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x))) +for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n') + + +opt_conv = 1-sapply(mods, function(x) x$opt$convergence) +ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0) +df.mods$conv <- as.logical(opt_conv) +df.mods$pdHess <- as.logical(ok_sdrep) +df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3)) + +not_conv <- !df.mods$conv | !df.mods$pdHess +mods2 <- mods +mods2[not_conv] <- NULL +df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab) +df.aic <- df.aic.tmp[FALSE,] +ct = 1 +for(i in 1:n.mods){ + if(not_conv[i]){ + df.aic[i,] <- rep(NA,5) + } else { + df.aic[i,] <- df.aic.tmp[ct,] + ct <- ct + 1 + } +} +df.mods <- cbind(df.mods, df.aic) +df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),] +df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---" +rownames(df.mods) <- NULL + +top4 <- df.mods |> + dplyr::filter(ecovdat %in% ("logmean-est_1")) |> + dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC) + +``` + +```{r} + +flextable::flextable(top4) |> + flextable::set_header_labels(ecov_process = "Process", + ecov_how = "Rec. link") |> + flextable::set_table_properties(layout = "autofit") +``` + +### Sep-Feb small copepods in herring fall survey area covariate: logscale rw diagnostics + + +![smCopeSepFeb2_fallstrat m2 fit](WHAMfits/mm192_smcopeSepFeb2_fallstrata/m2/plots_png/diagnostics/smCopeSepFeb2_diagnostic.png) + +![smCopeSepFeb2_fallstrat m2 osa](WHAMfits/mm192_smcopeSepFeb2_fallstrata/m2/plots_png/diagnostics/OSA_resid_ecov_4panel_smCopeSepFeb2.png) + +### Sep-Feb small copepods in herring fall survey area covariate: logscale rw recruitment + + +![smCopeSepFeb2_fallstrat m2 rec](WHAMfits/mm192_smcopeSepFeb2_fallstrata/m2/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + +![smCopeSepFeb2_fallstrat m6 rec](WHAMfits/mm192_smcopeSepFeb2_fallstrata/m6/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +```{r} +m2par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2_fallstrata/m2/res_tables/parameter_estimates_table.RDS")) + +m6par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2_fallstrata/m6/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m2par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m6par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; smcopeSepFeb2 beta_1 is `r round(m6par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",1],3)`, CI `r round(m6par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",3],3)`, `r round(m6par["stock 1 Recruitment Ecov: smCopeSepFeb2 $\\beta_1$",4],3)` + + +## September - December duration of larval optimal temperature + +*Both rw and ar1 models worked with covariates, rw better fit* + +```{r, results='hide'} + +config <- "LarvalTempDuration" + +# name the model directory +name <- paste0("mm192_", config) + +write.dir <- here::here(sprintf("WHAMfits/%s/", name)) + +# larger set of ecov setups to compare +# larger set of ecov setups to compare +df.mods <- data.frame(Recruitment = c(rep(2, 8)), + ecov_process = c(rep("rw",4),rep("ar1",4)), + ecov_how = rep(c("none","controlling-lag-1-linear"), 4), + ecovdat = c(rep("mean-est_1", 2),rep("logmean-est_1",2)), + stringsAsFactors=FALSE) +n.mods <- dim(df.mods)[1] +df.mods$Model <- paste0("m",1:n.mods) +df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col + + +mod.list <- paste0(write.dir, df.mods$Model,".rds") +mods <- lapply(mod.list, readRDS) + +n.mods <- length(mod.list) + +vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x))) +for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n') + + +opt_conv = 1-sapply(mods, function(x) x$opt$convergence) +ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0) +df.mods$conv <- as.logical(opt_conv) +df.mods$pdHess <- as.logical(ok_sdrep) +df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3)) + +not_conv <- !df.mods$conv | !df.mods$pdHess +mods2 <- mods +mods2[not_conv] <- NULL +df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab) +df.aic <- df.aic.tmp[FALSE,] +ct = 1 +for(i in 1:n.mods){ + if(not_conv[i]){ + df.aic[i,] <- rep(NA,5) + } else { + df.aic[i,] <- df.aic.tmp[ct,] + ct <- ct + 1 + } +} +df.mods <- cbind(df.mods, df.aic) +df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),] +df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---" +rownames(df.mods) <- NULL + +top4 <- df.mods |> + dplyr::filter(ecovdat %in% ("logmean-est_1")) |> + dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC) + +``` + +```{r} + +flextable::flextable(top4) |> + flextable::set_header_labels(ecov_process = "Process", + ecov_how = "Rec. link") |> + flextable::set_table_properties(layout = "autofit") +``` + +### Duration of optimal larval temperature, Sept-Dec: logscale rw diagnostics + + +![LarvalTempDuration/m3 fit](WHAMfits/mm192_LarvalTempDuration/m3/plots_png/diagnostics/LarvalTempDuration_diagnostic.png) + + +![LarvalTempDuration/m3 osa](WHAMfits/mm192_LarvalTempDuration/m3/plots_png/diagnostics/OSA_resid_ecov_4panel_LarvalTempDuration.png) + + + +### Duration of optimal larval temperature, Sept-Dec: logscale rw recruitment + + +![LarvalTempDuration/m3 rec](WHAMfits/mm192_LarvalTempDuration/m3/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + +![LarvalTempDuration/m4 rec](WHAMfits/mm192_LarvalTempDuration/m4/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png) + + +```{r} +m3par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration/m3/res_tables/parameter_estimates_table.RDS")) + +m4par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration/m4/res_tables/parameter_estimates_table.RDS")) + +``` + +Without covariate, recruitment variance is `r round(m3par["stock 1 NAA $\\sigma$ (age 1)",1],3)`, and with is `r round(m4par["stock 1 NAA $\\sigma$ (age 1)",1],3)`; LarvalTempDuration beta_1 is `r round(m4par["stock 1 Recruitment Ecov: LarvalTempDuration $\\beta_1$",1],3)`, CI `r round(m4par["stock 1 Recruitment Ecov: LarvalTempDuration $\\beta_1$",3],3)`, `r round(m4par["stock 1 Recruitment Ecov: LarvalTempDuration $\\beta_1$",4],3)` # Discussion +The duration of optimal larval temperature covariate resulted in slightly better model fits and reduced recruitment variability relative to the model without covariates. The effect of optimal larval temperature on recruitment was postive, as hypothesized, with fewer days of optimal temperature resulting in a lower recruitment scaling parameter. However, the confidence interval of the effect included 0. + +While some of the zooplankton time series also resulted in slightly better model fits and reduced recruitment variability relative to the model without covariates, the effects of zooplankton were negative on recruitment. This relationship is opposite the hypothesized relationship between herring and food, which was expected to be positive. + + + # References diff --git a/docs/WP_RecCovariates_Gaichas.html b/docs/WP_RecCovariates_Gaichas.html new file mode 100644 index 0000000..f66089c --- /dev/null +++ b/docs/WP_RecCovariates_Gaichas.html @@ -0,0 +1,2308 @@ + + + + + + + + + + + + + + + +Working Paper: Recruitment covariate testing in WHAM + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + +
+

1 Introduction

+

The Atlantic herring research track assessment working group (WG) prioritized investigation of recruitment drivers as potential stock assessment model covariates, because low recruitment in recent years is an important issue for the stock and for fishery management.

+

The WG used a boosted regression tree analysis (Molina 2024) to identify zooplankton indices that best explained patterns in herring recruitment. These indices included large copepods in spring (influencing growth of herring postlarvae and juveniles), small copeopods in fall (influencing survival of herring larvae over the winter), haddock egg predation (influencing egg mortality), and temperature (influencing larval and juvenile survival).

+

In this working paper, we evaluate each of these indices as potential recruitment covariates in the Atlantic herring assessment implemented during the research track in the Woods Hole Assessment Model (WHAM, Stock and Miller (2021)).

+
+
+

2 Methods

+

We are using the devel version of WHAM: https://github.com/timjmiller/wham/tree/devel

+

Model mm192 is our starting point.

+

Haddock egg predation, Zooplankton, and temperature indices were explored as covariates on herring recruitment. Methods for developing each indicator are in separate working papers.

+

Recruitment is modeled as deviations from the “recruitment scaling parameter”, leaving one option for modeling effects of covariates on recruitment: “controlling”.

+

A “controlling” recruitment covariate results in a time-varying recruitment scaling parameter.

+

[merge with Jon’s ecisting text on the haddock egg predation testing here… ]

+

We explored indices with different zooplankton groups, seasons, and regions according to herring life history and results from the boosted regression tree:

+
    +
  • Jan-Jun (Spring) large copepods in spring herring BTS strata with lag-0 to represent food for pre-recruit juveniles
    • +
  • Jul-Dec (Fall) small copepods in fall herring BTS strata with lag-1 to represent food for larvae in general
    • +
  • Sep-Feb small copepods in herring larval area with lag-1 to represent food for larvae more specifically
    • +
  • Combinations of large and small copepod covariates above
    • +
+

In addition, we explored an index of optimal temperature duration (days) during the fall larval season, September-December

+

We evaluated

+
    +
  • Options for covariate input (millions of cells vs. log(cells), VAST estimated SE vs. WHAM estimated SE)
    +
    • +
  • Options for covariate observation model (“rw” vs. “ar1”)
    +
    • +
  • Options for recruitment link (“none” vs. “controlling-linear” with lag-0 for large copepods, lag-1 for small copepods, and lag-1 for larval temperature duration)
    +
    • +
  • No attempts were made to fit polynomial effects although this is possible in WHAM
    • +
+

The main diagnostics we used to determine if the model was improved by covariates were:

+
    +
  • model converged and Hessian matrix invertable
    • +
  • dAIC lowest
    • +
  • recruitment sigma reduced (how much?)
    • +
  • estimated covariate effect CI does not include 0
    • +
  • direction of covariate effect is sensible
    • +
+

Sensitivity runs

+

Turn off NAA RE and then try fitting with the most robust recruitment covariates.

+
+
+

3 Results

+

Short story:

+
    +

  • Models with covariates input on the log scale generally converged

    • +

  • Models with WHAM estimated covariate SE (“est_1”) generally converged

    • +

  • Under the above conditions, most models with and without the recruitment link converged

    • +

  • Models with the Jan-Jun (Spring) large copepods covariate also converged with input as millions of cells and VAST estimated SE

    • +
+

Way too much detail including false starts

+

Summary table with each model, recruitment sigma compared with base, covariate beta with CI

+

We show diagnostics and results from the best fit models for each covariate in sections below.

+

Diagnostics include:

+
    +
  • WHAM’s fit to the covariate time series
    • +
  • One step ahead residuals for WHAM’s fit
    • +
  • Estimated time varying recruitment scaling parameter with N age 1
    • +
+
+

3.1 Spring large copepods

+

Models with no covariates had slightly better AIC than comparable models with recruitment links

+
config <- "lgcopeSpring2"
+
+# name the model directory
+name <- paste0("mm192_", config)
+
+write.dir <- here::here(sprintf("WHAMfits/%s/", name))
+
+# larger set of ecov setups to compare
+df.mods <- data.frame(Recruitment = c(rep(2, 16)),
+                      ecov_process = c(rep("rw",8),rep("ar1",8)),
+                      ecov_how = c(rep("none",4), 
+                                     rep("controlling-lag-0-linear",4)),
+                      ecovdat = rep(c("logmean-logsig", 
+                                      "logmean-est_1",
+                                      "meanmil-logsigmil",
+                                      "meanmil-est_1"),4),
+                      stringsAsFactors=FALSE)
+n.mods <- dim(df.mods)[1]
+df.mods$Model <- paste0("m",1:n.mods)
+df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
+
+
+mod.list <- paste0(write.dir, df.mods$Model,".rds")
+mods <- lapply(mod.list, readRDS)
+
+n.mods <- length(mod.list)
+
+vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
+for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
+
+
+opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
+ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
+df.mods$conv <- as.logical(opt_conv)
+df.mods$pdHess <- as.logical(ok_sdrep)
+df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
+
+not_conv <- !df.mods$conv | !df.mods$pdHess
+mods2 <- mods
+mods2[not_conv] <- NULL
+df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
+df.aic <- df.aic.tmp[FALSE,]
+ct = 1
+for(i in 1:n.mods){
+  if(not_conv[i]){
+    df.aic[i,] <- rep(NA,5)
+  } else {
+    df.aic[i,] <- df.aic.tmp[ct,]
+    ct <- ct + 1
+  }
+}
+df.mods <- cbind(df.mods, df.aic)
+df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
+df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
+rownames(df.mods) <- NULL
+
+top4 <- df.mods |>
+  dplyr::filter(ecovdat %in% ("logmean-est_1")) |>
+  dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC)
+
flextable::flextable(top4) |>
+  flextable::set_header_labels(ecov_process = "Process",
+                               ecov_how = "Rec. link") |>
+  flextable::set_table_properties(layout = "autofit")
+

Model

Process

Rec. link

NLL

conv

pdHess

dAIC

AIC

m10

ar1

none

-1,793.611

TRUE

TRUE

0

-3325.2

m14

ar1

controlling-lag-0-linear

-1,794.325

TRUE

TRUE

0.6

-3324.6

m2

rw

none

-1,790.837

TRUE

TRUE

3.5

-3321.7

m6

rw

controlling-lag-0-linear

-1,791.227

TRUE

TRUE

4.7

-3320.5

+
+

3.1.1 Spring large copepods covariate: logscale ar1 diagnostics

+

The WHAM ar1 fit looks strange

+
+lgCopeSp2 m10 fit +
lgCopeSp2 m10 fit
+
+
+lgCopeSp2 m10 osa +
lgCopeSp2 m10 osa
+
+
+
+

3.1.2 Spring large copepods covariate: logscale ar1 recruitment

+
+lgCopeSp2 m10 rec +
lgCopeSp2 m10 rec
+
+
+lgCopeSp2 m14 rec +
lgCopeSp2 m14 rec
+
+
m10par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m10/res_tables/parameter_estimates_table.RDS"))
+
+m14par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m14/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.797; lgCopeSpring2 beta_1 is -0.407, CI -1.063, 0.25

+
+
+

3.1.3 Spring large copepods covariate: logscale rw diagnostics

+
+lgCopeSp2 m2 fit +
lgCopeSp2 m2 fit
+
+
+lgCopeSp2 m2 osa +
lgCopeSp2 m2 osa
+
+
+
+

3.1.4 Spring large copepods covariate: logscale rw recruitment

+
+lgCopeSp2 m10 rec +
lgCopeSp2 m10 rec
+
+
+lgCopeSp2 m14 rec +
lgCopeSp2 m14 rec
+
+
m2par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m2/res_tables/parameter_estimates_table.RDS"))
+
+m6par <- readRDS(here::here("WHAMfits/mm192_lgcopeSpring2/m6/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.804; lgCopeSpring2 beta_1 is -0.45, CI -1.438, 0.538

+
+
+
+

3.2 Fall small copepods

+

Models with no covariates had slightly better AIC than the ar1 model with recruitment links

+
config <- "smcopeFall2"
+
+# name the model directory
+name <- paste0("mm192_", config)
+
+write.dir <- here::here(sprintf("WHAMfits/%s/", name))
+
+# larger set of ecov setups to compare
+df.mods <- data.frame(Recruitment = c(rep(2, 16)),
+                      ecov_process = c(rep("rw",8),rep("ar1",8)),
+                      ecov_how = c(rep("none",4), 
+                                     rep("controlling-lag-0-linear",4)),
+                      ecovdat = rep(c("logmean-logsig", 
+                                      "logmean-est_1",
+                                      "meanmil-logsigmil",
+                                      "meanmil-est_1"),4),
+                      stringsAsFactors=FALSE)
+n.mods <- dim(df.mods)[1]
+df.mods$Model <- paste0("m",1:n.mods)
+df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
+
+
+mod.list <- paste0(write.dir, df.mods$Model,".rds")
+mods <- lapply(mod.list, readRDS)
+
+n.mods <- length(mod.list)
+
+vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
+for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
+
+
+opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
+ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
+df.mods$conv <- as.logical(opt_conv)
+df.mods$pdHess <- as.logical(ok_sdrep)
+df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
+
+not_conv <- !df.mods$conv | !df.mods$pdHess
+mods2 <- mods
+mods2[not_conv] <- NULL
+df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
+df.aic <- df.aic.tmp[FALSE,]
+ct = 1
+for(i in 1:n.mods){
+  if(not_conv[i]){
+    df.aic[i,] <- rep(NA,5)
+  } else {
+    df.aic[i,] <- df.aic.tmp[ct,]
+    ct <- ct + 1
+  }
+}
+df.mods <- cbind(df.mods, df.aic)
+df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
+df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
+rownames(df.mods) <- NULL
+
+top4 <- df.mods |>
+  dplyr::filter(ecovdat %in% ("logmean-est_1")) |>
+  dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC)
+
flextable::flextable(top4) |>
+  flextable::set_header_labels(ecov_process = "Process",
+                               ecov_how = "Rec. link") |>
+  flextable::set_table_properties(layout = "autofit")
+

Model

Process

Rec. link

NLL

conv

pdHess

dAIC

AIC

m10

ar1

none

-1,797.177

TRUE

TRUE

0

-3332.4

m2

rw

none

-1,796.175

TRUE

TRUE

0.1

-3332.3

m14

ar1

controlling-lag-0-linear

-1,797.931

TRUE

TRUE

0.5

-3331.9

m6

rw

controlling-lag-0-linear

-1,794.370

TRUE

FALSE

---

---

+
+

3.2.1 Fall small copepods covariate: logscale ar1 diagnostics

+
+smCopeFall2 m10 fit +
smCopeFall2 m10 fit
+
+
+smCopeFall2 m10 osa +
smCopeFall2 m10 osa
+
+
+
+

3.2.2 Fall small copepods covariate: logscale ar1 recruitment

+
+smCopeFall2 m10 rec +
smCopeFall2 m10 rec
+
+
+smCopeFall2 m14 rec +
smCopeFall2 m14 rec
+
+
m10par <- readRDS(here::here("WHAMfits/mm192_smcopeFall2/m10/res_tables/parameter_estimates_table.RDS"))
+
+m14par <- readRDS(here::here("WHAMfits/mm192_smcopeFall2/m14/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.79; smcopeFall2 beta_1 is -1.013, CI -2.715, 0.689

+
+
+
+

3.3 September - February small copepods in the larval herring area

+

Model with ar1 covariate had slightly better AIC than models without recruitment links

+
config <- "smcopeSepFeb2"
+
+# name the model directory
+name <- paste0("mm192_", config)
+
+write.dir <- here::here(sprintf("WHAMfits/%s/", name))
+
+# larger set of ecov setups to compare
+df.mods <- data.frame(Recruitment = c(rep(2, 16)),
+                      ecov_process = c(rep("rw",8),rep("ar1",8)),
+                      ecov_how = c(rep("none",4), 
+                                   rep("controlling-lag-1-linear",4)),
+                      ecovdat = rep(c("logmean-logsig", 
+                                      "logmean-est_1",
+                                      "meanmil-logsigmil",
+                                      "meanmil-est_1"),4),
+                      stringsAsFactors=FALSE)
+n.mods <- dim(df.mods)[1]
+df.mods$Model <- paste0("m",1:n.mods)
+df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
+
+
+mod.list <- paste0(write.dir, df.mods$Model,".rds")
+mod.list <- mod.list[-16] # mod 16 didnt run
+mods <- lapply(mod.list, readRDS) 
+
+df.mods <- df.mods[-16,] # mod 16 didnt run
+
+n.mods <- length(mod.list)
+
+vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
+for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
+
+
+opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
+ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
+df.mods$conv <- as.logical(opt_conv)
+df.mods$pdHess <- as.logical(ok_sdrep)
+df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
+
+not_conv <- !df.mods$conv | !df.mods$pdHess
+mods2 <- mods
+mods2[not_conv] <- NULL
+df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
+df.aic <- df.aic.tmp[FALSE,]
+ct = 1
+for(i in 1:n.mods){
+  if(not_conv[i]){
+    df.aic[i,] <- rep(NA,5)
+  } else {
+    df.aic[i,] <- df.aic.tmp[ct,]
+    ct <- ct + 1
+  }
+}
+df.mods <- cbind(df.mods, df.aic)
+df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
+df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
+rownames(df.mods) <- NULL
+
+top4 <- df.mods |>
+  dplyr::filter(ecovdat %in% ("logmean-est_1")) |>
+  dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC)
+
flextable::flextable(top4) |>
+  flextable::set_header_labels(ecov_process = "Process",
+                               ecov_how = "Rec. link") |>
+  flextable::set_table_properties(layout = "autofit")
+

Model

Process

Rec. link

NLL

conv

pdHess

dAIC

AIC

m14

ar1

controlling-lag-1-linear

-1,791.900

TRUE

TRUE

0

-3319.8

m2

rw

none

-1,789.657

TRUE

TRUE

0.5

-3319.3

m10

ar1

none

-1,790.469

TRUE

TRUE

0.9

-3318.9

m6

rw

controlling-lag-1-linear

-1,785.746

TRUE

FALSE

---

---

+
+

3.3.1 Sep-Feb small copepods in herring larval area covariate: logscale ar1 diagnostics

+
+smCopeSepFeb2 m10 fit +
smCopeSepFeb2 m10 fit
+
+
+smCopeSepFeb2 m10 osa +
smCopeSepFeb2 m10 osa
+
+
+
+

3.3.2 Sep-Feb small copepods in herring larval area covariate: logscale ar1 recruitment

+
+smCopeSepFeb2 m10 rec +
smCopeSepFeb2 m10 rec
+
+
+smCopeSepFeb2 m14 rec +
smCopeSepFeb2 m14 rec
+
+
m10par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2/m10/res_tables/parameter_estimates_table.RDS"))
+
+m14par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2/m14/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.77; smcopeSepFeb2 beta_1 is -0.85, CI -1.883, 0.182

+
+
+
+

3.4 Alternate: Sep-Feb small copepods in fall herring survey strata covariate

+

Both rw and ar1 models worked with covariates, rw better fit

+
config <- "smcopeSepFeb2_fallstrata"
+
+# name the model directory
+name <- paste0("mm192_", config)
+
+write.dir <- here::here(sprintf("WHAMfits/%s/", name))
+
+# larger set of ecov setups to compare
+df.mods <- data.frame(Recruitment = c(rep(2, 16)),
+                      ecov_process = c(rep("rw",8),rep("ar1",8)),
+                      ecov_how = c(rep("none",4), 
+                                   rep("controlling-lag-1-linear",4)),
+                      ecovdat = rep(c("logmean-logsig", 
+                                      "logmean-est_1",
+                                      "meanmil-logsigmil",
+                                      "meanmil-est_1"),4),
+                      stringsAsFactors=FALSE)
+n.mods <- dim(df.mods)[1]
+df.mods$Model <- paste0("m",1:n.mods)
+df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
+
+
+mod.list <- paste0(write.dir, df.mods$Model,".rds")
+mod.list <- mod.list[-16] # mod 16 didnt run
+mods <- lapply(mod.list, readRDS) 
+
+df.mods <- df.mods[-16,] # mod 16 didnt run
+
+n.mods <- length(mod.list)
+
+vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
+for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
+
+
+opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
+ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
+df.mods$conv <- as.logical(opt_conv)
+df.mods$pdHess <- as.logical(ok_sdrep)
+df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
+
+not_conv <- !df.mods$conv | !df.mods$pdHess
+mods2 <- mods
+mods2[not_conv] <- NULL
+df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
+df.aic <- df.aic.tmp[FALSE,]
+ct = 1
+for(i in 1:n.mods){
+  if(not_conv[i]){
+    df.aic[i,] <- rep(NA,5)
+  } else {
+    df.aic[i,] <- df.aic.tmp[ct,]
+    ct <- ct + 1
+  }
+}
+df.mods <- cbind(df.mods, df.aic)
+df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
+df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
+rownames(df.mods) <- NULL
+
+top4 <- df.mods |>
+  dplyr::filter(ecovdat %in% ("logmean-est_1")) |>
+  dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC)
+
flextable::flextable(top4) |>
+  flextable::set_header_labels(ecov_process = "Process",
+                               ecov_how = "Rec. link") |>
+  flextable::set_table_properties(layout = "autofit")
+

Model

Process

Rec. link

NLL

conv

pdHess

dAIC

AIC

m6

rw

controlling-lag-1-linear

-1,791.981

TRUE

TRUE

0

-3322

m14

ar1

controlling-lag-1-linear

-1,792.464

TRUE

TRUE

1.1

-3320.9

m2

rw

none

-1,790.196

TRUE

TRUE

1.6

-3320.4

m10

ar1

none

-1,790.734

TRUE

TRUE

2.5

-3319.5

+
+

3.4.1 Sep-Feb small copepods in herring fall survey area covariate: logscale rw diagnostics

+
+smCopeSepFeb2_fallstrat m2 fit +
smCopeSepFeb2_fallstrat m2 fit
+
+
+smCopeSepFeb2_fallstrat m2 osa +
smCopeSepFeb2_fallstrat m2 osa
+
+
+
+

3.4.2 Sep-Feb small copepods in herring fall survey area covariate: logscale rw recruitment

+
+smCopeSepFeb2_fallstrat m2 rec +
smCopeSepFeb2_fallstrat m2 rec
+
+
+smCopeSepFeb2_fallstrat m6 rec +
smCopeSepFeb2_fallstrat m6 rec
+
+
m2par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2_fallstrata/m2/res_tables/parameter_estimates_table.RDS"))
+
+m6par <- readRDS(here::here("WHAMfits/mm192_smcopeSepFeb2_fallstrata/m6/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.762; smcopeSepFeb2 beta_1 is -1.01, CI -2.11, 0.089

+
+
+
+

3.5 September - December duration of larval optimal temperature

+

Both rw and ar1 models worked with covariates, rw better fit

+
config <- "LarvalTempDuration"
+
+# name the model directory
+name <- paste0("mm192_", config)
+
+write.dir <- here::here(sprintf("WHAMfits/%s/", name))
+
+# larger set of ecov setups to compare
+# larger set of ecov setups to compare
+df.mods <- data.frame(Recruitment = c(rep(2, 8)),
+                      ecov_process = c(rep("rw",4),rep("ar1",4)),
+                      ecov_how = rep(c("none","controlling-lag-1-linear"), 4),
+                      ecovdat = c(rep("mean-est_1", 2),rep("logmean-est_1",2)),
+                      stringsAsFactors=FALSE)
+n.mods <- dim(df.mods)[1]
+df.mods$Model <- paste0("m",1:n.mods)
+df.mods <- dplyr::select(df.mods, Model, tidyselect::everything()) # moves Model to first col
+
+
+mod.list <- paste0(write.dir, df.mods$Model,".rds")
+mods <- lapply(mod.list, readRDS)
+
+n.mods <- length(mod.list)
+
+vign2_conv <- lapply(mods, function(x) capture.output(check_convergence(x)))
+for(m in 1:n.mods) cat(paste0("Model ",m,":"), vign2_conv[[m]], "", sep='\n')
+
+
+opt_conv = 1-sapply(mods, function(x) x$opt$convergence)
+ok_sdrep = sapply(mods, function(x) if(x$na_sdrep==FALSE & !is.na(x$na_sdrep)) 1 else 0)
+df.mods$conv <- as.logical(opt_conv)
+df.mods$pdHess <- as.logical(ok_sdrep)
+df.mods$NLL <- sapply(mods, function(x) round(x$opt$objective,3))
+
+not_conv <- !df.mods$conv | !df.mods$pdHess
+mods2 <- mods
+mods2[not_conv] <- NULL
+df.aic.tmp <- as.data.frame(compare_wham_models(mods2, table.opts=list(sort=FALSE, calc.rho=T))$tab)
+df.aic <- df.aic.tmp[FALSE,]
+ct = 1
+for(i in 1:n.mods){
+  if(not_conv[i]){
+    df.aic[i,] <- rep(NA,5)
+  } else {
+    df.aic[i,] <- df.aic.tmp[ct,]
+    ct <- ct + 1
+  }
+}
+df.mods <- cbind(df.mods, df.aic)
+df.mods <- df.mods[order(df.mods$dAIC, na.last=TRUE),]
+df.mods[is.na(df.mods$AIC), c('dAIC','AIC','rho_R','rho_SSB','rho_Fbar')] <- "---"
+rownames(df.mods) <- NULL
+
+top4 <- df.mods |>
+  dplyr::filter(ecovdat %in% ("logmean-est_1")) |>
+  dplyr::select(Model, ecov_process, ecov_how, NLL, conv, pdHess, dAIC, AIC)
+
flextable::flextable(top4) |>
+  flextable::set_header_labels(ecov_process = "Process",
+                               ecov_how = "Rec. link") |>
+  flextable::set_table_properties(layout = "autofit")
+

Model

Process

Rec. link

NLL

conv

pdHess

dAIC

AIC

m4

rw

controlling-lag-1-linear

-1,827.543

TRUE

TRUE

0

-3393.1

m3

rw

none

-1,826.370

TRUE

TRUE

0.4

-3392.7

m8

ar1

controlling-lag-1-linear

-1,826.627

TRUE

TRUE

3.8

-3389.3

m7

ar1

none

-1,825.511

TRUE

TRUE

4.1

-3389

+
+

3.5.1 Duration of optimal larval temperature, Sept-Dec: logscale rw diagnostics

+
+LarvalTempDuration/m3 fit +
LarvalTempDuration/m3 fit
+
+
+LarvalTempDuration/m3 osa +
LarvalTempDuration/m3 osa
+
+
+
+

3.5.2 Duration of optimal larval temperature, Sept-Dec: logscale rw recruitment

+
+LarvalTempDuration/m3 rec +
LarvalTempDuration/m3 rec
+
+
+LarvalTempDuration/m4 rec +
LarvalTempDuration/m4 rec
+
+
m3par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration/m3/res_tables/parameter_estimates_table.RDS"))
+
+m4par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration/m4/res_tables/parameter_estimates_table.RDS"))
+

Without covariate, recruitment variance is 0.823, and with is 0.79; LarvalTempDuration beta_1 is 2.066, CI -0.657, 4.789

+
+
+
+
+

4 Discussion

+

The duration of optimal larval temperature covariate resulted in slightly better model fits and reduced recruitment variability relative to the model without covariates. The effect of optimal larval temperature on recruitment was postive, as hypothesized, with fewer days of optimal temperature resulting in a lower recruitment scaling parameter. However, the confidence interval of the effect included 0.

+

While some of the zooplankton time series also resulted in slightly better model fits and reduced recruitment variability relative to the model without covariates, the effects of zooplankton were negative on recruitment. This relationship is opposite the hypothesized relationship between herring and food, which was expected to be positive.

+
+
+

References

+
+
+Stock, B.C., and Miller, T.J. 2021. The Woods Hole Assessment Model (WHAM): A general state-space assessment framework that incorporates time- and age-varying processes via random effects and links to environmental covariates. Fisheries Research 240: 105967. doi:10.1016/j.fishres.2021.105967. +
+
+
+ + + +
+
+ +
+ + + + + + + + + + + + + + + + + diff --git a/docs/index.Rmd b/docs/index.Rmd index ae15a82..91c7c1c 100644 --- a/docs/index.Rmd +++ b/docs/index.Rmd @@ -33,3 +33,6 @@ output: html_document Draft zooplankton index WP - Nov 2024 + Draft zoop and temp covariates WP - Nov 2024 + + diff --git a/docs/index.html b/docs/index.html index 8ddf0ec..824af70 100644 --- a/docs/index.html +++ b/docs/index.html @@ -383,6 +383,8 @@

Analysis pages

Code to make the SOE indicator datasets - Oct 2024

Draft zooplankton index WP - Nov 2024

+

+Draft zoop and temp covariates WP - Nov 2024