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

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
wham.dir <- find.package("wham")
knitr::opts_knit$set(root.dir = file.path(wham.dir,"extdata"))
library(knitr)
library(kableExtra)

In this vignette we show how to use the do.fit = FALSE flag to debug a WHAM model.

# devtools::install_github("timjmiller/wham", dependencies=TRUE)
library(wham)

Load the ASAP3 data file.

wham.dir <- find.package("wham")
asap3 <- read_asap3_dat(file.path(wham.dir,"extdata","ex7_SNEMAYT.dat"))

Prepare the WHAM model (m3 from example 1).

input <- prepare_wham_input(asap3, recruit_model=2, model_name="Ex 1: SNEMA Yellowtail Flounder",
                              NAA_re = list(sigma="rec+1", cor="iid"))

Try to fit the model... uh oh.

mod <- fit_wham(input, do.osa = F, do.retro = F)

What's wrong? It looks like the likelihood function is returning NaN. It is often easier to diagnose problems like this using the unoptimized model, i.e. look at everything using the initial parameter values. WHAM includes a do.fit = F flag in fit_wham to return the unoptimized model object returned by TMB::MakeADFun. Let's see how it works.

mod <- fit_wham(input, do.fit = F)

This runs without an error. The optimization failed because the likelihood was NaN, and now we can see which of the likelihood components was responsible. To do so, we need to look at the REPORTed objects with "nll" in their name. We can use the $report() function from TMB.

therep = mod$report()

See all of the REPORTed objects.

names(therep)

Now just get the objects with "nll" in their name, and sum over all individual values.

sapply(grep("nll",names(therep),value=T), function(x) sum(therep[[x]]))

The likelihood components that are equal to 0 are not used in the model (no random effects on M, Ecov, selectivity, etc.). nll is the total likelihood and is NaN. The one troublesome component is nll_agg_catch. This is the total (aggregate) catch from the fishery in each year. It could be an issue with the catch data or the model-predicted catch, since both are in the likelihood calculation.

Search the WHAM .cpp file for "nll_agg_catch". We find that on line 888, this depends on agg_catch (the catch data) and pred_catch (model-predicted catch).

The catch data looks ok.

input$data$agg_catch

The predicted catch is calculated on lines 879-880. pred_catch depends on NAA (numbers-at-age), FAA (F at age), ZAA (Z at age), and waa (weight-at-age data).

NAA looks ok.

therep$NAA

FAA looks ok - no issues with F or selectivity.

therep$FAA[,1,] # middle dim is n.fleets = 1

ZAA looks ok - no issues with M either.

therep$ZAA

Ah, here is the problem. The weight at age data has an entry of -99. This means that pred_catch is negative on line 880, and we take the log of a negative number on line 883.