Add option to expose samples from derivatives about gratia HOT 6 CLOSED

If I follow correctly, derivatives() with a simultaneous intervals doesn't actually generate samples of the smooth at locations x and shifted values x + eps (say) and from each sample compute the derivative. It just simulates deviations from expected value (derivative) based on the uncertainty in the parameters. If that's sufficient for your use (but it doesn't sound like it) I could see how to expose those.

If you really want samples of the derivatives, would response_derivatives() work (suitably hacked to return the samples)? It uses fitted_samples() to generate draws from the posterior of the expected values of the response (posterior_samples() draws new response values from the posterior distribution) at x and x + eps (or other shift pattern for the different types of derivative) and computes the derivative as usual. You'd just need to hack it (temporarily) to stop it doing the summarise step where we lose the samples and just return the median and some upper (lower) tail quantiles. This is the part to delete/change.

If that is what you want, I will see about adding a derivatives_samples() function that does most of what response_derivatives() does but without the summary step.

If you want these for a specific smooth (which is what derivatives() does but without samples), you can use the ... argument to response_derivatives() to pass arguments to predict.gam(), which if you use terms = "s(x)" you'll get expected values just for the smooth with label "s(x)", excluding the intercept with recent versions of {mgcv}.

I can't recall if response derivatives() is in a released version of {gratia}, but I suspect not, so be sure to install the development version if you want to fiddle with all this.

from gratia.

Assuming I understood what you wanted, will this work:

library("ggplot2")
load_mgcv()
df <- data_sim("eg1", dist = "negbin", scale = 0.25, seed = 42)

# fit the GAM (note: for execution time reasons using bam())
m <- bam(y ~ s(x0) + s(x1) + s(x2) + s(x3),
    data = df, family = nb(), method = "fREML", discrete = TRUE)

# data slice through data along x2 - all other covariates will be set to
# typical values (value closest to median)
ds <- data_slice(m, x2 = evenly(x2, n = 200))

# samples from posterior of derivatives
fd_samp <- derivative_samples(m, data = ds, type = "central",
    focal = "x2", eps = 0.01, seed = 21)

# plot the first 20 posterior draws
fd_samp |>
    filter(draw <= 20) |>
    ggplot(aes(x = x2, y = derivative, group = draw)) +
    geom_line(alpha = 0.5)

where fd_samp looks like this:

> fd_samp
# A tibble: 2,000,000 × 8
    .row focal  draw derivative      x2    x0    x1    x3
   <int> <chr> <int>      <dbl>   <dbl> <dbl> <dbl> <dbl>
 1     1 x2        1       21.0 0.00359 0.488 0.477 0.488
 2     2 x2        1       21.9 0.00860 0.488 0.477 0.488
 3     3 x2        1       22.8 0.0136  0.488 0.477 0.488
 4     4 x2        1       23.9 0.0186  0.488 0.477 0.488
 5     5 x2        1       25.1 0.0236  0.488 0.477 0.488
 6     6 x2        1       26.6 0.0286  0.488 0.477 0.488
 7     7 x2        1       28.3 0.0336  0.488 0.477 0.488
 8     8 x2        1       30.4 0.0386  0.488 0.477 0.488
 9     9 x2        1       32.9 0.0436  0.488 0.477 0.488
10    10 x2        1       35.8 0.0486  0.488 0.477 0.488
# ℹ 1,999,990 more rows
# ℹ Use `print(n = ...)` to see more rows

I haven't quite settled on the correct names just yet, but .row indexes the row of the object passed to data in the derivative_samples() call, draw groups posterior draws together, and the other the other variables are the estimated derivative and the data values at which the derivative was evaluated.

If you just want the derivative of the smooth of $f(\mathtt{x2})$, s(x2), use

fd_samp <- derivative_samples(m, data = ds, type = "central", 
    focal = "x2", eps = 0.01, seed = 21, terms = "s(x2)")

Using terms and exclude you can select any combination of terms you want the derivatives for. I'm very tempted to add an argument include which would be passed to the terms argument of predict.[g|b]am() rather than just allow terms to be passed on. I think using include or exclude is more consistent. but for now, that's what I have added to the development version of {gratia}.

from gratia.

@gavinsimpson this is exactly what I was hoping for! Thank you! Have tested the development branch and works as expected. My use case is essentially a y ~ s(t, by=r) + s(r, bs='re') so limited need for the different term bits - I want to see the variation over t to extract when the maximum rate is.

Spotted the docs may need to be updated as the values are for example the same as the interval approach:

from gratia.

Would it be best to do row rather than .row to keep consistency with posterior_samples?

from gratia.

Thanks; glad it is working OK for you. I'm fixing the documentation periodically now that I have fixed on a naming convention...

Would it be best to do row rather than .row to keep consistency with posterior_samples?

I heading in the other direction; converging on using the prefix . on variables generated in {gratia} functions where we mix user-data/variables and {gratia}-generated ones.

from gratia.

Sounds good! I suppose . matches how brms and other tidybayes does this.

from gratia.