Issue with plotting using visreg_delta each model

[LSupatcha]

Hi! I’m currently getting an error when I try to plot for each model using visreg_delta. I’m using sdmTMB_0.6.0.

Here’s my code:
df$fyear_season <- as.factor(df$year_season)
fit <- sdmTMB(
list(density ~ fyear_season+factor(Fleet) + poly(cov1, 2) + poly(cov2, 2) + poly(cov3, 2),
density ~ fyear_season+factor(Fleet) + poly(cov2, 2) + poly(cov3, 2)),
data = df,
mesh = mesh,
time = "year_season",
family = delta_gamma(link1 = "logit", link2 = "log"),
spatial = list("on", "on"),
spatiotemporal = list("ar1", "ar1"),
control = sdmTMBcontrol(eval.max = 10000, iter.max = 10000),
anisotropy = TRUE,
silent = FALSE
)

Using the following commands, I can get the plots for all three covariates for model 1:
visreg_delta(fit, xvar = "cov1", model = 1, gg = TRUE)
visreg_delta(fit, xvar = "cov2", model = 1, gg = TRUE)
visreg_delta(fit, xvar = "cov3", model = 1, gg = TRUE)
But when I tried to plot the covariates in model 2 using this command…
visreg_delta(fit, xvar = "cov2", model = 2, gg = TRUE)
visreg_delta(fit, xvar = "cov3", model = 2, gg = TRUE)

I get this error: “Error in eval(predvars, data, env) : object “cov1” not found” even if I already removed cov1 in fitting model 2.

May I kindly ask for your suggestions on how to resolve this issue?

[seananderson]

Thanks, I can confirm that issue. It's related to the predict code always iterating over building the model matrices for both linear predictors but visreg prepares a data frame here based only on the 2nd linear predictor. I.e., it's going to happen if the 2 linear predictors have different formulas.

Here's an example of making a conditional effects plot yourself:

library(sdmTMB)
library(ggplot2)
library(dplyr)

df <- pcod
df$fyear <- as.factor(df$year)
df$cov1 <- df$depth
df$cov2 <- df$year

fit <- sdmTMB(
  list(density ~ poly(cov1, 2) + poly(cov2, 2),
    density ~ poly(cov2, 2)),
  data = df,
  time = "year",
  family = delta_gamma(),
  spatiotemporal = "off",
  spatial = "off"
)

nd <- expand.grid(cov1 = mean(df$cov1), year = df$year[1],
  cov2 = seq(min(df$cov2), max(df$cov2), length.out = 100))
head(nd)
#>       cov1 year     cov2
#> 1 190.4401 2003 2003.000
#> 2 190.4401 2003 2003.141
#> 3 190.4401 2003 2003.283
#> 4 190.4401 2003 2003.424
#> 5 190.4401 2003 2003.566
#> 6 190.4401 2003 2003.707

p <- predict(fit, newdata = nd, se_fit = TRUE, re_form = NA, model = 2)
p <- mutate(p, upr = est + 1.96 * est_se, lwr = est - 1.96 * est_se)

ggplot(p, aes(cov2, exp(est), ymin = exp(lwr), ymax = exp(upr))) +
  geom_ribbon(fill = "grey60") + geom_line()

^{Created on 2024-06-20 with reprex v2.1.0}

In terms of visreg, it's a major pain to accommodate within sdmTMB and I think I'd rather put the time into getting delta models working properly with ggpredict and ggeffect in the ggeffects package and deprecate visreg support.

ggeffects:
https://pbs-assess.github.io/sdmTMB/articles/web_only/ggeffects.html

There aren't examples of ggpredict in their yet, but you can see examples here:
https://github.com/strengejacke/ggeffects/blob/master/tests/testthat/test-sdmTMB.R

[LSupatcha]

Thank you so much for your response. I have tried the one that you suggested. It works with a few covariates, but when there are many covariates, expand.grid becomes very large and requires high computation, causing it to stop running on my PC.

By the way, I decided to run the model based on separate components to get the conditional effect using the visreg function. but I'm wondering why, after trying it, the results between the built-in delta_gamma model and the combined model using binomial and gamma models are different, even with the same assumptions.

[seananderson]

I have tried the one that you suggested. It works with a few covariates, but when there are many covariates, expand.grid becomes very large and requires high computation, causing it to stop running on my PC.

