DHARMa residuals with delta models

[seananderson]

Question by email:

I was wondering is it possible to use the DHARMa tools when running a delta-gamma model? It seems like it should be but I can’t seem to figure out how to specify first the binomial and then the gamma.

[seananderson]

They should just work for the full delta-gamma model prediction. If you want to look at the 2 parts separately, you'll need to build the DHARMa object yourself though because dharma_residuals() always works with the full prediction. See below.

Also, note these DHARMa residuals for spatial models can often show patterns even when everything is just fine. The MCMC residuals are probably best (residuals(fit, type = "mle-mcmc")) see ?residuals.sdmTMB and https://pbs-assess.github.io/sdmTMB/articles/residual-checking.html but they can be slow to calculate.

library(sdmTMB)
pcod_spde <- make_mesh(pcod, c("X", "Y"), cutoff = 15)

# Look at using dharma_residuals on aggregate
fit <- sdmTMB(density ~ 1,
  data = pcod, mesh = pcod_spde, spatial = "on",
  spatiotemporal = list("off", "off"), time = "year",
  family = delta_gamma()
)
p <- predict(fit, type = "response")

s <- simulate(fit, nsim = 200)
dharma_residuals(s, fit)

# Split the fit into looking at residuals for each submodel -- first presence
s1 <- simulate(fit, nsim = 200, model = 1)

res <- DHARMa::createDHARMa(
  simulatedResponse = s1,
  observedResponse = pcod$present,
  fittedPredictedResponse = p$est1
)
plot(res)

# Second positive
s2 <- simulate(fit, nsim = 200, model = 2)

pos <- pcod$density > 0

res <- DHARMa::createDHARMa(
  simulatedResponse = s2[pos, ],
  observedResponse = pcod$density[pos],
  fittedPredictedResponse = p$est2[pos]
)
plot(res)

^{Created on 2022-08-23 with reprex v2.0.2}