Make something similar to plot.gam to quickly visualize smoothers

[seananderson]

The smooth predictions are already separate from the rest. You'd want what is currently called proj_smooth_i in the .cpp file. The only problem right now, is it gets reset with every smoother and not recorded. So, record the value of proj_smooth_i for each smoother.

Make a function that calls predict with re_form = NA and se_fit = TRUE with a data frame that has the smooth predictors going from min to max value and everything else set to zero and then plot the relevant proj_smooth_i predictions.

@joenomiddlename suggestion

[seananderson]

There's a working version in the sim2 branch now, but it includes uncertainty from all parameters including the intercept. It's likely worth doing a version that just plots proj_smooth_i for each smoother, which would better match mgcv. It would be relatively easy to do, but we'd have to record proj_smooth_i for each smoother separately and run ADREPORT on it if needed.

library(sdmTMB)
d <- subset(pcod, year >= 2000 & density > 0)
pcod_spde <- make_mesh(d, c("X", "Y"), cutoff = 30)
m <- sdmTMB(
  data = d,
  formula = log(density) ~ s(depth_scaled) + s(year, k = 5),
  mesh = pcod_spde
)
plot_smooth(m)

plot_smooth(m, ggplot = TRUE, level = 0.8)

plot_smooth(m, ggplot = TRUE, select = 2, level = 0.5)

^{Created on 2021-12-13 by the reprex package (v2.0.1)}

[seananderson]

@joenomiddlename points out these are called marginal and conditional smooth effect plots in brms: https://rdrr.io/cran/brms/man/conditional_smooths.brmsfit.html http://paul-buerkner.github.io/brms/reference/marginal_smooths.html

and it's nice to be able to do both.

[ericward-noaa]

I'm working on a couple things here, trying to get sdmTMB to play nice with 'visreg' and 'gratia'. First visreg: I have a couple tiny changes that I'll put in a pull request to visreg -- but in doing so, realized it might be easier if we changed est to fit and est_se to se.fit first. And perhaps the se_fit to se.fit in the arg list?

My changes to visreg are here: https://github.com/ericward-noaa/visreg/commit/0191c3a901b3d5b82da076daf9d15a08353b3ca0 Note: if we change names from predict.sdmTMB(), these would need to be changed before PR submitted

[ericward-noaa]

Example of this working in action:

pcod_2011$year = as.factor(pcod_2011$year)
pcod_spde <- make_mesh(pcod_2011, c("X", "Y"), cutoff = 25)
m <- sdmTMB(density ~ 0 + s(depth_scaled) + year,
            data = pcod_2011, mesh = pcod_spde, family = tweedie(link = "log"))
visreg::visreg(m,"depth_scaled")

In a couple early tests, I got warnings via this block: https://github.com/pbs-assess/sdmTMB/blob/5de93521ddfaabf42e2c51e8ecfed5d29b89d89a/R/predict.R#L290 Doesn't seem to be an issue now, but may revisit in the future

[ericward-noaa]

On the compatibility with gratia: that package only works with smooths constructed in mgcv. I think that's ok though, because we use the same mgcv::smoothCon() function. This comment above was initially thinking about doing it in TMB -- but I was thinking an equivalent approach would be to return the betas and do the predictions in R, e.g.

https://www.rdocumentation.org/packages/mgcv/versions/1.8-39/topics/smoothCon

Or maybe I'm missing something obvious about why this won't work?

[Lewis-Barnett-NOAA]

looks great and simple, thanks Eric.

[seananderson]

Very cool on visreg! That stop statement on mismatched time elements shouldn't be triggered if pop_pred is TRUE, which I assume it will be for these standard conditional effects plots that don't include the random fields.

Regarding smooths, we could recalculate them in R after the fact. I think the only downside would be that we couldn't get at the standard errors. It would also be fairly simple to keep track of them within TMB. The problem is just that they aren't recorded right now: https://github.com/pbs-assess/sdmTMB/blob/5de93521ddfaabf42e2c51e8ecfed5d29b89d89a/src/sdmTMB.cpp#L467-L470

