Example with built-in delta model

Here, we will show an example of fitting using the built-in delta functionality, as well as how to build each model component separately and then combine. The built-in approach is convenient, allows for parameters to be shared across components, and allows for calculation of derived quantities such as standardized indexes (get_index()) with internally calculated standard errors.

We will use a dataset built into the sdmTMB package: trawl survey data for Pacific Cod in Queen Charlotte Sound, British Columbia, Canada. The density units are kg/km². Here, X and Y are coordinates in UTM zone 9.

We will first create a mesh that we will use for all the models.

pcod_mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)

Then we can fit a model of Pacific cod density using a delta-gamma model, including a smoothed effect of depth.

fit_dg <- sdmTMB(density ~ 1 + s(depth),
  data = pcod,
  mesh = pcod_mesh,
  time = "year",
  family = delta_gamma()
)

The default in built-in delta models is for the formula, spatial and spatio-temporal structure, and anisotropy to be shared between the two model components. However, some elements (formula, spatial, spatiotemporal, and share_range) can also be specified independently for each model using a list format within the function argument (see examples below). The first element of the list is for the binomial component and the second element is for the positive component (e.g., Gamma). Some elements must be shared for now (e.g., smoothers, spatially varying coefficients, and time-varying coefficients).

To specify the settings for spatial and spatio-temporal effects in each model component, create a list of settings within the spatial and spatiotemporal arguments. For example, spatial = list("on", "off"), spatiotemporal = list("off", "rw").

We could similarly specify a different formula for each component of the model, using list(y ~ x1, y ~ x2) For instance, we could include the effect of depth for only the positive model, and remove it for the presence-absence model.

However, there are currently limitations if specifying separate formulas for each model component. The two formulas cannot have:

• smoothers
• threshold effects
• random intercepts

For now, these must be specified through a single formula that is shared across the two models.

Each model component can similarly have separate settings for share_range, which determines whether there is a shared spatial and spatiotemporal range parameter (TRUE) or independent range parameters (FALSE), by using a list.

Lastly, whether or not anisotropy is included in the model is determined with the logical argument anisotropy (i.e., TRUE or FALSE), and cannot be separately specified for each model. If anisotropy is included, it is by default shared across the two model components. However it can be made unique in each model component by using sdmTMBcontrol(map = ...) and adding the argument control when fitting the model. This ‘maps’ the anisotropy parameters be unique across model components.

Once we fit the delta model, we can evaluate and plot the output, diagnostics, and predictions similar to other sdmTMB models.

The printed model output will show estimates and standard errors of parameters for each model separately.

print(fit_dg)

Using the tidy() function will turn the sdmTMB model output into a data frame, with the argument model=1 or model=2 to specify which model component to extract as a dataframe. See tidy.sdmTMB() for additional arguments and options.

tidy(fit_dg) # model = 1 is default
tidy(fit_dg, model = 1)
tidy(fit_dg, model = 1, "ran_pars", conf.int = TRUE)
tidy(fit_dg, model = 2)
tidy(fit_dg, model = 2, "ran_pars", conf.int = TRUE)

We can use predict() to predict density in either log space (the default) or transformed into the response scale (type = "response").

For built-in delta models, the default function will return estimated response and parameters for each grid cell for each model separately, notated with a 1 (for the presence/absence model) or 2 (for the positive catch model) in the column name. See predict.sdmTMB() for a description of values in the data frame.

p <- predict(fit_dg, newdata = qcs_grid)
str(p)

p_response <- predict(fit_dg, newdata = qcs_grid, type = "response")

We can use the model argument to return only predictions from model 1 or model 2.

To return combined predictions from both components on the response scale, we can use model=NA

We can use predictions from the built-in delta model (making sure that return_tmb_object=TRUE) to get the index values using the get_index() function. This can be used with predictions that include both the first and second models (i.e., using the default and specifying no model argument) or with predictions generated using model=NA. The get_index() function will automatically combine predictions from the first and second model in calculating the index values. For more on modelling for the purposes of creating an index see the vignette on Index standardization with sdmTMB.

p2 <- predict(fit_dg, newdata = qcs_grid, return_tmb_object = TRUE)
ind_dg <- get_index(p2, bias_correct = FALSE)

We can plot conditional effects of covariates (such as depth in the example model) using the package visreg by specifying each model component with model=1 for the presence-absence model or model=2 for the positive catch model. Currently, plotting effects of built-in delta models with ggeffects is not supported. See the vignette on using visreg with sdmTMB for more information.

