See the residual-checking vignette: browseVignettes("sdmTMB") or on the documentation site. See notes about types of residuals in 'Details' section below.

# S3 method for sdmTMB
residuals(
  object,
  type = c("mle-laplace", "mle-mcmc", "mvn-laplace", "response"),
  mcmc_iter = 500,
  mcmc_warmup = 250,
  print_stan_model = FALSE,
  stan_args = NULL,
  model = c(1, 2),
  ...
)

Arguments

object

An sdmTMB() model

type

Type of residual. See details.

mcmc_iter

Iterations for MCMC residuals. Will take the last one.

mcmc_warmup

Warmup for MCMC residuals.

print_stan_model

Print the Stan model from MCMC residuals?

stan_args

A list of arguments that will be passed to rstan::sampling().

model

Which delta/hurdle model component?

...

Passed to residual function. Only n works for binomial.

Value

A vector of residuals.

Details

Types of residuals currently supported:

"mle-laplace" refers to randomized quantile residuals (Dunn & Smyth 1996), which are also known as probability integral transform (PIT) residuals (Smith 1985). Under model assumptions, these should be distributed as standard normal with the following caveat: the Laplace approximation used for the latent/random effects can cause these residuals to deviate from the standard normal assumption even if the model is consistent with the data (Thygesen et al. 2017). Therefore, these residuals are fast to calculate but can be unreliable.

"mle-mcmc" refers to randomized quantile residuals where the fixed effects are fixed at their MLE (maximum likelihoood estimate) values and the random effects are sampled with MCMC via tmbstan/Stan. As proposed in Thygesen et al. (2017) and used in Rufener et al. (2021). Under model assumptions, these should be distributed as standard normal. These residuals are theoretically preferred over the regular Laplace approximated randomized-quantile residuals, but will be considerably slower to calculate.

Ideally MCMC is run until convergence and then the last iteration can be used for residuals. MCMC samples are defined by mcmc_iter - mcmc_warmup. The Stan model can be printed with print_stan_model = TRUE to check. The defaults may not be sufficient for many models.

"mvn-laplace" is the same as "mle-laplace" except the parameters are based on simulations drawn from the assumed multivariate normal distribution (using the joint precision matrix).

"response" refers to response residuals: observed minus predicted.

References

Dunn, P.K. & Smyth, G.K. (1996). Randomized Quantile Residuals. Journal of Computational and Graphical Statistics, 5, 236–244.

Smith, J.Q. (1985). Diagnostic checks of non-standard time series models. Journal of Forecasting, 4, 283–291.

Rufener, M.-C., Kristensen, K., Nielsen, J.R., and Bastardie, F. 2021. Bridging the gap between commercial fisheries and survey data to model the spatiotemporal dynamics of marine species. Ecological Applications. e02453. doi:10.1002/eap.2453

Thygesen, U.H., Albertsen, C.M., Berg, C.W., Kristensen, K., and Nielsen, A. 2017. Validation of ecological state space models using the Laplace approximation. Environ Ecol Stat 24(2): 317–339. doi:10.1007/s10651-017-0372-4

Examples

if (inla_installed() &&
    require("tmbstan", quietly = TRUE) &&
    require("rstan", quietly = TRUE)) {

  mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 10)
  fit <- sdmTMB(
    present ~ as.factor(year) + poly(depth, 3),
    data = pcod_2011, mesh = mesh,
    family = binomial()
  )

  # response residuals will be not be normally distributed unless
  # the family is Gaussian:
  r0 <- residuals(fit, type = "response")
  qqnorm(r0)
  qqline(r0)

  # quick but can have issues because of Laplace approximation:
  r1 <- residuals(fit, type = "mle-laplace")
  qqnorm(r1)
  qqline(r1)

  # MCMC-based with fixed effects at MLEs; best but slowest:
  set.seed(2938)
  r2 <- residuals(fit, type = "mle-mcmc", mcmc_iter = 101, mcmc_warmup = 100)
  qqnorm(r2)
  qqline(r2)

  # Example of passing control arguments to rstan::sampling():
  # 11 iterations used for a quick example; don't do this normally
  stan_args <- list(control = list(adapt_delta = 0.9, max_treedepth = 12))
  r3 <- residuals(
    fit, type = "mle-mcmc", mcmc_iter = 11, mcmc_warmup = 10,
    stan_args = stan_args
  )
}


#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001329 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.29 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: There aren't enough warmup iterations to fit the
#> Chain 1:          three stages of adaptation as currently configured.
#> Chain 1:          Reducing each adaptation stage to 15%/75%/10% of
#> Chain 1:          the given number of warmup iterations:
#> Chain 1:            init_buffer = 15
#> Chain 1:            adapt_window = 75
#> Chain 1:            term_buffer = 10
#> Chain 1: 
#> Chain 1: Iteration:   1 / 101 [  0%]  (Warmup)
#> Chain 1: Iteration:  10 / 101 [  9%]  (Warmup)
#> Chain 1: Iteration:  20 / 101 [ 19%]  (Warmup)
#> Chain 1: Iteration:  30 / 101 [ 29%]  (Warmup)
#> Chain 1: Iteration:  40 / 101 [ 39%]  (Warmup)
#> Chain 1: Iteration:  50 / 101 [ 49%]  (Warmup)
#> Chain 1: Iteration:  60 / 101 [ 59%]  (Warmup)
#> Chain 1: Iteration:  70 / 101 [ 69%]  (Warmup)
#> Chain 1: Iteration:  80 / 101 [ 79%]  (Warmup)
#> Chain 1: Iteration:  90 / 101 [ 89%]  (Warmup)
#> Chain 1: Iteration: 100 / 101 [ 99%]  (Warmup)
#> Chain 1: Iteration: 101 / 101 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 3.51462 seconds (Warm-up)
#> Chain 1:                0.013746 seconds (Sampling)
#> Chain 1:                3.52837 seconds (Total)
#> Chain 1: 

#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.001417 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.17 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: WARNING: No variance estimation is
#> Chain 1:          performed for num_warmup < 20
#> Chain 1: 
#> Chain 1: Iteration:  1 / 11 [  9%]  (Warmup)
#> Chain 1: Iteration:  2 / 11 [ 18%]  (Warmup)
#> Chain 1: Iteration:  3 / 11 [ 27%]  (Warmup)
#> Chain 1: Iteration:  4 / 11 [ 36%]  (Warmup)
#> Chain 1: Iteration:  5 / 11 [ 45%]  (Warmup)
#> Chain 1: Iteration:  6 / 11 [ 54%]  (Warmup)
#> Chain 1: Iteration:  7 / 11 [ 63%]  (Warmup)
#> Chain 1: Iteration:  8 / 11 [ 72%]  (Warmup)
#> Chain 1: Iteration:  9 / 11 [ 81%]  (Warmup)
#> Chain 1: Iteration: 10 / 11 [ 90%]  (Warmup)
#> Chain 1: Iteration: 11 / 11 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.301117 seconds (Warm-up)
#> Chain 1:                0.028479 seconds (Sampling)
#> Chain 1:                0.329596 seconds (Total)
#> Chain 1: