pbs-assess / sdmTMB

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

Delta model with spatially varying factor predictor and no spatiotemporal field returns error #237

Closed maxlindmark closed 1 year ago

maxlindmark commented 1 year ago

Hi! I've been playing around with models that have a spatially varying factor predictor (from the zeta-intercept branch) and it's worked really well. Until I tried a delta model with spatiotemporal="off". Here's an example:

# Using dev branch of sdmTMB
# with_libpaths(new = "/Users/maxlindmark/Dropbox/Max work/R/sdmTMB-zeta-intercept/",
#               install_github("pbs-assess/sdmTMB", ref = "zeta-intercept"))

library(sdmTMB, lib.loc = "/Users/maxlindmark/Dropbox/Max work/R/sdmTMB-zeta-intercept")

# Add in fake quarter to mimic the simple example in my data
pcod_q2 <- pcod
pcod_q1 <- pcod

pcod_q1$quarter <- as.factor(1)
pcod_q2$quarter <- as.factor(2)

pcod_q2$density <- pcod_q2$density + rnorm(10, 20, n = nrow(pcod)) # just adding some difference between quarters..

pcod2 <- rbind(pcod_q1, pcod_q2)

# Fit delta model with spatially varying quarter effect
mesh <- make_mesh(pcod2, c("X", "Y"), cutoff = 15)

m <- sdmTMB(density ~ 0 + as.factor(year) + quarter,
            data = pcod2,
            mesh = mesh,
            family = delta_gamma(link1 = "logit", link2 = "log"),
            spatiotemporal = "off",
            spatial = "off", # since spatially varying predictor is a factor
            spatial_varying = ~0 + quarter,
            time = "year"
            )

This model returns the following error:

Error in if (spatiotemporal[m] == "off" && spatial[m] == "off") tmb_map$ln_kappa[,  : 
  missing value where TRUE/FALSE needed

I don't get this error if I change to tweedie or if I add a spatiotemporal field (e.g., IID). I understand both model components need to share the same structure for spatial_varying, but I'm not sure this is intended.

Thanks for the help! Max

seananderson commented 1 year ago

Thanks for the bug report and reproducible example! I've pushed a fix and added a unit test. Note that the 'zeta-intercept' branch was merged into the main branch a while back so there's no need to use it now.

maxlindmark commented 1 year ago

Ah I see. Amazing, thanks a lot!