pbs-assess / sdmTMB

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

combined predictions from delta/hurdle model not working as described #160

Closed ecophilina closed 1 year ago

ecophilina commented 1 year ago

predict(fit_delta, model = NA) does not give the output described for the function:

"model = Type of prediction if a delta/hurdle model: NA returns the combined prediction from both components on the link scale for the positive component; 1 or 2 return the first or second model component only on the link or response scale depending on the argument type."

Instead it reports the same way as predict(fit_delta).

predict(fit_delta, model = NA, se_fit = TRUE)

also returns the same dataframe as

predict(fit_delta, model = 1, se_fit = TRUE)

Does model = NA only apply when using nsim?

maxlindmark commented 1 year ago

I was just about to post a question on the scale of the combined prediction for a delta model and saw this. Now I wonder if these issues are related and am unsure if I should make a separate post or comment here (let me know if you'd prefer that!).

This code shows that when I specify type="link" and type="response", the separate model components do get different scales, but not the total est

pcod_mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 15)

fit_dg <- sdmTMB(density ~ 1 + s(depth),
                 data = pcod,
                 mesh = pcod_mesh,
                 time = "year",
                 spatial = "off", 
                 spatiotemporal = "off",
                 family = delta_gamma()
)

nd <- data.frame(
  depth = seq(min(pcod$depth),
              max(pcod$depth),
              length.out = 20
  ),
  X = 0,
  Y = 0,
  year = 2015L # a chosen year
)

p1 <- predict(fit_dg,
              newdata = nd,
              re_form = NULL,
              re_form_iid = NULL,
              se_fit = TRUE,
              type = "link")

p2 <- predict(fit_dg,
              newdata = nd,
              re_form = NULL,
              re_form_iid = NULL,
              se_fit = TRUE,
              type = "response")

head(p1 %>% dplyr::select(est1, est2, est), 2)
> head(p1 %>% dplyr::select(est1, est2, est), 2)
       est1     est2       est
1 -2.582749 3.973598 -2.582749
2 -1.241710 4.118488 -1.241710

head(p2 %>% dplyr::select(est1, est2, est), 2)
> head(p2 %>% dplyr::select(est1, est2, est), 2)
        est1     est2       est
1 0.07025697 53.17552 -2.582749
2 0.22413842 61.46625 -1.241710

P.S. Thanks for an amazing package!

seananderson commented 1 year ago

I checked: yes it only affects nsim > 0 or a Bayesian model. I've updated the docs.

@maxlindmark I'll link to your example over in the issue on type = link/response.