I am fitting GAMs in R using the mgcv R package. I then want to perform posterior predictive checks to compare observed and simulated outcomes. A toy example using the performance R package:

library(mgcv)
library(performance)

df <- mtcars |>
  tibble::as_tibble() |>
  dplyr::mutate(carb = as.integer(carb))

model <- mgcv::gam(
  carb~s(wt, k=5),
  data = df,
  family = mgcv::ocat(R=8)
)

plot(check_model(model))

That does not work:

Error in d[, !(names(data) %in% all.varying), drop = FALSE] : 
  incorrect number of dimensions

I am trying to get something similar to this:

The only issue seems to be that when the performance calls simulate, what it gets is not supported. I can obtain simulated outcomes using mgcViz:

library(ggplot2)
sims <- mgcViz::simulate.gam(model, newdata=df, nsim=10) |> tibble::as_tibble() |>
    tidyr::pivot_longer(dplyr::everything())
ggplot(sims, aes(x = value)) +
    geom_point(aes(y = after_stat(count)), stat = "count")

The problem is that I only obtain one point, not the jitter of points as in the image above. How can I achieve that?