Residual checking with sdmTMB

If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under ‘Articles’.

library(sdmTMB)

We will start with some data simulated from scratch. We will simulate from an NB2 negative binomial observation model, a spatial random field, an intercept, and one predictor named ‘a1’ that will have a linear effect on the observed data.

set.seed(1)
predictor_dat <- data.frame(X = runif(1000), Y = runif(1000), a1 = rnorm(1000))
mesh <- make_mesh(predictor_dat, xy_cols = c("X", "Y"), cutoff = 0.1)
dat <- sdmTMB_simulate(
  formula = ~ 1 + a1,
  data = predictor_dat,
  mesh = mesh,
  family = nbinom2(link = "log"),
  phi = 0.2,
  range = 0.4,
  sigma_O = 0.3,
  seed = 1,
  B = c(0.2, 0.5) # B0 = intercept, B1 = a1 slope
)

Next, we will fit versions with various responses and predictors. The first model will use the Poisson instead of the NB2. The 2nd model will match the simulated data. The third model is missing the ‘a1’ predictor. We’ll use a PC prior on the Matérn parameters to aid in estimation.

pc <- pc_matern(range_gt = 0.1, sigma_lt = 1)

fit_pois <- sdmTMB(observed ~ 1 + a1, data = dat, family = poisson(), mesh = mesh,
  priors = sdmTMBpriors(matern_s = pc))
fit_pois

fit_nb2 <- sdmTMB(observed ~ 1 + a1, data = dat, family = nbinom2(), mesh = mesh,
  priors = sdmTMBpriors(matern_s = pc))
fit_nb2

fit_nb2_miss <- sdmTMB(observed ~ 1, data = dat, family = nbinom2(), mesh = mesh,
  priors = sdmTMBpriors(matern_s = pc))
fit_nb2_miss

We can see just by looking at these fits that the Poisson model inflates the spatial random field standard deviation (SD) compared to the truth. The model missing the ‘a1’ predictor does so to a lesser degree.

Here are randomized quantile residuals at fixed effect MLEs (Maximum Likelihood Estimates) and random effects that maximize the log likelihood at estimated fixed effects:

rq_res <- residuals(fit_pois)
rq_res <- rq_res[is.finite(rq_res)] # some Inf
qqnorm(rq_res);qqline(rq_res)

rq_res <- residuals(fit_nb2)
qqnorm(rq_res);qqline(rq_res)

These use the approach from Dunn and Smyth (1996). They are also known as PIT (probability-integral-transform) residuals. They apply randomization to integer response values, transform the residuals using the distribution function (e.g., pnorm()), simulate from a uniform distribution, and transform the samples such that they would be Gaussian if consistent with the model. You can see the source code at https://github.com/pbs-assess/sdmTMB/blob/master/R/residuals.R

We can see here that there are likely issues with the Poisson model in the tails.

These types of residuals are known to have statistical issues for state-space models; even if the model is the ‘correct’ model, the QQ plot may appear to have problems (Thygesen et al. 2017).

One-step-ahead residuals (Thygesen et al. 2017) are one option to fix this problem (although slow to calculate). Another option is to take a draw from the posterior with MCMC (e.g., Rufener et al. 2021). Also see https://kaskr.github.io/adcomp/_book/Validation.html

Here we will draw MCMC predictions and calculate residuals. The fixed effects will be fixed at their maximum likelihood estimates (MLE) and the random effects will be sampled. We will only take a single draw for speed:

mcmc_res <- residuals(fit_nb2, type = "mle-mcmc", mcmc_iter = 201, mcmc_warmup = 200)
qqnorm(mcmc_res);qqline(mcmc_res)

We can see these look a bit better. Remember, this is the ‘correct’ model.

We can take simulations from the fitted model to use with simulation-based residuals:

s_pois <- simulate(fit_pois, nsim = 500)
s_nb2_miss <- simulate(fit_nb2_miss, nsim = 500)
s_nb2 <- simulate(fit_nb2, nsim = 500)

