MCMC residuals error about length(mcmc_samples)

[seananderson]

I received an email about the following error:

library(sdmTMB)
pcod_mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)
fit_dg <- sdmTMB(
  density ~ s(depth),
  data = pcod,
  mesh = pcod_mesh,
  family = delta_gamma()
)
r <- residuals(fit_dg, type = "mle-mcmc", mcmc_iter = 101, mcmc_warmup = 100)

#> Error: length(mcmc_samples) not equal to nrow(object$data)

This will happen on the CRAN 0.3.0 version of sdmTMB and is due to that functionality moving to the sdmTMBextra package and a not-very-helpful error message. I'll expand in an answer below.

[seananderson]

This is due to a recent change we had to make to get back on CRAN. See the residuals vignette for working code:
https://pbs-assess.github.io/sdmTMB/articles/web_only/residual-checking.html

Search for "mcmc".

I just added a check for this to provide a more useful warning.

Basically, install the sdmTMBextra package and use sdmTMBextra::predict_mle_mcmc() first. The rstan/MCMC part was causing memory issues on some obscure CRAN checks that I could not solve otherwise.

Also note that for delta models, the default is the first component of the model. There’s also a model = 2 option (in this case in the predict_mle_mcmc() call). I have just added a message about that if you don't specify the model argument.

Here's a full example, first with a binomial model:

library(sdmTMB)
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()
)

# MCMC-based with fixed effects at MLEs; best but can be slow:
set.seed(2938)
samp <- sdmTMBextra::predict_mle_mcmc(fit, mcmc_iter = 400, mcmc_warmup = 200)
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000771 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 7.71 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 400 [  0%]  (Warmup)
#> Chain 1: Iteration:  40 / 400 [ 10%]  (Warmup)
#> Chain 1: Iteration:  80 / 400 [ 20%]  (Warmup)
#> Chain 1: Iteration: 120 / 400 [ 30%]  (Warmup)
#> Chain 1: Iteration: 160 / 400 [ 40%]  (Warmup)
#> Chain 1: Iteration: 200 / 400 [ 50%]  (Warmup)
#> Chain 1: Iteration: 201 / 400 [ 50%]  (Sampling)
#> Chain 1: Iteration: 240 / 400 [ 60%]  (Sampling)
#> Chain 1: Iteration: 280 / 400 [ 70%]  (Sampling)
#> Chain 1: Iteration: 320 / 400 [ 80%]  (Sampling)
#> Chain 1: Iteration: 360 / 400 [ 90%]  (Sampling)
#> Chain 1: Iteration: 400 / 400 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.91997 seconds (Warm-up)
#> Chain 1:                1.5093 seconds (Sampling)
#> Chain 1:                4.42927 seconds (Total)
#> Chain 1:
#> Some Stan warnings...
r <- residuals(fit, type = "mle-mcmc", mcmc_samples = samp)
qqnorm(r)
qqline(r)

^{Created on 2023-02-14 with reprex v2.0.2}

Delta-gamma example:

library(sdmTMB)
fit_dg <- sdmTMB(
  density ~ 1,
  data = pcod_2011,
  mesh = pcod_mesh_2011,
  family = delta_gamma()
)
set.seed(123)
p1 <- sdmTMBextra::predict_mle_mcmc(fit_dg,
  mcmc_iter = 201, mcmc_warmup = 200,
  model = 1
)
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000506 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.06 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.4385 seconds (Warm-up)
#> Chain 1:                0.00509 seconds (Sampling)
#> Chain 1:                2.44359 seconds (Total)
#> Chain 1:
set.seed(123)
p2 <- sdmTMBextra::predict_mle_mcmc(fit_dg,
  mcmc_iter = 201, mcmc_warmup = 200,
  model = 2
)
#> 
#> SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000498 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.98 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
#> Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
#> Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
#> Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
#> Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
#> Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
#> Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
#> Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
#> Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
#> Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
#> Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
#> Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 2.42512 seconds (Warm-up)
#> Chain 1:                0.005024 seconds (Sampling)
#> Chain 1:                2.43014 seconds (Total)
#> Chain 1:

set.seed(42)
r1 <- residuals(fit_dg, type = "mle-mcmc", mcmc_samples = p1)
set.seed(42)
r2 <- residuals(fit_dg, type = "mle-mcmc", mcmc_samples = p2)

^{Created on 2023-02-14 with reprex v2.0.2}