tests/testthat/test_ex2_CPI_SNEMAYT.R

In timjmiller/wham: Woods Hole Assessment Model (WHAM)

# Test WHAM example 2: Cold Pool Index effect on SNEMA Yellowtail Flounder recruitment
# Replicate Miller et al 2016 results
#   adds environmental covariate (CPI, treated as rw)
#   uses 5 indices (ex 1: only 2 indices)
#   only fit to 1973-2011 data (ex 1: 1973-2016)
#   age compositions = 5, logistic normal pool zero obs (ex 1: 7, logistic normal missing zero obs)
#   selectivity = logistic (ex 1: age-specific)

# devtools::load_all()
# btime <- Sys.time(); devtools::test("/home/bstock/Documents/wham"); etime <- Sys.time(); runtime = etime - btime;
# btime <- Sys.time(); testthat::test_file("/home/bstock/Documents/wham/tests/testthat/test_ex2_CPI_SNEMAYT.R"); etime <- Sys.time(); runtime = etime - btime;
# 12 min

context("Ex 2: CPI on yellowtail recruitment")

test_that("Ex 2 works",{
# get results to check NLL and par estimates
path_to_examples <- system.file("extdata", package="wham")
ex2_test_results <- readRDS(file.path(path_to_examples,"ex2_test_results.rds"))

asap3 <- read_asap3_dat(file.path(path_to_examples,"ex2_SNEMAYT.dat"))
env.dat <- read.csv(file.path(path_to_examples,"CPI.csv"), header=T)

df.mods <- data.frame(Recruitment = c(2,2,3,3,3,3,4),
                      Ecov_process = c(rep("rw",4),rep("ar1",3)),
                      Ecov_how = c(0,1,0,2,2,1,1), 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

tmp.dir <- tempdir(check=TRUE)
mods <- list()
for(m in 1:n.mods){
  ecov <- list(
    label = "CPI",
    mean = as.matrix(env.dat$CPI),
    logsigma = as.matrix(log(env.dat$CPI_sigma)),
    year = env.dat$Year,
    use_obs = matrix(1, ncol=1, nrow=dim(env.dat)[1]), # use all obs (=1)
    lag = 1, # CPI in year t affects recruitment in year t+1
    process_model = df.mods$Ecov_process[m], # "rw" or "ar1"
    where = c("none","recruit")[as.logical(df.mods$Ecov_how[m])+1], # CPI affects recruitment
    how = df.mods$Ecov_how[m]) # 0 = no effect (but still fit Ecov to compare AIC), 1 = controlling (dens-indep mortality), 2 = limiting (carrying capacity), 3 = lethal (threshold), 4 = masking (metabolism/growth), 5 = directive (behavior)

  # (not used in this vignette) can set Ecov = NULL to fit model without Ecov data
  if(is.na(df.mods$Ecov_process[m])) Ecov = NULL

  # Generate wham input from ASAP3 and Ecov data
  input <- prepare_wham_input(asap3, recruit_model = df.mods$Recruitment[m],
                              model_name = "Ex 2: SNEMA Yellowtail Flounder with CPI effects on R",
                              ecov = ecov,
                              NAA_re = list(sigma="rec+1", cor="iid"),
                              age_comp = "logistic-normal-pool0") # logistic normal pool 0 obs

  # Selectivity = logistic, not age-specific as in ex1
  #   2 pars per block instead of n.ages
  #   sel pars of indices 4/5 fixed at 1.5, 0.1 (specified via neg phase in ex2_SNEMAYT.dat)
  input$par$logit_selpars[1:4,7:8] <- 0 # last 2 rows will not be estimated (mapped to NA)

  # Fit model
  mods[[m]] <- suppressWarnings(fit_wham(input, do.retro=T, do.osa=F, MakeADFun.silent = TRUE, retro.silent = TRUE))
  suppressWarnings(plot_wham_output(mod=mods[[m]], dir.main=tmp.dir))
}
# mod.list <- paste0("/home/bstock/Documents/wham/sandbox/ex2/",grep(".rds",list.files("/home/bstock/Documents/wham/sandbox/ex2"),value=TRUE))
# mods <- lapply(mod.list, readRDS)
# ex2_test_results <- list()
# ex2_test_results$pars <- lapply(mods, function(x) as.numeric(x$opt$par))
# ex2_test_results$nll <- sapply(mods, function(x) x$opt$obj)
# saveRDS(ex2_test_results, file="/home/bstock/Documents/wham/inst/extdata/ex2_test_results.rds")

# hard to see which model fails bc they're indexed by m
# print out each one by one
mcheck <- check_convergence(mods[[1]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[1]]$opt$par), ex2_test_results$pars[[1]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[1]]$opt$obj), ex2_test_results$nll[1], tolerance=1e-6) # nll

mcheck <- check_convergence(mods[[2]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[2]]$opt$par), ex2_test_results$pars[[2]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[2]]$opt$obj), ex2_test_results$nll[2], tolerance=1e-6) # nll

mcheck <- check_convergence(mods[[3]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[3]]$opt$par), ex2_test_results$pars[[3]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[3]]$opt$obj), ex2_test_results$nll[3], tolerance=1e-6) # nll

mcheck <- check_convergence(mods[[4]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[4]]$opt$par), ex2_test_results$pars[[4]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[4]]$opt$obj), ex2_test_results$nll[4], tolerance=1e-6) # nll

mcheck <- check_convergence(mods[[5]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[5]]$opt$par), ex2_test_results$pars[[5]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[5]]$opt$obj), ex2_test_results$nll[5], tolerance=1e-6) # nll

# m6 does not converge now
# mcheck <- check_convergence(mods[[6]], ret=TRUE)
# expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
# expect_false(mcheck$na_sdrep) # sdrep should succeed
# expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
# expect_equal(as.numeric(mods[[6]]$opt$par), ex2_test_results$pars[[6]], tolerance=1e-3) # parameter values
# expect_equal(as.numeric(mods[[6]]$opt$obj), ex2_test_results$nll[6], tolerance=1e-6) # nll

mcheck <- check_convergence(mods[[7]], ret=TRUE)
expect_equal(mcheck$convergence, 0) # opt$convergence should be 0
expect_false(mcheck$na_sdrep) # sdrep should succeed
expect_lt(mcheck$maxgr, 1e-5) # maximum gradient should be < 1e-06
expect_equal(as.numeric(mods[[7]]$opt$par), ex2_test_results$pars[[7]], tolerance=1e-3) # parameter values
expect_equal(as.numeric(mods[[7]]$opt$obj), ex2_test_results$nll[7], tolerance=1e-6) # nll
})

# # remove files created during testing
# teardown(unlink(tmp.dir, recursive=TRUE))