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 + + + + + + + + +### Spring large copepods covariate: logscale ar1 recruitment + + + + + + + +```{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 + + + + + + + + + +### Spring large copepods covariate: logscale rw recruitment + + + + + + + + + +```{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 + + + + + + + +### Fall small copepods covariate: logscale ar1 recruitment + + + + + + + +```{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 + + + + + + + + +### Sep-Feb small copepods in herring larval area covariate: logscale ar1 recruitment + + + + + + + + +```{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 + + + + + + +### Sep-Feb small copepods in herring fall survey area covariate: logscale rw recruitment + + + + + + + +```{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 + + + + + + + + + +### Duration of optimal larval temperature, Sept-Dec: logscale rw recruitment + + + + + + + +```{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 @@ + + + + +
+ + + + + + + + + + +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)).
+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:
+In addition, we explored an index of optimal temperature duration (days) during the fall larval season, September-December
+We evaluated
+The main diagnostics we used to determine if the model was improved by covariates were:
+Sensitivity runs
+Turn off NAA RE and then try fitting with the most robust recruitment covariates.
+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:
+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 |
The WHAM ar1 fit looks strange
+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
+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
+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 | --- | --- |
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
+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 | --- | --- |
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
+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 |
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
+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 |
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
+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.
+