PPC Visualizations for Categorical Data

This notebook reproduces and expands on the examples of predictive checks for categorical data, shown in Section 5 of the article.

library(ggplot2)
library(bayesplot)

This is bayesplot version 1.11.1

- 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(caret)

Loading required package: lattice

source("code/R/helpers.R")
source("code/R/helpers_kde.R")
source("code/R/pava_plots.R")

SEED <- 2451
set.seed(SEED)

theme_set(
  ppc_paper_theme() +
  theme(
    panel.grid = theme_minimal()$panel.grid
  )
)

bayesplot_theme_set(ppc_paper_theme())

We create 500 observations where each observation has independent uniform weights for belonging to each class and then normalize these to obtain class probabilities. We add a 30% chance of switching the probabilities of classes B and C in the process, to induce confusion between the observations and the predictive probabilities.

We then draw predictive samples using the normalized class weights without confusion.

n_obs <- 500
probs <-
  matrix(runif(3 * n_obs),
         ncol = 3,
         dimnames = list(NULL, c("pA", "pB", "pC")))
probs <- data.frame(probs / rowSums(probs))

with(probs,
  y <<- as.factor(sapply(1:n_obs,
    function(idx) {
      
      sample(c("A", "B", "C"),
             size = 1,
             # We confuse B and C.
             prob = c(pA[idx], pB[idx], pC[idx])[c(1, 1 + sample.int(2,prob = c(.7,.3)))]
             )
  }))
)
yrep <- t(replicate(100,
                   with(probs,
  as.factor(sapply(1:n_obs,
    function(idx) {
      sample(c("A", "B", "C"),
             size = 1,
             prob = c(pA[idx], pB[idx], pC[idx])
             )
  }))
)))

In the bar graph, we can’t tell if there is anything wrong with the predictions.

ppc_bars(
  as.integer(y),
  matrix(as.integer(as.factor(yrep)), ncol = ncol(yrep)),
  freq = FALSE) +
  scale_x_continuous(breaks = 1:3, labels = c("A", "B", "C")) +
  theme(legend.position = "none")

Scale for x is already present.
Adding another scale for x, which will replace the existing scale.

We then proceed to make one-vs-others binned calibration and PAV-calibration plots for the predictions.

For A, everything looks fine.

ggplot(calibration(
  y ~ prob,
  data = data.frame(prob = probs$pA, y = y),
  cuts = 10,
  class = "A"
)$data |> dplyr::filter(Count > 0), aes(x = midpoint / 100, y = Percent / 100)) +
  geom_abline(slope = 1, intercept = 0, col = "black", lty = 2, alpha = .3) +
  geom_point(color = paper_colors["dark_highlight"]) +
  geom_errorbar(aes(ymin = Lower / 100, ymax = Upper / 100), width = .02, color = paper_colors["dark"]) +
  labs(
    x = "Predicted probability",
    y = "Observed rate",
    title = "A vs. others"
  ) +
  coord_equal(xlim = c(0,1.02), ylim = c(0,1.02), expand = FALSE) +
  ppc_paper_theme() +
  theme(
    panel.grid = theme_minimal()$panel.grid
  )

with(probs, {
  ppc_calibration_pava(
    y = as.numeric(y == "A"),
    p = pA,
    fill_alpha = .4,
    cep_line_color = paper_colors["orange"]) +
        labs(title = "A vs. others")
}
)

Loading required package: reliabilitydiag

Loading required package: ggdist

Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

For B and C below, we clearly see the confusion in the predictions.

ggplot(calibration(
  y ~ prob,
  data = data.frame(prob = probs$pB, y = y),
  cuts = 10,
  class = "B"
)$data |> dplyr::filter(Count > 0), aes(x = midpoint / 100, y = Percent / 100)) +
  geom_abline(slope = 1, intercept = 0, col = "black", lty = 2, alpha = .3) +
  geom_point(color = paper_colors["dark_highlight"]) +
  geom_errorbar(aes(ymin = Lower / 100, ymax = Upper / 100), width = .02, color = paper_colors["dark"]) +
  labs(
    x = "Predicted probability",
    y = "Observed rate",
    title = "B vs. others"
  ) +
  coord_equal(xlim = c(0,1.02), ylim = c(0,1.02), expand = FALSE) +
  ppc_paper_theme() +
  theme(
    panel.grid = theme_minimal()$panel.grid
  )

with(probs, {
  ppc_calibration_pava(
    y = as.numeric(y == "B"),
    p = pB,
    fill_alpha = .4,
    cep_line_color = paper_colors["orange"]) +
    labs(
      title = "B vs. others"
    )
})

ggplot(calibration(
  y ~ prob,
  data = data.frame(prob = probs$pC, y = y),
  cuts = 10,
  class = "C"
)$data |> dplyr::filter(Count > 0), aes(x = midpoint / 100, y = Percent / 100)) +
  geom_abline(slope = 1, intercept = 0, col = "black", lty = 2, alpha = .3) +
  geom_point(color = paper_colors["dark_highlight"]) +
  geom_errorbar(aes(ymin = Lower / 100, ymax = Upper / 100), width = .02, color = paper_colors["dark"]) +
  labs(
    x = "Predicted probability",
    y = "Observed rate",
    title = "C vs. others"
  ) +
  coord_equal(xlim = c(0,1.02), ylim = c(0,1.02), expand = FALSE) +
  ppc_paper_theme() +
  theme(
    panel.grid = theme_minimal()$panel.grid
  )

with(probs, {
  ppc_calibration_pava(
    y = as.numeric(y == "C"),
    p = pC,
    fill_alpha = .4,
    cep_line_color = paper_colors["orange"]) +
    labs(
      title = "C vs. others"
    )
})