Is there a way to produce a pp_check density overlay plot with cmdstanr objects?
I have only found info on hacky ways to do this (e.g. turning your object into a brmsfit object and then working with that).
There is also this closed issue on the bayesplot repository.

Please help, I am losing faith in my ability to google…

You can’t use the pp_check function (like you can with brms), but you can use the underlying bayesplot functions that pp_check uses. In this case I think you’re referring to the plot made by bayesplot::ppc_dens_overlay(). Here’s an example of using it with cmdstanr:

// example.stan
data {
int<lower=0> N;
array[N] real y;
}
parameters {
real mu;
real<lower=0> sigma;
}
model {
y ~ normal(mu, sigma);
}
generated quantities {
array[N] real y_rep;
for (n in 1:N) y_rep[n] = normal_rng(mu, sigma);
}

N <- 100
y <- rnorm(N, mean = 0, sd = 1)
data_list <- list(N = N, y = y)
mod <- cmdstan_model("example.stan")
fit <- mod$sample(data = data_list)
y_rep <- fit$draws("y_rep", format = "matrix")
ppc_dens_overlay(y = data_list$y, yrep = y_rep[1:100, ])