How does sdmTMB calculate the variance within get_index()?

[seananderson]

Question by email:

I was wondering if you could shed some light on how the get_index() function estimates variance for the summed predictions?
[...] my (very limited) understanding is that the predictions are summed and then the Delta method is used on the backend to estimate the variance of the summed predictions. Is that sort of on the right track?

[seananderson]

The standard error on the index is calculated using the generalized delta method as implemented in TMB. Yes, you’re on the right track, but the ‘generalized’ version is used because of the random effects. That is described in ?TMB::sdreport and probably in the TMB paper. There’s also some bias on the mean introduced because of the exp() transformation before summing up the biomass/abundance. That can be corrected in a variety of ways. If bias_correct = TRUE then this method is implemented. This can be quite slow and memory intensive, but less so with recent versions of sdmTMB where it only bias corrects the necessary values and uses an approach Kasper recently implemented in TMB that gets the same result.

Another option we’ve used sometimes is with the get_index_sims() function. That needs output from predict() with the nsim argument > 0… perhaps something like 500 or 1000. There’s an example at the end of the index vignette, I believe. This takes draws from the joint precision matrix for all the parameters in the internal transformed space (assumes an MVN parameter distribution), and calculates the index with all those draws, and then summarizes that. This usually agrees with the bias_correct = TRUE option, but it’s possible it will fail in some edge cases, so the bias_correct = TRUE would be better to use for more formal applications/advice. See also ?gather_sims() and ?spread_sims().

One last option is to use the MCMC samples from tmbstan and pass those to predict(..., tmbstan_model = ...) and then get_index_sims(). There's an example in ?extract_mcmc.