visreg_delta(fit_dg, xvar = "depth", model = 1, gg = TRUE)

visreg_delta(fit_dg, xvar = "depth", model = 2, gg = TRUE)

The built-in delta models can also be evaluated with the residuals() functions in sdmTMB. Similarly to generating predictions, we can specify which of the model components we want to return residuals for using the model argument and specifying =1 or =2. See residuals.sdmTMB() for additional options for evaluating residuals in sdmTMB models.

We can also simulate new observations from a fitted delta model. As in other functions, we can specify which model to simulate from using the argument model=1 for only presence/absence, model=2 for only positive catches, or model=NA for combined predictions. See simulate.sdmTMB() for more details on simulation options.

simulations <- simulate(fit_dg, nsim = 5, seed = 5090, model = NA)

Delta models by fitting two components separately and combining predictions

Next, we will show an example of how to implement a delta-gamma model in sdmTMB, but with each component fit separately and then combined. This approach gives maximum flexibility for each model and lets you develop them on at a time. It has limitations if you are calculating an index of abundance or if you want to share some parameters.

glimpse(pcod)

mesh1 <- make_mesh(pcod, c("X", "Y"), cutoff = 20) # coarse for vignette speed

It is not necessary to use the same mesh for both models, but one can do so by updating the first mesh to match the reduced data frame as shown here:

dat2 <- subset(pcod, density > 0)
mesh2 <- make_mesh(dat2,
  xy_cols = c("X", "Y"),
  mesh = mesh1$mesh
)

This delta-gamma model is similar to the Tweedie model in the Intro to modelling with sdmTMB vignette, except that we will use s() for the depth effect.

m1 <- sdmTMB(
  formula = present ~ 0 + as.factor(year) + s(depth, k = 3),
  data = pcod,
  mesh = mesh1,
  time = "year", family = binomial(link = "logit"),
  spatiotemporal = "iid",
  spatial = "on"
)
m1

One can use different covariates in each model, but in this case we will just let the depth effect be more wiggly by not specifying k = 3.

m2 <- sdmTMB(
  formula = density ~ 0 + as.factor(year) + s(depth),
  data = dat2,
  mesh = mesh2,
  time = "year",
  family = Gamma(link = "log"),
  spatiotemporal = "iid",
  spatial = "on"
)
m2

Next, we need some way of combining the predictions across the two models. If all we need are point predictions, we can just multiply the predictions from the two models after applying the inverse link:

pred <- qcs_grid # use the grid as template for saving our predictions
p_bin <- predict(m1, newdata = qcs_grid)
p_pos <- predict(m2, newdata = qcs_grid)
p_bin_prob <- m1$family$linkinv(p_bin$est)
p_pos_exp <- m2$family$linkinv(p_pos$est)
pred$est_exp <- p_bin_prob * p_pos_exp

But if a measure of uncertainty is required, we can simulate from the joint parameter precision matrix using the predict() function with any number of simulations selected (e.g., sims = 500). Because the predictions come from simulated draws from the parameter covariance matrix, the predictions will become more consistent with a larger number of draws. However, a greater number of draws takes longer to calculate and will use more memory (larger matrix), so fewer draws (~100) may be fine for experimentation. A larger number (say ~1000) may be appropriate for final model runs.

set.seed(28239)
p_bin_sim <- predict(m1, newdata = qcs_grid, nsim = 100)
p_pos_sim <- predict(m2, newdata = qcs_grid, nsim = 100)
p_bin_prob_sim <- m1$family$linkinv(p_bin_sim)
p_pos_exp_sim <- m2$family$linkinv(p_pos_sim)
p_combined_sim <- p_bin_prob_sim * p_pos_exp_sim

p_combined_sim is a matrix with a row for each row of data that was predicted on and width nsim. You can process this matrix however you would like. We can save median predictions and upper and lower 95% confidence intervals:

pred$median <- apply(p_combined_sim, 1, median)
plot(pred$est_exp, pred$median)

ggplot(subset(pred, year == 2017), aes(X, Y, fill = median)) +
  geom_raster() +
  coord_fixed() +
  scale_fill_viridis_c(trans = "sqrt")

And we can calculate spatial uncertainty:

pred$cv <- apply(p_combined_sim, 1, function(x) sd(x) / mean(x))
ggplot(subset(pred, year == 2017), aes(X, Y, fill = cv)) + # 2017 as an example
  geom_raster() +
  coord_fixed() +
  scale_fill_viridis_c(trans = "log10")