search

pbs-assess / sdmTMB

:earth_americas: An R package for spatial and spatiotemporal GLMMs with TMB
https://pbs-assess.github.io/sdmTMB/
182 stars 26 forks source link

`mvn-laplace` residual estimation fails due to `predict.sdmTMB()` offset issue #326

Closed sebpardo closed 4 months ago

sebpardo commented 4 months ago

Estimating residuals using type = "mvn-laplace" throws an error due to problems with extra_time accounting in the offset.

library(sdmTMB)

mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 10)

pcod$os <- rep(log(0.01), nrow(pcod)) # offset

m <- sdmTMB(
  data = pcod,
  formula = density ~ poly(log(depth), 2),
  mesh = mesh,
  offset = pcod$os,
  family = tweedie(link = "log"),
  spatial = "on",
  time = "year",
  extra_time = c(2006, 2008, 2010, 2012, 2014, 2016), 
  spatiotemporal = "AR1"
)

resids_mvn <- residuals(m, type = "mvn-laplace")

The culprit in this situation seems to be here, where in this case the vector offset already has the extra_time values appended to it before this line, so extra_time is being added twice.

https://github.com/pbs-assess/sdmTMB/blob/c01b13adb9da13a2ebde91e8f71331d4b4725db3/R/predict.R#L386

Output from browser() inside predict.sdmTMB():

Browse[5]> dim(pcod)
[1] 2143   13
Browse[5]> dim(newdata)
[1] 2149   19
Browse[5]> nrow(proj_X_ij[[1]])
[1] 2149
Browse[5]> length(offset)
[1] 2155