Visualizing marginal effects in sdmTMB models with ggeffects
Julia Indivero, Sean Anderson, Lewis Barnett, Philina English, Eric Ward
If the code in this vignette has not been evaluated, a rendered version is available on the documentation site under ‘Articles’.
The package ggeffects can be used to plot marginal effects of given predictor variables in sdmTMB models.
A advantage to this approach over visreg, is that ggeffects calculates marginal effects with the effects package using the parameter covariance matrix. This is nearly instant compared to visreg, which has to calculate conditional effects by calculating predictions with TMB.
A disadvantage to using ggeffects is that it will only work for
regular linear effects in the main model formula. I.e., it will not work
with smoothers (internally these are random effects) or breakpoint
Another important distinction is that
ggeffects::ggeffect() is plotting marginal
effects. This means the effects are “marginalized” or “averaged” over
the other fixed effects.
visreg::visreg() is plotting
conditional effects. This means they are conditional on
the other predictors being set to certain values.
ggeffects::ggpredict() also does conditional effects, but
this has not yet been set up in sdmTMB using the CRAN version of
Example with Pacific cod presence
To start, we will use the Pacific cod example data. We will fit a model of fish presence/absence with covariates of depth and a fixed effect of year using a Tweedie distribution.
pcod$fyear <- as.factor(pcod$year) mesh <- make_mesh(pcod, c("X", "Y"), cutoff = 20) fit <- sdmTMB(present ~ poly(depth, 2) + fyear, data = pcod, mesh = mesh, spatial = "on", family = binomial() )
We can then use
ggeffects::ggeffect() to see the effect
of depth on the probability of Pacific cod being present. We can control
what range and interval of depths are predicted within the function
We can also plot the effects of each year.
We can add in data points
plot(g, add.data = TRUE)
We can also use
ggeffect to plot multiple variables by
listing them in
terms = c(), with the first term listed
indicating the variable to be plotted on the x-axis, and the remaining
listed terms (up to four total) indicating the groups. Adding
facet = TRUE will show each year as a separate plot,
instead of overlain on one plot.
facet = TRUE will show each year as a separate
plot, instead of overlain on one plot.
plot(dat, facet = TRUE)
We can also use make our own ggplot plot by calling the ggeffects
dat as the data frame.
Plotting using with a continuous response (here
rather than presence-only is similar. For instance:
We can fit a model with an interaction of two continuous variables:
For plotting two continuous variables,
make the non-target (2nd) variable discrete by selecting different
plot(g5, facet = TRUE)
To specify the levels rather than letting
choose them, use brackets with the selected values within the term list,