We would have to record eta_smooth_i in a matrix and report it.

Regarding fit and se.fit over est and est_se. Hmm, yes, fit and se.fit would be more standard. The problem is just all the existing code out there based on est. I assume not as many people are using est_se yet. That said, sdmTMB still isn't even on CRAN, so perhaps we could change it and issue a warning for 6 months or so. We could depreciate se_fit to se.fit. Would be nice to consistently use underscores, but maybe consistency with other packages is more important.

For the data frame issue, I wonder if we could detect if predict is being called by whatever the visreg function is and return a list of fit and se.fit? https://stackoverflow.com/a/15621405 If modifying visreg, I imagine we'll need to be on CRAN first.

[seananderson]

I'm working on reporting and optionally adreporting an eta_smooth matrix right now...

[ericward-noaa]

I like the fact that we can detect which function's calling predict -- that'd minimize changes to visreg. I'll work on changing that later. Then there's no name changes -- and it should be easy to add packages alongside visreg if this comes up in the future.

I understand the eta_smooth_i thing now & uncertainty -- hadn't thought that through entirely. For gratia, we'll still want to use smoothCon to construct the smooths I think but I'll work on that later -- next weekish

[seananderson]

Is there a need to record the each smooth contribution within TMB then? We'd end up with a matrix of proj_smooth_i, one column per smoother and one row per row or predicted data, with optional standard errors.

[ericward-noaa]

Not sure -- I need to look into smoothCon more, and see if the objects can be created via the matrix & SEs

[ericward-noaa]

Ok, think it's fixed now with the new version being changes discussed above to predict() rather than on the visreg side. Working here:

[ericward-noaa]

Working off of the plot_smooth() function, I added a new function pred_smooth() on the gratia branch,

https://github.com/pbs-assess/sdmTMB/blob/d46aab3bc98a5fbdcef132a88fba2c01b09a6fb8/R/plot.R#L214

We can't pass sdmTMB objects to gratia::draw() directly, but pred_smooth() creates an object that gratia will recognize + plot. Usage is then

m <- sdmTMB(data = d, formula = log(density) ~ s(depth_scaled), mesh = pcod_spde)
ps <- pred_smooth(m, "s(depth_scaled)")
gratia::draw(ps)

1. Right now pred_smooth() is plotting est and est_se -- but if the eta_smooth_i and standard errors are being returned, we can switch to showing those instead
2. I made no changes to plot_smooth(), but it could be tweaked slightly to work with the df returned by pred_smooth()

[seananderson]

Awesome! We should deprecate and ultimately remove plot_smooth then, because it only works with simple 1D smooths. I’ll add those other prediction standard errors soon, which would match what is usually plotted with mgcv.

[seananderson]

I think plot_smooth() will be deprecated, because sdmTMB should be compatible with visreg in the delta (soon to be main) branch. Note that the residuals are randomized quantile residuals ones not the usual deviance ones.

library(sdmTMB)
pcod_2011$fyear <- as.factor(pcod_2011$year)
fit <- sdmTMB(
  density ~ s(depth_scaled, k = 5) + fyear,
  data = pcod_2011, mesh = pcod_mesh_2011,
  family = tweedie()
)
visreg::visreg(fit, xvar = "depth_scaled")

visreg::visreg(fit, xvar = "fyear")

visreg::visreg(fit, xvar = "depth_scaled", scale = "response")

fit_dg <- sdmTMB(
  density ~ s(depth_scaled, year, k = 8),
  data = pcod_2011, mesh = pcod_mesh_2011,
  spatial = "off",
  family = delta_gamma()
)
visreg_delta(fit_dg, xvar = "depth_scaled", model = 1, gg = TRUE)
#> Loading required namespace: ggplot2

visreg_delta(fit_dg, xvar = "depth_scaled", model = 2, scale = "response")

visreg2d_delta(fit_dg,
  xvar = "depth_scaled", yvar = "year",
  model = 2, scale = "response"
)

^{Created on 2022-05-03 by the reprex package (v2.0.1)}