How do you obtain and display a spline in the sdmTMB model object?

[dpanzeri]

I have a model like this:
fit <- sdmTMB(
N_km2 ~ 0+ s(X.utm)+s(depth)+s(tmp_bot)+s(oxy)+s(npp)+fyear,
data = my.data,
mesh = mesh,
time = "YEAR",
family = delta_gamma(link1='logit',link2 = "log"),
spatial = "on",
spatiotemporal = 'iid'
)

the visreg(fit, "depth") doesn't work.... (error: Error in formula.default(object, env = baseenv()) )
How can plot my spline? Or to obtain spline values and display a spline?

[seananderson]

Here's an example below.

Note the example is about the 'delta' or hurdle families. For non-delta families you can skip the visreg_delta() part and use visreg::visreg() directly and skip set_delta_model() and use ggeffects::ggpredict() directly.

Also note that these examples require the latest sdmTMB (0.4.1 or the GitHub version).

library(sdmTMB)
library(ggplot2)

fit_dg <- sdmTMB(
  density ~ s(depth_scaled, k = 3),
  data = pcod_2011,
  spatial = "off",
  family = delta_gamma()
)

# with visreg:
visreg_delta(fit_dg, xvar = "depth_scaled", model = 1, gg = TRUE)
#> These are residuals for delta model component 1. Use the `model` argument to
#> select the other component.

visreg_delta(fit_dg, xvar = "depth_scaled", model = 2, gg = TRUE)
#> These are residuals for delta model component 2. Use the `model` argument to
#> select the other component.

# with the latest ggeffects from GitHub:
# https://github.com/strengejacke/ggeffects/
# nicest option going forward!?
# remotes::install_github("strengejacke/ggeffects")

# combined:
set_delta_model(fit_dg, model = NA) |> 
  ggeffects::ggpredict("depth_scaled [all]") |>
  plot()

# 1st part:
set_delta_model(fit_dg, model = 1) |> 
  ggeffects::ggpredict("depth_scaled [all]") |>
  plot()

# 2nd part:
set_delta_model(fit_dg, model = 2) |> 
  ggeffects::ggpredict("depth_scaled [all]") |>
  plot()

# by hand:
nd <- data.frame(depth_scaled = seq(
  min(pcod_2011$depth_scaled), 
  max(pcod_2011$depth_scaled), length.out = 100)
)

# combined:
p <- predict(fit_dg, newdata = nd, se_fit = TRUE, re_form = NA, model = NA)
ggplot(p, aes(depth_scaled, exp(est),
  ymin = exp(est - 2 * est_se), 
  ymax = exp(est + 2 * est_se))) + 
  geom_line() +
  geom_ribbon(alpha = 0.5)

# 1st part:
p <- predict(fit_dg, newdata = nd, se_fit = TRUE, re_form = NA, model = 1)
ggplot(p, aes(depth_scaled, plogis(est), 
  ymin = plogis(est - 2 * est_se), 
  ymax = plogis(est + 2 * est_se))) + 
  geom_line() +
  geom_ribbon(alpha = 0.5)

# 2nd part:
p <- predict(fit_dg, newdata = nd, se_fit = TRUE, re_form = NA, model = 2)
ggplot(p, aes(depth_scaled, exp(est), 
  ymin = exp(est - 2 * est_se), 
  ymax = exp(est + 2 * est_se))) + 
  geom_line() +
  geom_ribbon(alpha = 0.5)

^{Created on 2023-11-03 with reprex v2.0.2}

[dpanzeri]

Thank you. Not all the code works (ggeffects does not work), but 'visreg' does.
And if I want to understand the significance of each covariate and the relative variance explained by each of the covariates?
I can use summary in mgcv, but not with sdmTMB.

[seananderson]

Re: ggeffects, yes, there was a note that for now you'll need the latest GitHub version until ggeffects next gets updated on CRAN:

remotes::install_github("strengejacke/ggeffects")

Re: significance of the smooths, no, but there are 2 places you can look for related information:

library(sdmTMB)

fit <- sdmTMB(
  density ~ s(depth_scaled),
  data = pcod,
  spatial = "off",
  family = tweedie()
)

fit
#> Model fit by ML ['sdmTMB']
#> Formula: density ~ s(depth_scaled)
#> Mesh: NULL (isotropic covariance)
#> Data: pcod
#> Family: tweedie(link = 'log')
#>  
#>               coef.est coef.se
#> (Intercept)       2.93    0.06
#> sdepth_scaled    -2.18    1.21
#> 
#> Smooth terms:
#>                   Std. Dev.
#> sds(depth_scaled)      6.47
#> 
#> Dispersion parameter: 14.85
#> Tweedie p: 1.62
#> ML criterion at convergence: 6553.533

^{Created on 2023-11-04 with reprex v2.0.2}

This part is the linear component:

              coef.est coef.se
sdepth_scaled    -2.18    1.21

This part is the SD of the penalized basis function coefficients:

Smooth terms:
                  Std. Dev.
sds(depth_scaled)      6.47

So, as this approaches zero, the smooth approaches a straight line with slope sdepth_scaled, here.

The approach uses mgcv::smooth2random.

Also, it's unfinished, but soon there may be an R^2 function, which can be broken down by covariate/smooth.