predict() function issues

[ericward-noaa]

1. For delta model, filter out NAs for positive component internally

2. Model argument doesn't seem to be changing anything for delta models -- predictions for both components are returned. I think the fix here is just drop 'model' as an argument

3. When re_form = NA, predict() doesn't respond to the "type" argument and returns identical predictions. Investigating this with delta models, but same applies to non-delta

[seananderson]

1. was intentional so the two parts always line up
2. was only affecting nsim > 0 or Stan models; fixed the docs
3. this one is complex and a mess. A few questions in the delta version below (also need to test with se_fit = TRUE and nsim > 0). The complexity of this makes me wonder if we should rethink how this whole chunk is coded.

library(sdmTMB)
fit <- sdmTMB(
  density ~ 1,
  family = tweedie(),
  data = pcod_2011, 
  mesh = pcod_mesh_2011
)

p <- predict(fit, type = "link")
mean(p$est)
#> [1] 3.10099
# OK

p <- predict(fit, type = "link", re_form = NA)
mean(p$est)
#> [1] 3.001329
# OK

p <- predict(fit, type = "response")
mean(p$est)
#> [1] 38.18027
# OK

p <- predict(fit, type = "response", re_form = NA)
mean(p$est)
#> [1] 3.001329
# should be fixed to match

fit_delt <- sdmTMB(
  density ~ 1,
  family = delta_gamma(),
  data = pcod_2011, 
  mesh = pcod_mesh_2011
)

p <- predict(fit_delt, type = "link")
mean(p$est) 
#> Warning: Unknown or uninitialised column: `est`.
#> Warning in mean.default(p$est): argument is not numeric or logical: returning NA
#> [1] NA
# NA because there is no 'est' in link space (but we could inverse link it)

p <- predict(fit_delt, type = "link", re_form = NA)
mean(p$est)
#> [1] -0.1672364
# responding to `model` arg, which is getting set to 1 if not visreg; fix

p <- predict(fit_delt, type = "response")
mean(p$est)
#> [1] 37.31266
# OK

p <- predict(fit_delt, type = "response", re_form = NA)
mean(p$est)
#> [1] -0.1672364
# should be fixed to match

^{Created on 2022-10-28 with reprex v2.0.2}

[seananderson]

Patching the code that's currently there, how does this logic look?

library(sdmTMB)
fit <- sdmTMB(
  density ~ 1,
  family = tweedie(),
  data = pcod_2011,
  mesh = pcod_mesh_2011
)

p <- predict(fit, type = "link")
mean(p$est)
#> [1] 3.10099

p <- predict(fit, type = "link", re_form = NA)
mean(p$est)
#> [1] 3.001329

p <- predict(fit, type = "response")
mean(p$est)
#> [1] 38.18027

p <- predict(fit, type = "response", re_form = NA)
mean(p$est)
#> [1] 20.11224

fit_delt <- sdmTMB(
  density ~ 1,
  family = delta_gamma(),
  data = pcod_2011,
  mesh = pcod_mesh_2011
)

p <- predict(fit_delt, type = "link")
mean(p$est)
#> Warning: Unknown or uninitialised column: `est`.
#> Warning in mean.default(p$est): argument is not numeric or logical: returning NA
#> [1] NA
mean(p$est1)
#> [1] -0.2229376
mean(p$est2)
#> [1] 4.261005

p <- predict(fit_delt, type = "link", re_form = NA)
mean(p$est)
#> Warning: Unknown or uninitialised column: `est`.
#> argument is not numeric or logical: returning NA
#> [1] NA
mean(p$est1)
#> [1] -0.1672364
mean(p$est2)
#> [1] 4.282224

p <- predict(fit_delt, type = "response")
mean(p$est)
#> [1] 37.31266
mean(p$est1)
#> [1] 0.4602225
mean(p$est2)
#> [1] 77.35623

p <- predict(fit_delt, type = "response", re_form = NA)
mean(p$est)
#> [1] 33.18066
mean(p$est1)
#> [1] 0.4582881
mean(p$est2)
#> [1] 72.40131

^{Created on 2022-10-28 with reprex v2.0.2}

[ericward-noaa]

Yep - looks good