WP_RecCovariates_Gaichas.Rmd

---
title: "Working Paper: Recruitment covariate testing in WHAM"
author: "Sarah Gaichas and Jon Deroba"
date: "`r Sys.Date()`"
output:
  bookdown::pdf_document2:
    includes: 
       in_header: latex/header1.tex
    keep_tex: true
  bookdown::html_document2: 
    toc: true
    toc_float: true
    code_fold: hide
  bookdown::word_document2:
    toc: true
link-citations: yes
csl: "canadian-journal-of-fisheries-and-aquatic-sciences.csl"
bibliography: zoopindex.bib
urlcolor: blue
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE,
                      message = FALSE,
                      warning = FALSE)

library(tidyverse)
library(here)
library(DT)
library(patchwork)
library(wham)

```

# 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_woods_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. 


# Methods

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

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


# 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](https://noaa-edab.github.io/zooplanktonindex/WHAMcovariate_tests.html)

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

## 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)`


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

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

```{r, results='hide'}

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)

```

```{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_NAA_REoff/m3/plots_png/diagnostics/LarvalTempDuration_diagnostic.png)


![LarvalTempDuration/m3 osa](WHAMfits/mm192_LarvalTempDuration_NAA_REoff/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_NAA_REoff/m3/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png)

![LarvalTempDuration/m4 rec](WHAMfits/mm192_LarvalTempDuration_NAA_REoff/m4/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png)


```{r}
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 `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. 

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

# References