This notebook reproduces and expands on the examples of predictive checks for ordinal data, shown in Section 6 of the article.
library("bayesplot")
library("cmdstanr")
library("patchwork")
library("ggplot2")
library("khroma")
library("dplyr")
# Source for the modified reliability plot
source("code/R/helpers.R")
source("code/R/pava_plots.R")
theme_set(ppc_paper_theme())
bayesplot_theme_set(ppc_paper_theme())
SEED <- 236543
set.seed(SEED)
SAVE_FITS <- TRUEIn assessing the calibration, we use the ordinal nature of the data and use the cumulative posterior predictive mass function.
We generate data by drawing 1500 observations from a mixture of 5 Gaussians with means 1, 2, 3, 4, and 5, and standard deviation \(0.5\).
\[\begin{align} x_n &\sim \mathcal N\!\!\left(k, 0.5^2\right), &\text{for } n \in\{1,\dots,N\}\\ k &\sim \text{Categorical}(\theta_k),&\\ \theta_k &= \frac 1 K, &\text{for } k \in \{1, \dots, K\}. \end{align}\]
data.frame(
x = rep(seq(min(standata_gmm$x), max((standata_gmm$x)), .01), K),
d = 1:K |> sapply(\(c) dnorm(seq(min(standata_gmm$x), max((standata_gmm$x)), .01), c, sigma)) |> c(),
c = rep(1:K, each = length(seq(min(standata_gmm$x), max((standata_gmm$x)), .01)))
) |>
group_by(x) |>
mutate(percentage = d / sum(d)) |>
ggplot() +
geom_line(
aes(
x = x,
y = percentage,
color = as.factor(c),
group = c
),
linewidth = 2,
key_glyph = "rect"
) +
scale_color_manual(
aesthetics = c("color"),
values = unname(paper_colors[1:6])
) +
labs(color = "Value", x = expression("x"), y = "Proportion") +
coord_equal(ratio = 2, xlim = range(standata_gmm$x), ylim = c(0, 1), expand = FALSE) +
theme(legend.position = "top")
We fit two models to the data, both structured to first normalize the data and then fit a Gaussian mixture model (GMM) with K = 5` mixture components. The first model replicates a simple implementation error where standard deviation of the mixture components was incorrectly scaled.
data {
int<lower=2> K; // Number of classes
int<lower=0> N; // Total number of observations
array[N] int<lower=1, upper=K> y; // Target classes
vector[N] x; // Observed predictor values
real<lower=0> sigma; // User supplied standard deviation
int correct_sigma; // How to handle sigma, see below.
}
transformed data{
vector[N] x_st;
real Sigma;
// Standardize data
x_st = (x - mean(x)) / sd(x);
// Maybe remember to scale sigma accordingly
if (correct_sigma == 1) {
Sigma = sigma / sd(x);
} else {
Sigma = sigma;
}
}
parameters {
// Inferred means.
ordered[K] c;
simplex[K] p_c;
}
model {
// Prior
c ~ normal(0,1);
p_c ~ dirichlet(rep_vector(1,K));
// Likelihood
for (n in 1:N) {
target += normal_lpdf(x_st[n] | c[y[n]], Sigma);
}
}
generated quantities {
// Posterior predictive sample
vector[N] yrep;
// For each observation, posterior predictive mass of classes
array[N] vector[K] ppm;
for (n in 1:N) {
for (k in 1:K) {
ppm[n, k] = normal_lpdf(x_st[n] | c[k], Sigma);
}
ppm[n, ] = softmax(ppm[n, ]);
yrep[n] = categorical_rng(ppm[n, ]);
}
}First, we use a bar graph to visualize the relative frequencies of the groups in the observation and the predictive draws. We don’t see much difference between the models, and all of the observations seem to fall within the credible intervals of the predictions of both models.
p2 <- ppc_bars(
y = as.numeric(c),
yrep = fit_2$draws(variables = "yrep", format = "matrix")
) +
ggtitle("Model 2") +
theme(
axis.title.y = element_blank(),
axis.text.y = element_blank()
)
p1 <- ppc_bars(
y = as.numeric(c),
yrep = fit_1$draws(variables = "yrep", format = "matrix")
) +
ggtitle("Model 1")
(p1 + p2) + plot_layout(guides = "collect") & theme(legend.position = "bottom")
When using the PAV-adjusted calibration plots, we see that the first model, with the initial implementation error is under confident (S-shaped calibration curves) in its predictions, and assigns probabilities too close to 0.5 to many cases that the other model is better able to separate towards the observed outcomes.
Loading required package: reliabilitydiagLoading required package: ggdist
