Plot periods of significant change from a GAM on the response scale using Gratia

[MasonWirtz]

I am curious whether there is a straightforward way to plot periods of significant change using generalised additive models (GAMs), similar to what Gavin proposed here (https://fromthebottomoftheheap.net/2014/05/15/identifying-periods-of-change-with-gams/). Using the code below, one can do this relatively easily for the derivatives, but I am unsure of how to plot periods of significant change on the response scale. One proposed method here (https://stats.stackexchange.com/questions/600783/test-statistics-for-significant-change-points-from-first-derivative-of-gam-smoot) demonstrates a work-around, but I wonder whether there is an easier way now after several years. Thanks for any suggestions!

df <- data_sim("eg1", seed = 42)

mod <- bam(y ~ s(x2), data = df)

df <- derivatives(mod, select = "s(x2)")

draw(df, select = "s(x2)", add_change = TRUE, change_col = "black",
     increase_col = "blue", decrease_col = "red")

[gavinsimpson]

With the big caveat that the derivatives here are on the scale of the linear predictor not the response (although in this specific example and the one in the blog post these are the same thing — if you want derivatives on the response scale, see gratia::response_derivatives() or marginaleffects::slopes()) here's how to do this.

library("ggplot2")
library("dplyr")
library("mgcv")
library("gratia")

# example data
df <- data_sim("eg1", seed = 42)

# model
m <- bam(y ~ s(x2), data = df)

# best to make sure everything uses the same data, so new data evenly over x2
ds <- m |>
  data_slice(
    x2 = evenly(x2)
  ) |>
  mutate(
    .row = row_number() # add a row number just so you can assure yourself everything in order
  )

# derivatives on linear predictor scale
d <- derivatives(m, select = "s(x2)", data = ds) |>
  mutate(
    .row = row_number() # row number for joining to fitted values
  ) |>
  add_sizer(type = "sizer") # add sizer indicators

# compute fitted values
fv <- fitted_values(m, data = ds) |>
  left_join( # join with the derivatives
    d |> select(.decrease, .increase, .row), # only want a few columns in `d`
    by = join_by(.row == .row) # match on row number
  ) |>
  mutate( # update info from derivatives to
    .decrease = case_when( # replace .decrease derivative with fitted value
      !is.na(.decrease) ~ .fitted,
      .default = NA_real_
    ),
    .increase = case_when( # and again but for increase
      !is.na(.increase) ~ .fitted,
      .default = NA_real_
    )
  )

fv |> # plot
  ggplot(
    aes(x = x2, y = .fitted)
  ) +
  geom_ribbon(
    aes(ymin = .lower_ci, ymax = .upper_ci), alpha = 0.2
  ) +
  geom_line() +
  geom_line(
    aes(y = .decrease), colour = "blue", linewidth = 2
  ) +
  geom_line(
    aes(y = .increase), colour = "red", linewidth = 2
  )

This produces

[MasonWirtz]

This is fantastic, very many thanks, Gavin! Exactly what I was looking for!