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

Fix mixture scaling for mixture families #225

Closed ericward-noaa closed 1 year ago

ericward-noaa commented 1 year ago

@seananderson -- not sure if you had other ways of thinking about this, but I made a change over on the mix-predict branch.

Here's a reproducible example showing more similar results to Kelli / Ian's example:

pcod_spde <- make_mesh(pcod, c("X", "Y"), n_knots = 60, type = "kmeans")
m <- sdmTMB(
  data = pcod,
  formula = density ~ 0 + as.factor(year),
  time = "year", mesh = pcod_spde, family = tweedie(link = "log")
)

nd <- replicate_df(qcs_grid, "year", unique(pcod$year))
predictions <- predict(m, newdata = nd, return_tmb_object = TRUE)
ind1 <- get_index(predictions)

m <- sdmTMB(
  data = pcod,
  formula = density ~ 0 + as.factor(year),
  time = "year", mesh = pcod_spde, family = delta_lognormal_mix()
)
predictions <- predict(m, newdata = nd, return_tmb_object = TRUE)
ind2 <- get_index(predictions)

ggplot(ind1, aes(year, est)) + 
  geom_point(col="blue") + 
  geom_point(data=ind2, aes(year, est), col="red")
seananderson commented 1 year ago

Thanks. I was going to condense the C++ case switches and pull out common parts in the observation likelihood part. Otherwise, this seems good. I’m not sure the residuals are correct right now. Can we just take the combined mean? Wouldn’t we expect that it could be bimodal? We could simulation test to check. Simulation style (dharma) residuals would be simple.

ericward-noaa commented 1 year ago

Thanks, I wasn't sure about best way to reduce the switches, but that'd be great.

I agree the residuals will be problematic -- if the data are at all multimodal (because of the mixture) and we use the combined mean, residuals should be multimodal. If we were calculating the probability of group assignment for each data point, then I think each point could have its own mean (and residuals would be well behaved) -- I think perhaps putting a warning in residuals() for models that get passed in with these families would be ok for now?

kellijohnson-NOAA commented 1 year ago

Sorry my tardiness on this issue, I am looking for the mix-predict branch and I don't see it.

seananderson commented 1 year ago

I think this is now fixed in the main branch. Residuals are not fixed, and I think will just error out right now (although DHARMa ones should be fine... or at least as fine as DHARMa residuals are for these types of models).

seananderson commented 1 year ago

For future reference, it's fixed mostly in this commit. This still isn't widely tested, so be extra critical of the results and let us know if something looks funny.

kellijohnson-NOAA commented 1 year ago

Will do, thanks @seananderson and @ericward-noaa!