Attempt recruitment covariates in herring WHAM model
Sarah Gaichas
2024-10-04
Introduction
Now that we have way too many copeopod time series, lets see if any of them are useful covariates in the Atlantic herring assessment model.
The primary target is a covariate on Atlantic herring recruitment.
Adelle’s boosted regression tree model shows both spring large copepods and fall small copepods as important factors explaining recruitment patterns estimated by the WHAM model.
We have also developed a herring-larvae specific small copepod index within the area with more than the 70th percentile of cumulative 1982-2022 herring larvae density.
So we can try at least three flavors of copepod index as a model covariate.
The code that takes the output of the VAST models and makes WHAM inputs is WHAMinputs.R
Also adding a zooplankton volume spring index to try in the place of the spring large copepods.
Methods
Installed multiWHAM from the devel branch: https://github.com/timjmiller/wham/tree/devel
devtools::install_github("timjmiller/wham", dependencies=TRUE, ref="devel")
Installation worked on the outdated mac os… so far.
Lets test the installation by seeing if I can get the same outputs as Jon did by running the model without any additions.
Modify Jon’s code from the Atlantic Herring RTWG google drive and now
in this directory, Example WHAM code to share.R which
sources his WHAMfxns.R and hindcast.R also now
in this directory.
First run the base model and compare with the stored to ensure my installation is ok:
#source some functions
source(here::here("WHAMfits/WHAMfxns.R"))
#source the hindcast function for MASE etc.
source(here::here("WHAMfits/hindcast.R"))
# starting with mm192, add suffix for different attempts
# test with no change
config <- "nochange"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
# copy the input dat file to model directory
asapRun2_87 <- here::here("WHAMfits/mm192/Run2_87.DAT")
file.copy(from=file.path(asapRun2_87), to=write.dir, overwrite=FALSE)
asapRun2dat<- read_asap3_dat(here::here(write.dir, "Run2_87.DAT"))
# need data files read in to run these--but don't,
# the full data object is already available in the mm192 rds file
#asapRun2dat=CVfxn(asapRun9dat) #empirical CVs
#asapRun2dat=ESSfxn(asapRun9dat,addnumcatch=150,addnumsurv=150)
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
m1 <- fit_wham(input, do.osa = F)
saveRDS(m1, file.path(write.dir, "mm192nochange.rds"))
check_convergence(m1)
# these should be identical; they are
test <- compare_wham_models(mods = list("orig" = mm192mod, "nochange" = m1), fdir = write.dir)
saveRDS(test, file.path(write.dir, "compare.rds"))The original model and my rerun using the same inputs (nochange) directly overlap:
knitr::include_graphics(here::here("WHAMfits/mm192_nochange/compare_png/compare_SSB_F_R.png"))
New explore environmental covariate (ecov) options. What format does the input need?
Following https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html columns Year, Index, Index_sigma
We are testing indices of bottom-up recruitment impacts–food for larvae, postlarvae, and juveniles
Format ecov indices; we’ll try at least 3 variants:
- Spring large copepods (food for postlarvae and juveniles, series tested in boosted regression tree)
- Fall small copepods (food for larvae, series tested in boosted regression tree)
- Sep-Feb small copepods in herring larval area (food for larvae, new series not tested in boosted regression tree but developed to better match timing and area of herring larvae)
The code in WHAMinputs.R was used to create csv files of the VAST derives zooplankton indices. Files are saved in the WHAMfits folder.
Spring large copepods covariate
Start with spring large copepods. Following the suggestions in WHAM vignette 2:
OK thats outmoded. Can just add an ecov object to existing input
using set_ecov:
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
# make it the same length as model ?
#env.dat <- env.dat[env.dat$Time > 1986,]
# test a single ecov setup with no effect first
ecov_on <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)), #"est_1",
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = "ar1", # "rw" or "ar1"
#where = "recruit",
recruitment_how = as.matrix("controlling-lag-0-linear")) # 0 = no effect, 1 = controlling, 2 = limiting
# modify current model input call adding ecov list
#input$call
ecovinput <- set_ecov(input, ecov=ecov_on)
# ecovinput <- prepare_wham_input(asap3 = asapdat, # correct structurr, asap3 is wrong structure
# model_name = "test", #df.mods[m, "Model"],
# ecov = ecov,
# recruit_model = 2,
# selectivity = list(fix_pars = list(7:8,
# 2,
# c(1:2, 4:8),
# 1:8,
# c(1, 3:8),
# NULL,
# c(1, 4:8),
# NULL)),
# M = list(re_model = as.matrix("none")), #as.matrix(df.mods[m, "M_re"])),
# NAA_re = list(sigma = "rec+1",
# cor = "ar1_a", #df.mods[m, "naa_cor"],
# N1_model = "age-specific-fe", #init_Nnum[m],
# decouple_recruitment = TRUE),
# age_comp = "logistic-normal-ar1-miss0",
# basic_info = list(bias_correct_process = FALSE,
# bias_correct_BRPs = FALSE,
# percentSPR = 40,
# percentFXSPR = 100,
# XSPR_R_avg_yrs = 6:35,
# XSPR_R_opt = 5))
ecovon_mod <- fit_wham(ecovinput, do.osa = F)
saveRDS(ecovon_mod, file.path(write.dir, "testecovon.rds"))
plot_wham_output(mod = ecovon_mod, dir.main = file.path(write.dir))
# test a single ecov setup with no effect first
ecov_off <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = "ar1", # "rw" or "ar1"
#where = "recruit",
recruitment_how = as.matrix("none")) # 0 = no effect, 1 = controlling, 2 = limiting
# modify current model input call adding ecov list
#input$call
ecovinput <- set_ecov(input, ecov=ecov_off)
mod <- fit_wham(ecovinput, do.osa = F)
saveRDS(mod, file.path(write.dir, "testecovoff.rds"))
#compare to original, not the same data
test <- compare_wham_models(mods = list("testecovoff" = mod, "testecovon" = ecovon_mod), fdir = write.dir)
saveRDS(test, file.path(write.dir, "testcompare.rds"))Test run was broken until I changed units of zooplankton to millions of cells.
Now everything runs, the covariate is included, and makes no difference in the model (which looks identical to the base model so thats a plus):
testcompare <- readRDS(here::here("WHAMfits/mm192_largecope/testcompare.rds"))
knitr::include_graphics(here::here("WHAMfits/mm192_largecope/plots_png/results/Ecov_1_lgCopeSpr.png"))
testcompare$tab## dAIC AIC rho_R rho_SSB rho_Fbar
## testecovoff 0.0 -2759.1 0.8682 0.5901 -0.2314
## testecovon 0.4 -2758.7 0.9081 0.5920 -0.2379testcompare$g[[1]]
Lets try a suite of models with the same recruitment process but different ways the large copepod series relates to it.
We’ll always use recruitment process 2, random about mean, because that is in mm192.
We’ll look at both random walk “rw” and AR1 “ar1” processes for the environmental covariate.
We’ll set the parameterrecruitment_how as one of the
following:
{= “none”}{no effect.}
{= “controlling”}{pre-recruit density-independent mortality.}
{= “limiting”}{ maximum recruitment, e.g. ecov determines amount of
suitable habitat)}
{= “lethal”}{threshold, i.e. R –> 0 at some ecov value.}
{= “masking”}{metabolic/growth, decreases dR/dS}
{= “directive”}{e.g. behavioral}
Going to try “none”, “controlling”, “limiting”, and “masking” (note that “masking” didn’t work well in the state space RTA simulations).
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),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
for(m in 5:n.mods){
ecov <- list(
label = "lgCopeSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Check all for convergence (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
config <- "largecope"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
# need this for html to work because block above running model not evaluated on knit
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
recruitment_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.31e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.47e-03
## Max gradient parameter: logit_q
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.96e-03
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.38e+07
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.33e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e-03
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this model
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-02
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.62e-03
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this modelCompare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.7 -2755.4 0.8594 0.5907 -0.2362
## m2 4.3 -2754.8 0.8104 0.5940 -0.2448
## m3 4.3 -2754.8 0.8162 0.6011 -0.2446
## m4 0.0 -2759.1 0.8682 0.5901 -0.2314
## m5 0.4 -2758.7 0.9081 0.5920 -0.2379
## m6 0.4 -2758.7 0.5715 0.5914 -0.2399
## m7 1.9 -2757.2 0.8908 0.5945 -0.2374df.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
df.mods## Model Recruitment ecov_process recruitment_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m6 2 ar1 controlling-lag-0-linear TRUE TRUE
## 3 m7 3 ar1 limiting-lag-0-linear TRUE TRUE
## 4 m8 3 ar1 masking-lag-0-linear TRUE TRUE
## 5 m1 2 rw none TRUE TRUE
## 6 m2 2 rw controlling-lag-0-linear TRUE TRUE
## 7 m3 3 rw limiting-lag-0-linear TRUE TRUE
## 8 m4 3 rw masking-lag-0-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1509.570 0 -2759.1 0.8682 0.5901 -0.2314
## 2 -1510.349 0.4 -2758.7 0.9081 0.592 -0.2379
## 3 -1510.349 0.4 -2758.7 0.5715 0.5914 -0.2399
## 4 -1509.592 1.9 -2757.2 0.8908 0.5945 -0.2374
## 5 -1506.709 3.7 -2755.4 0.8594 0.5907 -0.2362
## 6 -1507.424 4.3 -2754.8 0.8104 0.594 -0.2448
## 7 -1507.424 4.3 -2754.8 0.8162 0.6011 -0.2446
## 8 -1175.320 --- --- --- --- ---Outputs are nearly identical for SSB and F, but there are small differences in status determination.
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
The numbers in the plots are not the model numbers! They are in order by the first table above
lgcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(lgcope1, file.path(write.dir, "lgcope1compare.rds"))
lgcope1$g[[1]]
lgcope1$g[[8]]
lgcope1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Fall small copepods covariate
Similar approach for first pass at fall small copepods. These will lag 1.
Setup
config <- "smcopeFall"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-1-linear",
"limiting-lag-1-linear",
"masking-lag-1-linear"),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 colRun models
env.dat <- read.csv(here::here("WHAMfits/fallsmallcopeALLindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
for(m in 1:n.mods){
ecov <- list(
label = "msCopeFall",
mean = as.matrix(env.dat$her_fa_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_fa_SE/1000000)),
year = env.dat$Time,
use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Estimated small copepods all Fall covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_msCopeFall.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.20e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.65e-01
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.06e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.34e-08
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.39e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.34e+00
## Max gradient parameter: logit_selpars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.64e-03
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 1.15e+07
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 11.6 -2597.8 3.2038 0.6009 -0.2487
## m2 0.0 -2609.4 0.8624 0.5830 -0.2315df.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
df.mods## Model Recruitment ecov_process ecov_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m1 2 rw none TRUE TRUE
## 3 m2 2 rw controlling-lag-1-linear TRUE FALSE
## 4 m3 3 rw limiting-lag-1-linear TRUE FALSE
## 5 m4 3 rw masking-lag-1-linear TRUE FALSE
## 6 m6 2 ar1 controlling-lag-1-linear TRUE FALSE
## 7 m7 3 ar1 limiting-lag-1-linear TRUE FALSE
## 8 m8 3 ar1 masking-lag-1-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1434.708 0 -2609.4 0.8624 0.583 -0.2315
## 2 -1427.902 11.6 -2597.8 3.2038 0.6009 -0.2487
## 3 -1429.056 --- --- --- --- ---
## 4 -1429.056 --- --- --- --- ---
## 5 -1427.902 --- --- --- --- ---
## 6 -1426.200 --- --- --- --- ---
## 7 -1428.500 --- --- --- --- ---
## 8 -1431.880 --- --- --- --- ---smcope1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcope1, file.path(write.dir, "smcope1compare.rds"))
smcope1$g[[1]]
smcope1$g[[8]]
smcope1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Sept-Feb small copepods in larval area covariate
Similar approach for first pass at fall small copepods. These will lag 1.
Setup
config <- "smcopeSepFeb"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-1-linear",
"limiting-lag-1-linear",
"masking-lag-1-linear"),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 colRun models
env.dat <- read.csv(here::here("WHAMfits/sepfebsmallcopeALLlarvareaindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
for(m in 1:n.mods){
ecov <- list(
label = "smCopeSepFeb",
mean = as.matrix(env.dat$her_larv_Estimate/1000000),
logsigma = as.matrix(log(env.dat$her_larv_SE/1000000)),
year = env.dat$Time,
use_obs = use.obs, # matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input206 <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Estimated small copepods September-February in herring larval area covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_smCopeSepFeb.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-11
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.51e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.21e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.99e+00
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.46e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.91e-04
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.60e-09
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.9 -2627.3 0.8541 0.5879 -0.2311
## m2 0.0 -2631.2 0.8793 0.6035 -0.2472df.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
df.mods## Model Recruitment ecov_process ecov_how conv pdHess NLL
## 1 m5 2 ar1 none TRUE TRUE -1445.582
## 2 m1 2 rw none TRUE TRUE -1442.652
## 3 m2 2 rw controlling-lag-1-linear TRUE FALSE -1444.922
## 4 m3 3 rw limiting-lag-1-linear TRUE FALSE -1444.922
## 5 m4 3 rw masking-lag-1-linear TRUE FALSE -1444.215
## 6 m6 2 ar1 controlling-lag-1-linear TRUE FALSE -1447.832
## 7 m7 3 ar1 limiting-lag-1-linear TRUE FALSE -1439.374
## 8 m8 3 ar1 masking-lag-1-linear TRUE FALSE -1445.582
## dAIC AIC rho_R rho_SSB rho_Fbar
## 1 0 -2631.2 0.8793 0.6035 -0.2472
## 2 3.9 -2627.3 0.8541 0.5879 -0.2311
## 3 --- --- --- --- ---
## 4 --- --- --- --- ---
## 5 --- --- --- --- ---
## 6 --- --- --- --- ---
## 7 --- --- --- --- ---
## 8 --- --- --- --- ---smcopelarv1 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcopelarv1, file.path(write.dir, "smcopelarv1compare.rds"))
smcopelarv1$g[[1]]
smcopelarv1$g[[8]]
smcopelarv1$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Spring zooplankton volume covariate
Try the spring zooplankton with full zooplankton volume instead of just large copepods.
Setup
config <- "zoopvolSpring"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# larger set of ecov setups to compare
df.mods <- data.frame(Recruitment = c(rep(c(2,2,3,3), 2)),
ecov_process = c(rep("rw",4),rep("ar1",4)),
ecov_how = rep(c("none",
"controlling-lag-0-linear",
"limiting-lag-0-linear",
"masking-lag-0-linear"),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 colRun models
env.dat <- read.csv(here::here("WHAMfits/springzoopvolumeindex.csv"), header=T)
# units of zoop volume are ml so lets convert to liters?
for(m in 1:n.mods){
ecov <- list(
label = "zoopvolSpr",
mean = as.matrix(env.dat$her_sp_Estimate/1000),
logsigma = as.matrix(log(env.dat$her_sp_SE/1000)),
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
#lag = 0, # postlarval and juvenile food in year t affects recruitment in year t
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
#where = c("recruit")[as.logical(df.mods$ecov_how[m])+1],
recruitment_how = as.matrix(df.mods$ecov_how[m])
) # 0 = no effect, 1 = controlling, 2 = limiting
# this isn't working, not estimating a recruitment function
# for the best because it will break
# for reference, this code should be added to attempt
#input <- set_NAA(mm192$mm192$input, NAA_re=list(recruit_model=3,sigma="rec+1",cor="ar1_a",N1_model="age-specific-fe",decouple_recruitment=TRUE))
input$data$recruit_model <- df.mods$Recruitment[m]
input$options$basic_info$recruit_model <- df.mods$Recruitment[m]
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Estimated zooplankton volume covariate (“ar1”, m5):
knitr::include_graphics(here::here(write.dir, "m5/plots_png/results/Ecov_1_zoopvolSpr.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
Ignore results from the limiting and masking models as Bev-Holt parameters aren’t really being estimated.
Also nothing worked that had a recruitment effect
# only take those that are there, but doesnt work with rest of code
#mod.list <- paste0(list.dirs(write.dir, recursive = FALSE), ".rds")
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.40e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.20e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.41e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 4:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.12e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.14e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has NOT converged: mod$opt$convergence != 0
## Maximum gradient component: 5.52e+02
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 18.2 -2590.5 0.7097 0.7253 -0.3005
## m2 0.0 -2608.7 0.8652 0.5978 -0.2361df.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
df.mods## Model Recruitment ecov_process ecov_how conv pdHess
## 1 m5 2 ar1 none TRUE TRUE
## 2 m1 2 rw none TRUE TRUE
## 3 m2 2 rw controlling-lag-0-linear TRUE FALSE
## 4 m3 3 rw limiting-lag-0-linear TRUE FALSE
## 5 m4 3 rw masking-lag-0-linear FALSE FALSE
## 6 m6 2 ar1 controlling-lag-0-linear TRUE FALSE
## 7 m7 3 ar1 limiting-lag-0-linear FALSE FALSE
## 8 m8 3 ar1 masking-lag-0-linear FALSE FALSE
## NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 -1434.361 0 -2608.7 0.8652 0.5978 -0.2361
## 2 -1424.248 18.2 -2590.5 0.7097 0.7253 -0.3005
## 3 -1424.585 --- --- --- --- ---
## 4 -1424.585 --- --- --- --- ---
## 5 -1146.241 --- --- --- --- ---
## 6 -1434.730 --- --- --- --- ---
## 7 -1146.241 --- --- --- --- ---
## 8 -1146.241 --- --- --- --- ---zoopvol <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(zoopvol, file.path(write.dir, "zoopvolcompare.rds"))
zoopvol$g[[1]]
zoopvol$g[[8]]
zoopvol$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Alternative WHAM runs
- Try using log of zooplankton index and/or letting WHAM estimate the logsigma.
- Would different lags be appropriate if the food affects spawning adults rather than larvae? (accidentally done as lag1)
Spring large copepods covariate, herring spring bottom trawl survey strata
Setup
config <- "lgcopeSpring2"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# asapdat<- read_asap3_dat(here::here("WHAMfits/mm192/Run2_87.DAT"))
#
# asap3 <- input$asap3
# 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 colRun models with different data entry options. Don’t do it in the Rmd, do it with a script.
env.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
for(m in 1:n.mods){
ecovdat <- dplyr::case_when(df.mods$ecovdat[m] %in% c("logmean-logsig", "logmean-est_1") ~
as.matrix(log(env.dat$her_sp_Estimate)),
TRUE ~as.matrix(env.dat$her_sp_Estimate/1000000))
ecovsig <- if(df.mods$ecovdat[m] %in% c("logmean-logsig")){
as.matrix(log(env.dat$her_sp_SE))
}else if(df.mods$ecovdat[m] %in% c("meanmil-logsigmil")){
as.matrix(log(env.dat$her_sp_SE/1000000))
}else{"est_1"}
ecov <- list(
label = "lgCopeSpring2",
mean = ecovdat,
logsigma = ecovsig,
year = env.dat$Time,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
recruitment_how = as.matrix(df.mods$ecov_how[m])
)
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Four models run successfully with an ecov on. (none have better fits than with it off)
Estimated large copepods all Jan-June (Spring) covariate (log scale, “ar1”, WHAM “est_1” sigma, m10):
knitr::include_graphics(here::here(write.dir, "m10/plots_png/results/Ecov_1_lgCopeSpring2.png"))
Estimated large copepods all Jan-June (Spring) covariate (log scale,
“rw”, WHAM “est_1” sigma, m2):
knitr::include_graphics(here::here(write.dir, "m2/plots_png/results/Ecov_1_lgCopeSpring2.png"))
Estimated large copepods all Jan-June (Spring) covariate (natural scale in millions with VAST SE in millions, “ar1”, m11):
knitr::include_graphics(here::here(write.dir, "m11/plots_png/results/Ecov_1_lgCopeSpring2.png"))
Estimated large copepods all Jan-June (Spring) covariate (natural scale in millions with VAST SE in millions, “rw”, m3):
knitr::include_graphics(here::here(write.dir, "m3/plots_png/results/Ecov_1_lgCopeSpring2.png"))
Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: catch_paa_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.37e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.31e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.87e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.19e-07
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.45e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.47e-03
## Max gradient parameter: logit_q
## TMB:sdreport() was performed successfully for this model
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.05e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 9:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 10:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.58e-06
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 11:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.33e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 12:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.75e-08
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 13:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.31e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 14:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.89e-07
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 15:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e-03
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this model
##
## Model 16:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.69e-01
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.5 -3321.7 0.8683 0.5901 -0.2314
## m2 569.8 -2755.4 0.8594 0.5907 -0.2362
## m3 571.4 -2753.8 0.8624 0.5830 -0.2315
## m4 4.7 -3320.5 0.8844 0.5846 -0.2300
## m5 570.4 -2754.8 0.8104 0.5940 -0.2448
## m6 0.0 -3325.2 0.8684 0.5901 -0.2314
## m7 566.1 -2759.1 0.8682 0.5901 -0.2314
## m8 0.6 -3324.6 0.8999 0.5844 -0.2296
## m9 566.5 -2758.7 0.9081 0.5920 -0.2379df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat
## 1 m10 2 ar1 none logmean-est_1
## 2 m14 2 ar1 controlling-lag-0-linear logmean-est_1
## 3 m2 2 rw none logmean-est_1
## 4 m6 2 rw controlling-lag-0-linear logmean-est_1
## 5 m11 2 ar1 none meanmil-logsigmil
## 6 m15 2 ar1 controlling-lag-0-linear meanmil-logsigmil
## 7 m3 2 rw none meanmil-logsigmil
## 8 m7 2 rw controlling-lag-0-linear meanmil-logsigmil
## 9 m4 2 rw none meanmil-est_1
## 10 m1 2 rw none logmean-logsig
## 11 m5 2 rw controlling-lag-0-linear logmean-logsig
## 12 m8 2 rw controlling-lag-0-linear meanmil-est_1
## 13 m9 2 ar1 none logmean-logsig
## 14 m12 2 ar1 none meanmil-est_1
## 15 m13 2 ar1 controlling-lag-0-linear logmean-logsig
## 16 m16 2 ar1 controlling-lag-0-linear meanmil-est_1
## conv pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE TRUE -1793.611 0 -3325.2 0.8684 0.5901 -0.2314
## 2 TRUE TRUE -1794.325 0.6 -3324.6 0.8999 0.5844 -0.2296
## 3 TRUE TRUE -1790.837 3.5 -3321.7 0.8683 0.5901 -0.2314
## 4 TRUE TRUE -1791.227 4.7 -3320.5 0.8844 0.5846 -0.23
## 5 TRUE TRUE -1509.570 566.1 -2759.1 0.8682 0.5901 -0.2314
## 6 TRUE TRUE -1510.349 566.5 -2758.7 0.9081 0.592 -0.2379
## 7 TRUE TRUE -1506.709 569.8 -2755.4 0.8594 0.5907 -0.2362
## 8 TRUE TRUE -1507.424 570.4 -2754.8 0.8104 0.594 -0.2448
## 9 TRUE TRUE -1506.918 571.4 -2753.8 0.8624 0.583 -0.2315
## 10 TRUE FALSE -1021.262 --- --- --- --- ---
## 11 TRUE FALSE -1021.262 --- --- --- --- ---
## 12 TRUE FALSE -1502.889 --- --- --- --- ---
## 13 TRUE FALSE -1021.262 --- --- --- --- ---
## 14 TRUE FALSE -1509.190 --- --- --- --- ---
## 15 TRUE FALSE -1021.262 --- --- --- --- ---
## 16 TRUE FALSE -1499.112 --- --- --- --- ---lgcope2 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(lgcope2, file.path(write.dir, "lgCopeSpring2compare.rds"))
lgcope2$g[[1]]
lgcope2$g[[8]]
lgcope2$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Small fall (July-December) copepods in herring fall bottom trawl survey strata
Setup
config <- "smcopeFall2"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# 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 colRun models (not here, code for reference but using this script)
## Read environmental index
env.dat <- read.csv(here::here("WHAMfits/fallsmallcopeALLindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
## Run model
for(m in 1:n.mods){
ecovdat <- dplyr::case_when(df.mods$ecovdat[m] %in% c("logmean-logsig", "logmean-est_1") ~
as.matrix(log(env.dat$her_fa_Estimate)),
TRUE ~as.matrix(env.dat$her_fa_Estimate/1000000))
ecovsig <- if(df.mods$ecovdat[m] %in% c("logmean-logsig")){
as.matrix(log(env.dat$her_fa_SE))
}else if(df.mods$ecovdat[m] %in% c("meanmil-logsigmil")){
as.matrix(log(env.dat$her_fa_SE/1000000))
}else{"est_1"}
ecov <- list(
label = "smCopeFall2",
mean = ecovdat,
logsigma = ecovsig,
year = env.dat$Time,
use_obs = use.obs, #matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
recruitment_how = as.matrix(df.mods$ecov_how[m])
)
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
One model runs successfully with an ecov on (log scale covariate, “ar1”, WHAM “est_1” sigma, m14), but does not have a better fit than with it off:
knitr::include_graphics(here::here(write.dir, "m14/plots_png/results/Ecov_1_smCopeFall2.png"))
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.53e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.20e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.03e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.31e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.34e-06
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.65e-01
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.06e-01
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 9:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 10:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.56e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 11:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.39e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 12:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.73e-08
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 13:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.31e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 14:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.48e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 15:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.34e+00
## Max gradient parameter: logit_selpars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 16:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.91e-01
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 0.1 -3332.3 0.8682 0.5901 -0.2314
## m2 734.6 -2597.8 3.2038 0.6009 -0.2487
## m3 722.9 -2609.5 0.8682 0.5901 -0.2314
## m4 0.0 -3332.4 0.8682 0.5901 -0.2314
## m5 723.0 -2609.4 0.8624 0.5830 -0.2315
## m6 0.5 -3331.9 0.8486 0.6001 -0.2333df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat
## 1 m10 2 ar1 none logmean-est_1
## 2 m2 2 rw none logmean-est_1
## 3 m14 2 ar1 controlling-lag-1-linear logmean-est_1
## 4 m4 2 rw none meanmil-est_1
## 5 m11 2 ar1 none meanmil-logsigmil
## 6 m3 2 rw none meanmil-logsigmil
## 7 m1 2 rw none logmean-logsig
## 8 m5 2 rw controlling-lag-1-linear logmean-logsig
## 9 m6 2 rw controlling-lag-1-linear logmean-est_1
## 10 m7 2 rw controlling-lag-1-linear meanmil-logsigmil
## 11 m8 2 rw controlling-lag-1-linear meanmil-est_1
## 12 m9 2 ar1 none logmean-logsig
## 13 m12 2 ar1 none meanmil-est_1
## 14 m13 2 ar1 controlling-lag-1-linear logmean-logsig
## 15 m15 2 ar1 controlling-lag-1-linear meanmil-logsigmil
## 16 m16 2 ar1 controlling-lag-1-linear meanmil-est_1
## conv pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE TRUE -1797.177 0 -3332.4 0.8682 0.5901 -0.2314
## 2 TRUE TRUE -1796.175 0.1 -3332.3 0.8682 0.5901 -0.2314
## 3 TRUE TRUE -1797.931 0.5 -3331.9 0.8486 0.6001 -0.2333
## 4 TRUE TRUE -1434.768 722.9 -2609.5 0.8682 0.5901 -0.2314
## 5 TRUE TRUE -1434.708 723 -2609.4 0.8624 0.583 -0.2315
## 6 TRUE TRUE -1427.902 734.6 -2597.8 3.2038 0.6009 -0.2487
## 7 TRUE FALSE -964.219 --- --- --- --- ---
## 8 TRUE FALSE -964.219 --- --- --- --- ---
## 9 TRUE FALSE -1794.370 --- --- --- --- ---
## 10 TRUE FALSE -1429.056 --- --- --- --- ---
## 11 TRUE FALSE -1435.415 --- --- --- --- ---
## 12 TRUE FALSE -964.219 --- --- --- --- ---
## 13 TRUE FALSE -1433.763 --- --- --- --- ---
## 14 TRUE FALSE -964.219 --- --- --- --- ---
## 15 TRUE FALSE -1426.200 --- --- --- --- ---
## 16 TRUE FALSE -1355.745 --- --- --- --- ---smcope2 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcope2, file.path(write.dir, "smCopeFall2compare.rds"))
smcope2$g[[1]]
smcope2$g[[8]]
smcope2$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Small Sept-February copepods in herring larval area
Setup
config <- "smcopeSepFeb2"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# 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 colRun models (not here, code for reference but using this script)
## Read environmental index
env.dat <- read.csv(here::here("WHAMfits/sepfebsmallcopeALLlarvareaindex.csv"), header=T)
# 2020 is missing
env.dat[env.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(env.dat)[1])
use.obs[39,] <- 0
## Run model
for(m in 1:n.mods){
ecovdat <- dplyr::case_when(df.mods$ecovdat[m] %in% c("logmean-logsig", "logmean-est_1") ~
as.matrix(log(env.dat$her_larv_Estimate)),
TRUE ~as.matrix(env.dat$her_larv_Estimate/1000000))
ecovsig <- if(df.mods$ecovdat[m] %in% c("logmean-logsig")){
as.matrix(log(env.dat$her_larv_SE))
}else if(df.mods$ecovdat[m] %in% c("meanmil-logsigmil")){
as.matrix(log(env.dat$her_larv_SE/1000000))
}else{"est_1"}
ecov <- list(
label = "smCopeSepFeb2",
mean = ecovdat,
logsigma = ecovsig,
year = env.dat$Time,
use_obs = use.obs, #matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
recruitment_how = as.matrix(df.mods$ecov_how[m])
)
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
The model with ecov on (log scale covariate, “ar1”, WHAM “est_1” sigma, m14), fits best, better than with ecov off for both ar1 and rw:
knitr::include_graphics(here::here(write.dir, "m14/plots_png/results/Ecov_1_smCopeSepFeb2.png"))
mod.list <- paste0(write.dir, df.mods$Model,".rds")
mod.list <- mod.list[-16] # last model crashed no output saved
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.24e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.88e-11
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.88e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.31e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.25e-05
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.51e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.55e+00
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 9:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.28e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 10:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.10e-11
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 11:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.46e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 12:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.03e-09
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 13:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.31e-13
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 14:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.40e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 15:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.37e+00
## Max gradient parameter: F_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertibleopt_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 <- df.mods[-16,]
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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 0.5 -3319.3 0.8683 0.5901 -0.2314
## m2 692.5 -2627.3 0.8541 0.5879 -0.2311
## m3 693.7 -2626.1 0.8682 0.5901 -0.2314
## m4 0.9 -3318.9 0.8683 0.5901 -0.2314
## m5 688.6 -2631.2 0.8793 0.6035 -0.2472
## m6 0.0 -3319.8 0.8140 0.5966 -0.2327df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat
## 1 m14 2 ar1 controlling-lag-1-linear logmean-est_1
## 2 m2 2 rw none logmean-est_1
## 3 m10 2 ar1 none logmean-est_1
## 4 m11 2 ar1 none meanmil-logsigmil
## 5 m3 2 rw none meanmil-logsigmil
## 6 m4 2 rw none meanmil-est_1
## 7 m1 2 rw none logmean-logsig
## 8 m5 2 rw controlling-lag-1-linear logmean-logsig
## 9 m6 2 rw controlling-lag-1-linear logmean-est_1
## 10 m7 2 rw controlling-lag-1-linear meanmil-logsigmil
## 11 m8 2 rw controlling-lag-1-linear meanmil-est_1
## 12 m9 2 ar1 none logmean-logsig
## 13 m12 2 ar1 none meanmil-est_1
## 14 m13 2 ar1 controlling-lag-1-linear logmean-logsig
## 15 m15 2 ar1 controlling-lag-1-linear meanmil-logsigmil
## conv pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE TRUE -1791.900 0 -3319.8 0.814 0.5966 -0.2327
## 2 TRUE TRUE -1789.657 0.5 -3319.3 0.8683 0.5901 -0.2314
## 3 TRUE TRUE -1790.469 0.9 -3318.9 0.8683 0.5901 -0.2314
## 4 TRUE TRUE -1445.582 688.6 -2631.2 0.8793 0.6035 -0.2472
## 5 TRUE TRUE -1442.652 692.5 -2627.3 0.8541 0.5879 -0.2311
## 6 TRUE TRUE -1443.075 693.7 -2626.1 0.8682 0.5901 -0.2314
## 7 TRUE FALSE -965.948 --- --- --- --- ---
## 8 TRUE FALSE -965.948 --- --- --- --- ---
## 9 TRUE FALSE -1785.746 --- --- --- --- ---
## 10 TRUE FALSE -1444.922 --- --- --- --- ---
## 11 TRUE FALSE -1445.372 --- --- --- --- ---
## 12 TRUE FALSE -965.948 --- --- --- --- ---
## 13 TRUE FALSE -1442.064 --- --- --- --- ---
## 14 TRUE FALSE -965.948 --- --- --- --- ---
## 15 TRUE FALSE -1447.832 --- --- --- --- ---smcopelarv2 <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(smcopelarv2, file.path(write.dir, "smCopeSepFeb2compare.rds"))
smcopelarv2$g[[1]]
smcopelarv2$g[[8]]
smcopelarv2$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Only compare logmean-est_1 fits that worked
mod.list <- paste0(write.dir, c("m14", "m2", "m10"),".rds")
mods <- lapply(mod.list, readRDS)
names(mods) <- c("m14", "m2", "m10")
smcopelarv2_best4 <- compare_wham_models(mods, fdir = write.dir)## dAIC AIC rho_R rho_SSB rho_Fbar
## m14 0.0 -3319.8 0.8140 0.5966 -0.2327
## m2 0.5 -3319.3 0.8683 0.5901 -0.2314
## m10 0.9 -3318.9 0.8683 0.5901 -0.2314saveRDS(smcopelarv2_best4, file.path(write.dir, "smCopeSepFeb2compare_best3.rds"))
smcopelarv2_best4$g[[1]]
smcopelarv2_best4$g[[8]]
smcopelarv2_best4$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Small copepods Sep-Feb in the fall herring survey strata also worked, maybe better
Two models ran successfully with an ecov on (log scale covariate, “rw”, WHAM “est_1” sigma, m6), this fits best, better than with ecov off:
write.dir <- here::here("WHAMfits/mm192_smcopeSepFeb2_fallstrata")
knitr::include_graphics(here::here(write.dir, "m6/plots_png/results/Ecov_1_smCopeSepFeb2.png"))
and ecov on (log scale covariate, “ar1”, WHAM “est_1” sigma, m14), this fits second best, better than with ecov off:
knitr::include_graphics(here::here(write.dir, "m14/plots_png/results/Ecov_1_smCopeSepFeb2.png"))
Only compare logmean-est_1 fits that worked
mod.list <- paste0(write.dir, "/", c("m6", "m14", "m2", "m10"),".rds")
mods <- lapply(mod.list, readRDS)
names(mods) <- c("m6", "m14", "m2", "m10")
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.80e-11
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.83e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.78e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.49e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this modelopt_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)
NLL <- sapply(mods, function(x) round(x$opt$objective,3))
smcopelarv2_best4 <- compare_wham_models(mods, fdir = write.dir)## dAIC AIC rho_R rho_SSB rho_Fbar
## m6 0.0 -3322.0 0.7592 0.6023 -0.2327
## m14 1.1 -3320.9 0.7943 0.5971 -0.2329
## m2 1.6 -3320.4 0.8682 0.5901 -0.2314
## m10 2.5 -3319.5 0.8683 0.5901 -0.2314saveRDS(smcopelarv2_best4, file.path(write.dir, "smCopeSepFeb2compare_best4.rds"))
smcopelarv2_best4$g[[1]]
smcopelarv2_best4$g[[8]]
smcopelarv2_best4$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Try two covariates: lgCopeSpring lag 0 + smCopeSepFeb lag 1
Possibly insane. May not work? But these are the top two drivers in Adelle’s boosted regression tree.
The reason to try more than one is that probably just one factor doesn’t “control” recruitment.
Setup
config <- "lgcopeSpring_smcopeSepFeb"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# 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
# make smaller
df.mods <- data.frame(Recruitment = c(rep(2, 4)),
ecov_process = c(rep(c("rw","ar"),2)),
ecov_how = c(rep("none" ,2),
rep("controlling-lag-0-linear",2)),
ecovdat = rep("logmean-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
## Read environmental indices
lgcope.dat <- read.csv(here::here("WHAMfits/springlargecopeindex.csv"), header=T)
smcope.dat <- read.csv(here::here("WHAMfits/sepfebsmallcopeALLlarvareaindex.csv"), header=T)
# 2020 is missing for fall small copepods
smcope.dat[smcope.dat == 0] <- NA
# don't use the NA value
use.obs <- matrix(1, ncol=1, nrow=dim(smcope.dat)[1])
use.obs[39,] <- 0Run the models, the loop didnt work, only doing the 4 above, run from this script
## Try just running a model with both turned off and both turned on, log scale and est_1
## having issues getting the correct class of objects in this loop
#1: rw none logmean-est1 both
#2: ar1 none logmean-est1 both
#3: rw controlling-lag-0-linear and controlling-lag-1-linear logmean-est1 both
#4: ar1 controlling-lag-0-linear and controlling-lag-1-linear logmean-est1 both
ecov1 <- list(
label = c("lgCopeSpring2", "smCopeSepFeb2"),
mean = cbind(as.matrix(log(lgcope.dat$her_sp_Estimate)),
as.matrix(log(smcope.dat$her_larv_Estimate))),
logsigma = c("est_1", "est_1"),
year = lgcope.dat$Time,
use_obs = cbind(matrix(1, ncol=1, nrow=dim(lgcope.dat)[1]), use.obs), # use all obs (all = 1)
process_model = c("rw","rw"), # "rw" or "ar1"
recruitment_how = rbind(as.matrix("none"), as.matrix("none"))
)
ecov2 <- list(
label = c("lgCopeSpring2", "smCopeSepFeb2"),
mean = cbind(as.matrix(log(lgcope.dat$her_sp_Estimate)),
as.matrix(log(smcope.dat$her_larv_Estimate))),
logsigma = c("est_1", "est_1"),
year = lgcope.dat$Time,
use_obs = cbind(matrix(1, ncol=1, nrow=dim(lgcope.dat)[1]), use.obs), # use all obs (all = 1)
process_model = c("ar1","ar1"), # "rw" or "ar1"
recruitment_how = rbind(as.matrix("none"), as.matrix("none"))
)
ecov3 <- list(
label = c("lgCopeSpring2", "smCopeSepFeb2"),
mean = cbind(as.matrix(log(lgcope.dat$her_sp_Estimate)),
as.matrix(log(smcope.dat$her_larv_Estimate))),
logsigma = c("est_1", "est_1"),
year = lgcope.dat$Time,
use_obs = cbind(matrix(1, ncol=1, nrow=dim(lgcope.dat)[1]), use.obs), # use all obs (all = 1)
process_model = c("rw","rw"), # "rw" or "ar1"
recruitment_how = rbind(as.matrix("controlling-lag-0-linear"), as.matrix("controlling-lag-1-linear"))
)
ecov4 <- list(
label = c("lgCopeSpring2", "smCopeSepFeb2"),
mean = cbind(as.matrix(log(lgcope.dat$her_sp_Estimate)),
as.matrix(log(smcope.dat$her_larv_Estimate))),
logsigma = c("est_1", "est_1"),
year = lgcope.dat$Time,
use_obs = cbind(matrix(1, ncol=1, nrow=dim(lgcope.dat)[1]), use.obs), # use all obs (all = 1)
process_model = c("ar1","ar1"), # "rw" or "ar1"
recruitment_how = rbind(as.matrix("controlling-lag-0-linear"), as.matrix("controlling-lag-1-linear"))
)
# run mod 1
ecovinput <- set_ecov(input, ecov=ecov1)
mod <- fit_wham(ecovinput, do.osa = T)
saveRDS(mod, file.path(write.dir, "m1.rds"))
plot_wham_output(mod=mod, dir.main=file.path(write.dir,"m1"), out.type='html')
# run mod 2
ecovinput <- set_ecov(input, ecov=ecov2)
mod <- fit_wham(ecovinput, do.osa = T)
saveRDS(mod, file.path(write.dir, "m2.rds"))
plot_wham_output(mod=mod, dir.main=file.path(write.dir,"m2"), out.type='html')
# run mod 3
ecovinput <- set_ecov(input, ecov=ecov3)
mod <- fit_wham(ecovinput, do.osa = T)
saveRDS(mod, file.path(write.dir, "m3.rds"))
plot_wham_output(mod=mod, dir.main=file.path(write.dir,"m3"), out.type='html')
# run mod 4
ecovinput <- set_ecov(input, ecov=ecov4)
mod <- fit_wham(ecovinput, do.osa = T)
saveRDS(mod, file.path(write.dir, "m4.rds"))
plot_wham_output(mod=mod, dir.main=file.path(write.dir,"m4"), out.type='html')Compare mods
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.57e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.20e-06
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.58e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 7.05e-07
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed successfully for this modelopt_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 <- df.mods[-16,]
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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 3.2 -3269.5 0.8683 0.5901 -0.2314
## m2 0.0 -3272.7 0.8682 0.5901 -0.2314
## m3 3.4 -3269.3 0.8001 0.5976 -0.2321
## m4 0.0 -3272.7 0.8529 0.5922 -0.2313df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat conv
## 1 m2 2 ar none logmean-est_1 TRUE
## 2 m4 2 ar controlling-lag-0-linear logmean-est_1 TRUE
## 3 m1 2 rw none logmean-est_1 TRUE
## 4 m3 2 rw controlling-lag-0-linear logmean-est_1 TRUE
## pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE -1771.328 0.0 -3272.7 0.8682 0.5901 -0.2314
## 2 TRUE -1773.347 0.0 -3272.7 0.8529 0.5922 -0.2313
## 3 TRUE -1767.742 3.2 -3269.5 0.8683 0.5901 -0.2314
## 4 TRUE -1769.627 3.4 -3269.3 0.8001 0.5976 -0.2321Compare
lgandsmall <- compare_wham_models(mods, fdir = write.dir)## dAIC AIC rho_R rho_SSB rho_Fbar
## m2 0.0 -3272.7 0.8682 0.5901 -0.2314
## m4 0.0 -3272.7 0.8529 0.5922 -0.2313
## m1 3.2 -3269.5 0.8683 0.5901 -0.2314
## m3 3.4 -3269.3 0.8001 0.5976 -0.2321saveRDS(lgandsmall, file.path(write.dir, "smCopeSepFeb2compare_best3.rds"))
lgandsmall$g[[1]]
lgandsmall$g[[8]]
lgandsmall$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
The ar1 models have identical fits with and without the covariates.
Interesting. Maybe didnt work?
Temperature covariate test
Same format as for zooplankton indicators
Setup
config <- "LarvalTempDuration"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192/mm192_meanrecpar.rds"))
input <- mm192mod$input
# 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 colRun models (not here, code for reference but using this script)
## Read environmental index
env.dat <- read.csv(here::here("WHAMfits/Duration.Optimal.SST.Sept-Dec.csv"), header=T)
## Run model
for(m in 1:n.mods){
ecovdat <- dplyr::case_when(df.mods$ecovdat[m] %in% c("logmean-est_1") ~
as.matrix(log(env.dat$duration)),
TRUE ~as.matrix(env.dat$duration))
ecov <- list(
label = "LarvalTempDuration",
mean = ecovdat,
logsigma = "est_1",
year = env.dat$year,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
recruitment_how = as.matrix(df.mods$ecov_how[m])
)
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Compare models (from https://timjmiller.github.io/wham/articles/ex2_CPI_recruitment.html)
The model with ecov on (log scale covariate, “rw”, WHAM “est_1” sigma, m4), fits best, better than with ecov off for both ar1 and rw:
knitr::include_graphics(here::here(write.dir, "m4/plots_png/results/Ecov_1_LarvalTempDuration.png"))
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.78e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 5.31e-09
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.75e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 4.10e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.33e-08
## Max gradient parameter: Ecov_process_pars
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 9.75e-11
## Max gradient parameter: Ecov_beta_R
## TMB:sdreport() was performed successfully for this model
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 8.46e-12
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this model
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.63e-11
## Max gradient parameter: log_NAA_sigma
## TMB:sdreport() was performed successfully for this modelopt_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 <- df.mods[-16,]
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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 354.6 -3038.5 0.8594 0.5907 -0.2362
## m2 354.2 -3038.9 0.6473 0.5976 -0.2321
## m3 0.4 -3392.7 0.8682 0.5901 -0.2314
## m4 0.0 -3393.1 0.6543 0.5984 -0.2323
## m5 357.9 -3035.2 0.6643 0.5971 -0.2322
## m6 4.1 -3389.0 0.8682 0.5901 -0.2314
## m7 3.8 -3389.3 0.6714 0.5980 -0.2324df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat conv
## 1 m4 2 rw controlling-lag-1-linear logmean-est_1 TRUE
## 2 m3 2 rw none logmean-est_1 TRUE
## 3 m8 2 ar1 controlling-lag-1-linear logmean-est_1 TRUE
## 4 m7 2 ar1 none logmean-est_1 TRUE
## 5 m2 2 rw controlling-lag-1-linear mean-est_1 TRUE
## 6 m1 2 rw none mean-est_1 TRUE
## 7 m6 2 ar1 controlling-lag-1-linear mean-est_1 TRUE
## 8 m5 2 ar1 none mean-est_1 TRUE
## pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE -1827.543 0 -3393.1 0.6543 0.5984 -0.2323
## 2 TRUE -1826.370 0.4 -3392.7 0.8682 0.5901 -0.2314
## 3 TRUE -1826.627 3.8 -3389.3 0.6714 0.598 -0.2324
## 4 TRUE -1825.511 4.1 -3389 0.8682 0.5901 -0.2314
## 5 TRUE -1650.465 354.2 -3038.9 0.6473 0.5976 -0.2321
## 6 TRUE -1649.265 354.6 -3038.5 0.8594 0.5907 -0.2362
## 7 TRUE -1649.607 357.9 -3035.2 0.6643 0.5971 -0.2322
## 8 FALSE -1644.890 --- --- --- --- ---With NAA RE on, relatively weak impact of covariate
knitr::include_graphics(here::here(write.dir, "/m4/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png"))
larvtemp <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(larvtemp, file.path(write.dir, "larvtempduration.rds"))
larvtemp$g[[1]]
larvtemp$g[[8]]
larvtemp$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Larval temperature duration when NAA random effects are OFF
Setup
## Setup
config <- "LarvalTempDuration_NAA_REoff"
# name the model directory
name <- paste0("mm192_", config)
write.dir <- here::here(sprintf("WHAMfits/%s/", name))
if(!dir.exists(write.dir)) {
dir.create(write.dir)
}
mm192mod <- readRDS(here::here("WHAMfits/mm192_nonaa/mm192_nonaa.rds"))
input <- mm192mod$input
# 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 colRun mods (from this script)
## Read environmental index
env.dat <- read.csv(here::here("WHAMfits/Duration.Optimal.SST.Sept-Dec.csv"), header=T)
## Run model
for(m in 1:n.mods){
ecovdat <- dplyr::case_when(df.mods$ecovdat[m] %in% c("logmean-est_1") ~
as.matrix(log(env.dat$duration)),
TRUE ~as.matrix(env.dat$duration))
ecov <- list(
label = "LarvalTempDuration",
mean = ecovdat,
logsigma = "est_1",
year = env.dat$year,
use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (all = 1)
process_model = df.mods$ecov_process[m], # "rw" or "ar1"
recruitment_how = as.matrix(df.mods$ecov_how[m])
)
ecovinput <- set_ecov(input, ecov=ecov)
mod <- fit_wham(ecovinput, do.osa = T)
# Save model
saveRDS(mod, file.path(write.dir, paste0(df.mods$Model[m],".rds")))
# Plot output in new subfolder
plot_wham_output(mod=mod, dir.main=file.path(write.dir,df.mods$Model[m]), out.type='html')
}Compare models
The model with ecov on (log scale covariate, “rw”, WHAM “est_1” sigma, m4), fits best, better than with ecov off for both ar1 and rw:
knitr::include_graphics(here::here(write.dir, "/m4/plots_png/results/Ecov_1_LarvalTempDuration.png"))
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')## Model 1:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 3.30e-12
## Max gradient parameter: logit_q
## TMB:sdreport() was performed successfully for this model
##
## Model 2:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.56e-08
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 3:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.46e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 4:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.52e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 5:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.42e-12
## Max gradient parameter: logit_q
## TMB:sdreport() was performed successfully for this model
##
## Model 6:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.35e-03
## Max gradient parameter: logit_q
## TMB:sdreport() was performed for this model, but it appears hessian was not invertible
##
## Model 7:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 2.15e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this model
##
## Model 8:
## stats:nlminb thinks the model has converged: mod$opt$convergence == 0
## Maximum gradient component: 1.76e-12
## Max gradient parameter: F_pars
## TMB:sdreport() was performed successfully for this modelopt_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 <- df.mods[-16,]
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)## dAIC AIC rho_R rho_SSB rho_Fbar
## m1 380.6 -2886.0 3.5015 0.3662 -0.1421
## m2 354.1 -2912.5 1.5766 0.3325 -0.1206
## m3 26.4 -3240.2 3.5003 0.3664 -0.1424
## m4 0.0 -3266.6 1.6107 0.3339 -0.1210
## m5 384.3 -2882.3 3.5004 0.3673 -0.1437
## m6 30.1 -3236.5 3.5003 0.3664 -0.1424
## m7 4.2 -3262.4 1.6872 0.3341 -0.1213df.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
df.mods## Model Recruitment ecov_process ecov_how ecovdat conv
## 1 m4 2 rw controlling-lag-1-linear logmean-est_1 TRUE
## 2 m8 2 ar1 controlling-lag-1-linear logmean-est_1 TRUE
## 3 m3 2 rw none logmean-est_1 TRUE
## 4 m7 2 ar1 none logmean-est_1 TRUE
## 5 m2 2 rw controlling-lag-1-linear mean-est_1 TRUE
## 6 m1 2 rw none mean-est_1 TRUE
## 7 m5 2 ar1 none mean-est_1 TRUE
## 8 m6 2 ar1 controlling-lag-1-linear mean-est_1 TRUE
## pdHess NLL dAIC AIC rho_R rho_SSB rho_Fbar
## 1 TRUE -1762.281 0 -3266.6 1.6107 0.3339 -0.121
## 2 TRUE -1761.214 4.2 -3262.4 1.6872 0.3341 -0.1213
## 3 TRUE -1748.086 26.4 -3240.2 3.5003 0.3664 -0.1424
## 4 TRUE -1747.227 30.1 -3236.5 3.5003 0.3664 -0.1424
## 5 TRUE -1585.263 354.1 -2912.5 1.5766 0.3325 -0.1206
## 6 TRUE -1570.980 380.6 -2886 3.5015 0.3662 -0.1421
## 7 TRUE -1570.163 384.3 -2882.3 3.5004 0.3673 -0.1437
## 8 FALSE -1584.188 --- --- --- --- ---Without NAA RE on, bigger impact of covariate
knitr::include_graphics(here::here(write.dir, "/m4/plots_png/diagnostics/NAA_4panel_stock_1_region_1_age_1.png"))
larvtemp <- compare_wham_models(mods2, do.table= FALSE, fdir = write.dir)
saveRDS(larvtemp, file.path(write.dir, "larvtempduration.rds"))
larvtemp$g[[1]]
larvtemp$g[[8]]
larvtemp$g[[9]]
knitr::include_graphics(here::here(write.dir, "compare_png/compare_rel_status_kobe.png"))
Cumulative effects
Try the haddock predation plus the larval optimal temperature duration
Zooplankton indices are fitting, but they are fitting with the opposite direction we would expect; slope parameters are negative, suggesting an inverse relationship with herring recruitment. We expected food to have a positive effect on herring recruitment.
However, the influence estimated for haddock predation and larval optimal temperature were aligned with our hyotheses.
So, setup the model
Run models
using the script, copied here for referene
Compare models