For these types of plots, you'd typically be predicting on a data frame with a sequence of values for the variable of interest with all other variables held at a constant value. I.e., expand.grid() would only be needed if you were making a condition effects plot for interacting variables, e.g., a bivariate plot. So, regardless of how many covariates you have, I wouldn't expect the data frame to become large. expand.grid() wasn't really needed above. Internally, this is all that visreg or ggeffects::ggpredict() are doing.

By the way, I decided to run the model based on separate components to get the conditional effect using the visreg function. but I'm wondering why, after trying it, the results between the built-in delta_gamma model and the combined model using binomial and gamma models are different, even with the same assumptions.

In a default setup (e.g., no mapping parameters to be shared across linear predictors), and for families that can be fit independently (non "Poisson-link" families), they should match exactly. Is it possible it's because the mesh vertices aren't the same? The reasons for fitting them as a 'delta'/hurdle family are (1) convenience, (2) the ability to share parameters if desired, (3) the ability to calculate derived quantities (with standard errors) that depend on both linear predictors (e.g. the get_index() function). Here's an example of fitting with a delta_gamma() family vs. fitting binomial and gamma separately:

library(sdmTMB)
mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 10)
fit <- sdmTMB(
  density ~ depth_scaled,
  data = pcod, mesh = mesh,
  family = delta_gamma()
)

# Fit first linear predictor:
fit1 <- sdmTMB(
  present ~ depth_scaled,
  data = pcod, mesh = mesh,
  family = binomial()
)

# Fit second linear predictor:
pcod_pos <- subset(pcod, present == 1)
# force mesh to use same vertices:
mesh_pos <- make_mesh(pcod_pos, c("X", "Y"), mesh = mesh$mesh)
fit2 <- sdmTMB(
  density ~ depth_scaled,
  data = pcod_pos, mesh = mesh_pos,
  family = Gamma(link = "log")
)

tidy(fit, model = 1)
#> # A tibble: 2 × 5
#>   term         estimate std.error conf.low conf.high
#>   <chr>           <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)    -0.313     0.416    -1.13     0.503
#> 2 depth_scaled   -0.974     0.170    -1.31    -0.641

tidy(fit1)
#> # A tibble: 2 × 5
#>   term         estimate std.error conf.low conf.high
#>   <chr>           <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)    -0.313     0.416    -1.13     0.503
#> 2 depth_scaled   -0.974     0.170    -1.31    -0.641

tidy(fit, "ran_pars", model = 1)
#> # A tibble: 2 × 5
#>   term    estimate std.error conf.low conf.high
#>   <chr>      <dbl>     <dbl>    <dbl>     <dbl>
#> 1 range      30.5      6.06     20.6      45.0 
#> 2 sigma_O     2.22     0.222     1.82      2.70

tidy(fit1, "ran_pars")
#> # A tibble: 2 × 5
#>   term    estimate std.error conf.low conf.high
#>   <chr>      <dbl>     <dbl>    <dbl>     <dbl>
#> 1 range      30.5      6.06     20.6      45.0 
#> 2 sigma_O     2.22     0.222     1.82      2.70

tidy(fit, model = 2)
#> # A tibble: 2 × 5
#>   term         estimate std.error conf.low conf.high
#>   <chr>           <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)     4.14     0.0989    3.94      4.33 
#> 2 depth_scaled   -0.333    0.114    -0.556    -0.110

tidy(fit2)
#> # A tibble: 2 × 5
#>   term         estimate std.error conf.low conf.high
#>   <chr>           <dbl>     <dbl>    <dbl>     <dbl>
#> 1 (Intercept)     4.14     0.0989    3.94      4.33 
#> 2 depth_scaled   -0.333    0.114    -0.556    -0.110

tidy(fit, "ran_pars", model = 2)
#> # A tibble: 3 × 5
#>   term    estimate std.error conf.low conf.high
#>   <chr>      <dbl>     <dbl>    <dbl>     <dbl>
#> 1 range      5.17     7.44      0.308    86.8  
#> 2 phi        0.744    0.0303    0.687     0.805
#> 3 sigma_O    2.27     2.86      0.193    26.7

tidy(fit2, "ran_pars")
#> # A tibble: 3 × 5
#>   term    estimate std.error conf.low conf.high
#>   <chr>      <dbl>     <dbl>    <dbl>     <dbl>
#> 1 range      5.17     7.44      0.308    86.8  
#> 2 phi        0.744    0.0303    0.687     0.805
#> 3 sigma_O    2.27     2.86      0.193    26.7

^{Created on 2024-07-12 with reprex v2.1.0}