These return a matrix where each row represents a row of data and each column is a simulation draw:

dim(s_pois)

Test whether fitted models are consistent with the observed number of zeros:

sum(dat$observed == 0) / length(dat$observed)
sum(s_pois == 0)/length(s_pois)
sum(s_nb2 == 0)/length(s_nb2)

There are obviously too few zeros in the data simulated from the Poisson model.

Plot DHARMa residuals:

simulate(fit_pois, nsim = 300) %>% 
  dharma_residuals(fit_pois)

We could also do that manually, which lets us use other DHARMa tools:

# My reading of DHARMa documation is that the predicted response for the 
# residuals vs. fitted plot should ideally not include the random effects:
pred_fixed <- fit_pois$family$linkinv(predict(fit_pois)$est_non_rf)
r_pois <- DHARMa::createDHARMa(
  simulatedResponse = s_pois,
  observedResponse = dat$observed,
  fittedPredictedResponse = pred_fixed
)
plot(r_pois)
DHARMa::testResiduals(r_pois)
DHARMa::testSpatialAutocorrelation(r_pois, x = dat$X, y = dat$Y)
DHARMa::testZeroInflation(r_pois)

In the QQ residual plots we clearly see evidence of over dispersion compared to the Poisson. Note the values clumping near 1.0 on the observed axis and deviating downwards towards 0.0 observed. This is indicative of too many zeros and especially too many large values compared to the assumed Poisson distribution.

Lets try with the correct model:

simulate(fit_nb2, nsim = 300) %>% 
  dharma_residuals(fit_nb2)

Everything looks fine. But, again, the MCMC-based residuals above are likely the best approach.

What about the model where we were missing a predictor?

pred_fixed <- fit_nb2_miss$family$linkinv(predict(fit_nb2_miss)$est_non_rf)
r_nb2_miss <- DHARMa::createDHARMa(
  simulatedResponse = s_nb2_miss,
  observedResponse = dat$observed,
  fittedPredictedResponse = pred_fixed
)
plot(r_nb2_miss)

This looks fine so far, but the plot on the right represents simulated residuals against the prediction without the random effects, which here is just an intercept. Lets try plotting the residuals against the missing predictor:

DHARMa::plotResiduals(r_nb2_miss, form = dat$a1)

We can see a slight trend in the residuals against ‘a1’ since we have missed including it in the model.

We can also see the difference in the log likelihood or by using the AIC() method:

# negative log likelihood is lower; 
# i.e. log likelihood is higher, but we do have one more parameter
fit_nb2$model$objective
fit_nb2_miss$model$objective
AIC(fit_nb2_miss, fit_nb2) # AIC supports including the 'a1' predictor

The above used simulations with the parameters fixed at their Maximum Likelihood Estimate (MLE) and predictions conditional on the fitted random effects. Alternatively, we could simulate with the parameters drawn from their joint precision matrix to encapsulate uncertainty about the parameters. This may be a better test for residual analysis, but this is an open area of research as far as I can tell.

# simulate with the parameters drawn from the joint precision matrix:
s2 <- simulate(fit_nb2, nsim = 1, params = "MVN")

Or we could simulate with new random fields based on the estimated parameters governing the random fields (range and SD):

# simulate with new random fields:
s3 <- simulate(fit_nb2, nsim = 1, re_form = ~ 0)

We could, of course, combine those two options:

# simulate with new random fields and new parameter draws:
s4 <- simulate(fit_nb2, nsim = 500, params = "MVN", re_form = ~ 0)
pred_fixed <- fit_nb2$family$linkinv(predict(fit_nb2)$est_non_rf)
r_nb2 <- DHARMa::createDHARMa(
  simulatedResponse = s4,
  observedResponse = dat$observed,
  fittedPredictedResponse = pred_fixed
)
plot(r_nb2)

These also look OK.

For help interpreting the DHARMa residual plots, see vignette("DHARMa", package="DHARMa").