Visual predictive checking with overlaid plots

In this notebook we demonstrate the differences of quantile dot plots, histograms, and KDE density plots when overlaying predictive draws from the model on the observation.

library(ggplot2)
library(ggdist)
library(bayesplot)

This is bayesplot version 1.15.0

- Online documentation and vignettes at mc-stan.org/bayesplot

- bayesplot theme set to bayesplot::theme_default()

   * Does _not_ affect other ggplot2 plots

   * See ?bayesplot_theme_set for details on theme setting

library(cmdstanr)

This is cmdstanr version 0.9.0

- CmdStanR documentation and vignettes: mc-stan.org/cmdstanr

- CmdStan path: /home/teemu/.cmdstan/cmdstan-2.38.0

- CmdStan version: 2.38.0

source("code/R/helpers.R")
source("code/R/helpers_qdotplot.R")
source("code/R/helpers_kde.R")
source("code/R/helpers_histogram.R")

theme_set(ppc_paper_theme(28) +
  theme(
    axis.title = element_blank(),
    legend.position = "none",
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    axis.line.y = element_blank()
  ))

bayesplot_theme_set(ppc_paper_theme(28) +
  theme(
    axis.title = element_blank(),
    legend.position = "none",
    axis.ticks = element_blank(),
    axis.text = element_blank(),
    axis.line.y = element_blank()
  ))

color_scheme_set(unname(paper_colors[1:6]))
set.seed(9865875)

We start by creating a normally distributed sample.

n <- 1e3
nq <- 100
bw_dots <- sqrt(1 / nq)
y <- rnorm(n, sd = .4)

We fit the following Stan model to the sample and extract predictive draws.

model <- cmdstan_model(stan_file = "code/stan/02_05_normal.stan")

model$print()

data {
  int<lower=1> N;  // total number of observations
  vector[N] Y;  // response variable
}
parameters {
  real Intercept;  // temporary intercept for centered predictors
  real<lower=0> sigma;  // dispersion parameter
}
transformed parameters {
  real lprior = 0;  // prior contributions to the log posterior
  lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
  lprior += student_t_lpdf(sigma | 3, 0, 2.5)
    - 1 * student_t_lccdf(0 | 3, 0, 2.5);
}
model {
  // likelihood including constants
  // initialize linear predictor term
  vector[N] mu = rep_vector(0.0, N);
  mu += Intercept;
  target += normal_lpdf(Y | mu, sigma);
  // priors including constants
  target += lprior;
}
generated quantities {
  // actual population-level intercept
  real b_Intercept = Intercept;
  vector[N] yrep;
  yrep = to_vector(normal_rng(rep_vector(Intercept, N), sigma));
}

fit <- model$sample(data = list(N = n, Y = y), refresh = 0, thin = 50)

Running MCMC with 4 sequential chains...

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpJSVvrv/model-625484a274d61.stan', line 20, column 2 to column 39)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpJSVvrv/model-625484a274d61.stan', line 20, column 2 to column 39)

Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 1 

Chain 1 finished in 0.1 seconds.

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpJSVvrv/model-625484a274d61.stan', line 20, column 2 to column 39)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:

Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/tmp/RtmpJSVvrv/model-625484a274d61.stan', line 20, column 2 to column 39)

Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,

Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.

Chain 2 

Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.8 seconds.

yrep <- fit$draws(variables = "yrep", format = "matrix")
## For histogram binning:
yrep_h <- yrep
yrep_h[yrep_h < min(y)] <- min(y)
yrep_h[yrep_h > max(y)] <- max(y)

1 PPC with quantile dot plots

We plot the observation as quantile dot plot with 100 quantiles. Then, we overlay just the top dot of each stack to visualize the variation in the quantile dot plots of the predictive draws. As the dots from the same draw are not linked in the visualization, we can’t assess the variation in the overall shape between draws.

qdots <- ggplot(
  bin_dots(
    x = quantile(y, stats::ppoints(nq, a = 1 / 2), type = 5),
    y = 0,
    binwidth = bw_dots
)) +
  ggforce::geom_circle(
    aes(x0 = x, y0 = y, r = .5 * bw_dots),
    alpha = .2,
    fill = paper_colors["light_highlight"],
    data = apply(
      yrep[1:50,],
      1,
      \(row) bin_dots(
        x = quantile(row, stats::ppoints(nq, a = 1 / 2), type = 5),
        y = 0,
        binwidth = bw_dots
      ) |>
        dplyr::group_by(bin) |>
        dplyr::summarise(x = min(x), y = max(y)) |>
        dplyr::select(x, y, bin)
    ) |> dplyr::bind_rows(),
    linewidth = 0,
    n = 10
  ) +
    ggforce::geom_circle(
    aes(x0 = x, y0 = y, r = 0.5 * bw_dots),
    color = "black",
    linewidth = .5
  ) +
  coord_equal()

qdots

2 PPC with histogram

We compute the bin width with the observation and use this and the bin breaks also for visualizing the posterior predictive draws.

bw_hist <- 2 * IQR(y) / length(y)^(1 / 3)

Overlaying the predictive mean and credible interval of each bin only allows us to compare the predictions to the observation per individual bin. Like with quantile dot plots, we can’t see the variation in hte overall shape of the distribution between the predictive draws.

p_hist <- ggplot(data.frame(y = y)) +
  geom_histogram(
    aes(
      x = y,
      y = after_stat(density)
    ),
    fill = "white", #paper_colors["light_highlight"],
    color = paper_colors["dark_highlight"],#paper_colors["mid"],
    linewidth = .6,
    binwidth = bw_hist,
    center = mean(range(y))
  )

p_hist <- p_hist +
  geom_pointinterval(
    data = yrep_h |> apply(1, \(row) hist(
      row,
      breaks = c(layer_data(p_hist)$xmin[1], layer_data(p_hist)$xmax),
      plot = F
    )[c("mids", "density")]) |> dplyr::bind_rows() |> dplyr::group_by(mids) |> dplyr::summarise(
      y = mean(density),
      ymin = quantile(density, .05),
      ymax = quantile(density, .95)
    ),
    aes(
      x = mids,
      y = y,
      ymin = ymin,
      ymax = ymax
    ),
    color = paper_colors["mid"],#paper_colors["dark_highlight"],
    interval_color = paper_colors["mid"], #paper_colors["dark"],
    linewidth = 5.0,
    size = 5.0
  ) +
  coord_equal(
    xlim = layer_scales(qdots)$x$range$range,
    ylim = layer_scales(qdots)$y$range$range
  )

p_hist

3 PPC with overlaid KDE density plots

bayesplot already has overlaid KDE densities implemented. This is a common visual PPC plot and shows an overview of the observation and draws from the predictive distribution.

ppc_dens_overlay(y, yrep, alpha = 1) +
  coord_equal(
    xlim = layer_scales(qdots)$x$range$range,
    ylim = layer_scales(qdots)$y$range$range
  )

Coordinate system already present.
ℹ Adding new coordinate system, which will replace the existing one.