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

As a sensitivity test, we also evaluate the effects of potential recruitment covariates with NAA random effects turned off in the Atlantic herring assessment model.

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 osa

3.1.2 Spring large copepods covariate: logscale ar1 recruitment

lgCopeSp2 m10 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 osa

3.1.4 Spring large copepods covariate: logscale rw recruitment

lgCopeSp2 m10 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 osa

3.2.2 Fall small copepods covariate: logscale ar1 recruitment

smCopeFall2 m10 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 osa

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

smCopeSepFeb2 m10 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 osa

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

smCopeSepFeb2_fallstrat m2 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 osa

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

LarvalTempDuration/m3 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

3.6 Sensitivity: Remove NAA RE, add September - December duration of larval optimal temperature

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

config <- "LarvalTempDuration_NAA_REoff"

# 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,762.281	TRUE	TRUE	0	-3266.6
m8	ar1	controlling-lag-1-linear	-1,761.214	TRUE	TRUE	4.2	-3262.4
m3	rw	none	-1,748.086	TRUE	TRUE	26.4	-3240.2
m7	ar1	none	-1,747.227	TRUE	TRUE	30.1	-3236.5

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

LarvalTempDuration/m3 fit

LarvalTempDuration/m3 osa

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

LarvalTempDuration/m3 rec

LarvalTempDuration/m4 rec

m3par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration_NAA_REoff/m3/res_tables/parameter_estimates_table.RDS"))

m4par <- readRDS(here::here("WHAMfits/mm192_LarvalTempDuration_NAA_REoff/m4/res_tables/parameter_estimates_table.RDS"))

Without covariate, recruitment variance is 0.833, and with is 0.503; LarvalTempDuration beta_1 is 5.122, CI 2.709, 7.535

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.

The sensitivity run removing NAA RE resulted in a much stronger impact of the optimal larval temperature on recruitment.

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.