Calculate a population index via simulation from the joint precision matrix

Calculate a population index via simulation from the joint precision matrix. Compared to get_index(), this version can be dramatically faster if bias correction was turned on in get_index() while being approximately equivalent. This is an experimental function. We have yet to find a model where this function fails to provide a reasonable result, but make no guarantees.

get_index_sims(
  obj,
  level = 0.95,
  return_sims = FALSE,
  area = rep(1, nrow(obj)),
  est_function = stats::median,
  agg_function = function(x) sum(exp(x))
)

Arguments

obj

predict.sdmTMB() output with sims > 0.

level

The confidence level.

return_sims

Logical. Return simulation draws? The default (FALSE) is a quantile summary of those simulation draws.

area

A vector of grid cell/polyon areas for each year-grid cell (row of data) in obj. Adjust this if cells are not of unit area or not all the same area (e.g., some cells are partially over land/water). Note that the area vector is added as log(area) to the raw values in obj. In other words, the function assumes a log link, which typically makes sense.

est_function

Function to summarize the estimate (the expected value). mean() would be an alternative to median().

agg_function

Function to aggregate samples within each time slice. Assuming a log link, the sum(exp(x) * area) default makes sense.

Details

Can also be used to produce an index from a model fit with tmbstan.

This function does nothing more than summarize and reshape the matrix of simulation draws into a data frame.

See also

get_index()

Examples

if (inla_installed()) {

m <- sdmTMB(density ~ 0 + as.factor(year) + depth_scaled + depth_scaled2,
  data = pcod_2011, mesh = pcod_mesh_2011, family = tweedie(link = "log"),
  time = "year"
)
qcs_grid_2011 <- subset(qcs_grid, year >= 2011)
p <- predict(m, newdata = qcs_grid_2011, sims = 100)
x <- get_index_sims(p)
x_sims <- get_index_sims(p, return_sims = TRUE)

if (require("ggplot2", quietly = TRUE)) {
  ggplot(x, aes(year, est, ymin = lwr, ymax = upr)) +
    geom_line() +
    geom_ribbon(alpha = 0.4)
  ggplot(x_sims, aes(as.factor(year), .value)) +
    geom_violin()
}

}