Can I get p-values out of a fitted sdmTMB object

[ericward-noaa]

Question via email:

When I use summary() on a glmmTMB object, it returns p-values. But with sdmTMB we only get parameter means and SEs. Can similar functionality be included?

[ericward-noaa]

Currently, this isn't included in sdmTMB. It's unclear if this will be supported in the future, but for fixed effects, it's not terribly difficult to do this yourself with a fitted model, using the same approach taken with glmmTMB. These calculations are done in glmmTMB when the summary function is called here, https://github.com/glmmTMB/glmmTMB/blob/ba8575b22fd5b0b513fe08e5f48b95e94a3eabcc/glmmTMB/R/glmmTMB.R#L1641

Using the basic Pacific cod dataset, we'll construct a model with year as a factor and just a spatial field.

library(sdmTMB)

pcod$fyear <- as.factor(pcod$year)
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
fit <- sdmTMB(
  density ~ fyear,
  data = pcod, mesh = mesh,
  family = tweedie(link = "log")
)

The coefficients and SEs are summarized with tidy(), and the z-statistics are just the standardized coefficients

pars <- tidy(fit)

# calculate z-statistics
pars[["Z stat"]] <- pars$estimate/pars$std.error

Next we can calculate the probability of observing a test statistic more extreme than these values, using the standard normal CDF like glmmTMB

pars[[paste0("Pr(>|","z","|)")]] <- 2*pnorm(abs(pars$`Z stat`), 
                                            lower.tail = FALSE)

knitr::kable(pars, digits = 3)

[seananderson]

Or more directly with TMB::summary.sdreport(), but without the nicely named coefficients (look at the b_j in order and they should match print() or tidy().

library(sdmTMB)
pcod$fyear <- as.factor(pcod$year)
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 20)
fit <- sdmTMB(
  density ~ fyear,
  data = pcod, mesh = mesh,
  family = tweedie(link = "log")
)
summary(fit$sd_report, select = "fixed", p.value = TRUE)
#>             Estimate Std. Error    z value   Pr(>|z^2|)
#> b_j       2.86631816 0.33468294  8.5642793 1.087539e-17
#> b_j       0.71727094 0.19587369  3.6619056 2.503461e-04
#> b_j       0.66425861 0.19211819  3.4575520 5.451071e-04
#> b_j      -0.66291876 0.21140993 -3.1357031 1.714427e-03
#> b_j      -0.28448720 0.20808495 -1.3671686 1.715725e-01
#> b_j       0.39562364 0.19239170  2.0563447 3.974931e-02
#> b_j       0.28809870 0.19372977  1.4871163 1.369841e-01
#> b_j       0.49228332 0.19356839  2.5432010 1.098420e-02
#> b_j      -0.09669499 0.20590181 -0.4696170 6.386287e-01
#> ln_tau_O -1.17096807 2.39931937 -0.4880418 6.255203e-01
#> ln_kappa -1.47354854 1.30478326 -1.1293435 2.587529e-01
#> thetaf    0.46392655 0.04186171 11.0823601 1.527910e-28
#> ln_phi    2.70794406 0.03052901 88.7006864 0.000000e+00

^{Created on 2022-08-30 with reprex v2.0.2}

[ericward-noaa]

This is better than doing it manually!