How can I report the significance tests of p-values for factor variables with sdmTMB?

[ericward-noaa]

Question via email: I have several factors in a model, and I'd like to try to report the significance of each. sdmTMB objects aren't compatible with anova() -- but is there any other way I might be able to do this?

[ericward-noaa]

I think we can extact the log-likelihoods and degrees of freedom from fitted models, and generate p-values manually. Example with the pcod_2011 dataset below,

set.seed(123)
library(sdmTMB)
#> Warning: package 'sdmTMB' was built under R version 4.3.2
# Make the mesh
mesh <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 20)

# Create 2 dummy factor levels for the pcod_2011 dataset
pcod_2011$dummy1 <- sample(c("A","B"),size=nrow(pcod_2011), replace=T)
pcod_2011$dummy2 <- sample(c("C","D"),size=nrow(pcod_2011), replace=T)

# Full model -- additive effects of each, no interaction
fit <- sdmTMB(
  density ~ dummy1 + dummy2,
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)

# Fit the null model -- no covariates
null_model <- sdmTMB(density ~ 1, 
                     data = pcod_2011, 
                     , mesh = mesh,
                     family = tweedie(link = "log"))

# Perform Likelihood Ratio Test for the overall significance
# This is kind of a pain, but we can extract likelihoods manually
ll_full <- logLik(fit)
ll_null <- logLik(null_model)

# Calculate the test statistic
LRT_stat <- -2 * (ll_null - ll_full)

# Degrees of freedom is the difference in the number of parameters
df <- attr(ll_full, "df") - attr(ll_null, "df")

# Calculate the p-value
p_value <- pchisq(LRT_stat, df, lower.tail = FALSE)
print(p_value)
#> 'log Lik.' 0.08819866 (df=5)

# If we were interested in testing the significance of an individual
# factor, we could also do that. For example, to test the significance
# of factor1, we fit the model without that factor and recalculate the p-value
fit2 <- sdmTMB(
  density ~ dummy2,
  data = pcod_2011, mesh = mesh,
  family = tweedie(link = "log")
)
ll_full <- logLik(fit2)
ll_null <- logLik(null_model)

# Calculate the test statistic
LRT_stat <- -2 * (ll_null - ll_full)
# Degrees of freedom is the difference in the number of parameters
df <- attr(ll_full, "df") - attr(ll_null, "df")
p_value <- pchisq(LRT_stat, df, lower.tail = FALSE)
print(p_value)
#> 'log Lik.' 0.7903163 (df=5)

^{Created on 2024-08-02 with reprex v2.1.0}

[seananderson]

Also see: summary(fit$sd_report, p.value=TRUE) Similar with as.list() see as.list.sdreport https://rdrr.io/cran/TMB/man/as.list.sdreport.html And these methods: https://pbs-assess.github.io/sdmTMB/reference/emmeans.sdmTMB.html https://pbs-assess.github.io/sdmTMB/reference/Effect.sdmTMB.html