Primary Supplement for ‘When the Mean is Meaningless’

1 Introduction

This file is the online supplement for the paper “When the Mean is Misleading: A Guide to Ordered Regression for Meaningfully Modeling Ordinal Outcomes” (Jonathan R. Brauer & Jacob C. Day). Our primary aim with this paper is to show how and why to use ordered regression techniques like cumulative probit regression when analyzing ordinal outcomes (and perhaps count and continuous outcome too!). We illustrate via simulations and then by applying these methods in an analysis of data from Pickett, Graham, and Cullen’s 2022 study reported in Criminology in “The American racial divide in fear of the police” and shared as part of their supplemental materials on OSF.

This primary supplement contains all R code and output reported in the manuscript as well as other supplementary materials. In addition to this primary supplement, we have also provided a standalone tutorial in a separate Quarto file (Appendix A) that provides a step-by-step guide to our progression from frequentist and Bayesian versions of linear, to ordinal, to multilevel ordinal models.

1.1 Why this paper?

We picked Pickett, Graham, & Cullen’s paper for several reasons. First, we wanted to illustrate these methods with data on an important topic. That way, if the tutorial’s analysis improves our understanding at all, even if slightly, then our methods paper might also help us better understand a topic about that a lot of people care about.

Second, we wanted to analyze publicly available data. We searched a LOT of criminology papers that we thought would be good candidates and usually came up empty at this juncture. So, THANK YOU to Pickett, Graham, and Cullen for practicing open science and making this analysis possible!

Third, we wanted to pick a paper that we thought was already well done - that is, a paper that did a good job of carefully and thoroughly applying traditional tools to extract rich, meaningful signals from data. In doing so, it would make our job easier as we would already have a good sense of signals in the data and we could avoid taking too many methodological detours. This approach might also help readers better assess what, if anything, might be gained by spending valuable time learning new methods compared to traditional practices published in our top journals.

Finally, our main goal was to find data for a study in which researchers analyze Likert-type ordinal survey items as outcomes and combine them into composite scales so we could illustrate how ordinal methods might help us better understand signals in our data compared to traditional approaches in this common scenario. The methods can be - and we think often should be - applied to other types of data, including rare count outcomes (e.g., self-reported crime) and even continuous outcomes with thousands of levels! (See Frank Harrell, e.g. here and here).

2 Load libraries & settings

I will start by loading various libraries that I often use for such projects.

# Packages I need
library(here)
library(haven)
library(tidyverse)
library(gt)
library(janitor)
library(vtable)
library(modelsummary)
library(easystats)
library(marginaleffects)
library(ggplot2)
library(GGally)
library(ggeffects)
library(ggdist)
library(patchwork)
library(parallelly)
library(bayesplot)
library(ggdist)
library(correlation)
library(sjmisc)
library(forcats)
library(ggtext)
library(ggmosaic)
library(gtExtras)
library(tidybayes)
library(scales)
library(purrr)
library(viridis)
library(loo)
library(broom)  

#need MASS for frequentist polr but it masks common dplyr functions 
#use explicit dplyr::select() and dplyr::filter calls 
# or don't load MASS & call MASS::polr() instead 
library(MASS)  

library(extrafont)
library(stringr)
library(ragg)

# to use brms, do following:

# install rtools 4.4 

# install rstan
# https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started

# test ability to install packages from sources using Rtools: 
# install.packages("Rcpp", type="source")
  # it works

# remove any existing rstan installs
# remove.packages("rstan")
# if (file.exists(".RData")) file.remove(".RData")

# run the next line if you already have rstan installed
# remove.packages(c("StanHeaders", "rstan"))

# install rstan
# install.packages("rstan", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))

# test rstan install
# example(stan_model, package = "rstan", run.dontrun = TRUE)
  # it works 

library(rstan)

# install cmdstanr (alternative brms backend for interfacing w/stan)
# install.packages("cmdstanr", repos = c('https://stan-dev.r-universe.dev', getOption("repos")))

# newest developmental version
# install.packages("remotes")
# remotes::install_github("stan-dev/cmdstanr", force=TRUE)

# Sys.getenv("RTOOLS44_HOME")
# Sys.getenv("PATH")
# readLines("~/.Renviron")
  # had to delete old .Reviron fie pointing to Rtools40 from Documents folder 

# if prefer cmdstanr backend to brms, then: 
# cmdstanr::install_cmdstan()
# cmdstanr::check_cmdstan_toolchain(fix = TRUE)

# library(cmdstanr)


library(rstanarm)
library(brms)

print("T/F: Root 'here()' folder contains subfolder 'Models'")

[1] "T/F: Root 'here()' folder contains subfolder 'Models'"

#check for Models folder to save figures, create if one does not exist
ifelse(dir.exists(here("Models")), TRUE, dir.create(here("Models")))

[1] TRUE

I will also set some default settings.

#Fonts
# First time setup (run this once in console, not in your document):
# Import and register fonts (only needed once per system)
# Uncomment if you haven't imported fonts before
# extrafont::font_import()

# Load the fonts into the R session
extrafont::loadfonts(device = "win", quiet = TRUE)  # For Windows
# extrafont::loadfonts(device = "pdf", quiet = TRUE)  # For PDF output
# extrafont::loadfonts(device = "postscript", quiet = TRUE)  # For PostScript output

# If you only imported Times New Roman, load this instead: 
# Load the specific font you need for this session
# if(!("Times New Roman" %in% windowsFonts())) {
#   windowsFonts(`Times New Roman` = windowsFont("Times New Roman"))
# }

#To improve font rendering in viewer:
# Configure knitr to use the ragg device which handles fonts better
knitr::opts_chunk$set(
  dev = "ragg_png",
  dpi = 300
)

# Set the default theme for all ggplot2 plots
theme_set(theme_minimal(base_family = "Times New Roman"))


#set default knit & scientific notation options
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE, fig.align="center")
options(scipen = 999, digits = 2) #set global option to two decimals 

# identify & set # of available cores
  # replaced parallel::detectCores()
  # https://www.jottr.org/2022/12/05/avoid-detectcores/
nCoresphys <- parallelly::availableCores(logical = FALSE) 

options(mc.cores = parallelly::availableCores(logical=FALSE))
#options(mc.cores = 4, logical=FALSE)
  #alt: can specify exact number of cores desired here &/or in models/ppchecks   

color_scheme_set("viridisA")
# https://mc-stan.org/bayesplot/reference/bayesplot-colors.html

3 Simulation: Figure 1

Before getting into real data, let’s illustrate the underlying issue using a simulation.

I will simulate five different ordinal distributions with distinct shapes. Think of them as five groups, with N=250 observations per distributional group (N=1,250 total).

# Set seed for reproducibility
set.seed(1138)

# Number of observations per group
n <- 250

# Function to simulate ordinal data based on probabilities
simulate_ordinal <- function(probs, n) {
  sample(1:length(probs), size = n, replace = TRUE, prob = probs)
}

# Distribution 0: Symmetric Normal (peak at 3)
probs_symmetric <- c(0.10, 0.20, 0.40, 0.20, 0.10)  # Symmetric around 3
data_symmetric <- simulate_ordinal(probs_symmetric, n)

# Original distributions
probs_uniform <- rep(0.2, 5)
data_uniform <- simulate_ordinal(probs_uniform, n)

probs_bimodal <- c(0.3, 0.1333, 0.1333, 0.1333, 0.3)
data_bimodal <- simulate_ordinal(probs_bimodal, n)

probs_skew_left <- c(0.10, 0.20, 0.15, 0.25, 0.30)
data_skew_left <- simulate_ordinal(probs_skew_left, n)

probs_skew_right <- c(0.30, 0.25, 0.15, 0.20, 0.10)
data_skew_right <- simulate_ordinal(probs_skew_right, n)

# Combine data into a data frame
df <- data.frame(
  group = factor(rep(c("1. Symmetric Normal", "2. Uniform", "3. Bimodal", 
                      "4. Skewed Left", "5. Skewed Right"), 
                    each = n),
                levels = c("1. Symmetric Normal", "2. Uniform", "3. Bimodal", 
                          "4. Skewed Left", "5. Skewed Right")),
  response = c(data_symmetric, data_uniform, data_bimodal, data_skew_left, data_skew_right)
)

# Function to get normal curve points for each group separately
get_normal_curve <- function(data, group_name, x_range = seq(0.5, 5.5, length.out = 100)) {
  mean_val <- mean(data)
  sd_val <- sd(data)
  y_vals <- dnorm(x_range, mean = mean_val, sd = sd_val)
  scaling_factor <- (max(table(data))/length(data)) / max(y_vals)
  data.frame(
    x = x_range, 
    y = y_vals * scaling_factor * 100,
    group = group_name
  )
}

# Function to get ordered probit predicted probabilities with CIs
get_probit_probs <- function(data, group_name, n_boot = 1000) {
  model <- polr(ordered(response) ~ 1, data = data.frame(response = data), method = "probit")
  probs <- predict(model, newdata = data.frame(response = 1), type = "probs")
  
  boot_probs <- matrix(nrow = n_boot, ncol = 5)
  for(i in 1:n_boot) {
    boot_data <- sample(data, replace = TRUE)
    boot_model <- try(polr(ordered(boot_data) ~ 1, method = "probit"), silent = TRUE)
    if(!inherits(boot_model, "try-error")) {
      boot_probs[i,] <- predict(boot_model, newdata = data.frame(response = 1), type = "probs")
    }
  }
  
  ci_lower <- apply(boot_probs, 2, quantile, 0.025, na.rm = TRUE)
  ci_upper <- apply(boot_probs, 2, quantile, 0.975, na.rm = TRUE)
  
  data.frame(
    x = 1:5,
    y = as.numeric(probs) * 100,
    lower = ci_lower * 100,
    upper = ci_upper * 100,
    group = group_name
  )
}

# Get curves for each group
curves_df <- bind_rows(
  get_normal_curve(data_symmetric, "1. Symmetric Normal"),
  get_normal_curve(data_uniform, "2. Uniform"),
  get_normal_curve(data_bimodal, "3. Bimodal"),
  get_normal_curve(data_skew_left, "4. Skewed Left"),
  get_normal_curve(data_skew_right, "5. Skewed Right")
)

# Get probit probabilities with CIs for each group
probit_df <- bind_rows(
  get_probit_probs(data_symmetric, "1. Symmetric Normal"),
  get_probit_probs(data_uniform, "2. Uniform"),
  get_probit_probs(data_bimodal, "3. Bimodal"),
  get_probit_probs(data_skew_left, "4. Skewed Left"),
  get_probit_probs(data_skew_right, "5. Skewed Right")
)


# Update the stats_df creation to prepend A. through E.
stats_df <- df %>%
  group_by(group) %>%
  summarise(
    mean = round(mean(response), 2),
    sd = round(sd(response), 2)
  ) %>%
  mutate(
    # Create new labels with A. or B. etc. prepended
    new_label = case_when(
      group == "1. Symmetric Normal" ~ paste0("A. Normal-like\n(mean=", mean, ", sd=", sd, ")"),
      TRUE ~ paste0(
        letters[match(group, unique(group)) + 0] %>% toupper(), 
        ". ", 
        str_remove(group, "^\\d+\\.\\s+"), 
        "\n(mean=", mean, ", sd=", sd, ")"
      )
    )
  )


# Create two separate plots

# Plot 1: Normal-approx
p1 <- ggplot(df %>% dplyr::filter(group == "1. Symmetric Normal"), aes(x = response)) +
  geom_bar(aes(y = after_stat(prop)), fill = "#FEC488FF", color = "#FEC488FF") +
  geom_line(data = curves_df %>% dplyr::filter(group == "1. Symmetric Normal") %>% mutate(y = y/100), 
            aes(x = x, y = y), color = "#341069FF", linewidth = 1.5) +
  geom_pointrange(data = probit_df %>% dplyr::filter(group == "1. Symmetric Normal") %>% 
                   mutate(y = y/100, lower = lower/100, upper = upper/100), 
                 aes(x = x, y = y, ymin = lower, ymax = upper),
                 color = "#D3436EFF", size = 1.5, linewidth = 1, fatten = 2) +
  facet_wrap(~group, ncol = 1, labeller = as_labeller(setNames(stats_df$new_label, stats_df$group))) +
  scale_x_continuous(breaks = 1:5) +
  scale_y_continuous(limits = c(0, 0.45), 
                    expand = expansion(mult = c(0.02, 0.08))) +
  labs(x = "Response Category", y = "Proportion") +
  theme_minimal() +
  theme(
        text = element_text(family = "serif"),
        axis.title.x = element_text(size=12),  
        axis.title.y = element_text(size=12),  
        axis.text.x = element_text(size=12),  
        axis.text.y = element_text(size=12),  
        strip.text = element_text(hjust = 0, size=14), 
        panel.grid = element_blank()
        )

# Plot 2: Other distributions in 2x2 grid
p2 <- ggplot(df %>% dplyr::filter(group != "1. Symmetric Normal"), aes(x = response)) +
  geom_bar(aes(y = after_stat(prop)), fill = "#FEC488FF", color = "#FEC488FF") +
  geom_line(data = curves_df %>% dplyr::filter(group != "1. Symmetric Normal") %>% mutate(y = y/100), 
            aes(x = x, y = y), color = "#341069FF", linewidth = 1.5) +
  geom_pointrange(data = probit_df %>% dplyr::filter(group != "1. Symmetric Normal") %>% 
                   mutate(y = y/100, lower = lower/100, upper = upper/100), 
                 aes(x = x, y = y, ymin = lower, ymax = upper),
                 color = "#D3436EFF", size = 1.5, linewidth = 1, fatten = 2) +
  facet_wrap(~group, ncol = 2, labeller = as_labeller(setNames(stats_df$new_label, stats_df$group))) +
  scale_x_continuous(breaks = 1:5) +
  scale_y_continuous(limits = c(0, 0.45), 
                    expand = expansion(mult = c(0.02, 0.08))) +
  labs(x = "Response Category", y = "") +
  theme_minimal() +
  theme(
        text = element_text(family = "serif"),
        axis.title.x = element_text(size=12),  
        axis.title.y = element_text(size=12),  
        axis.text.x = element_text(size=12),  
        axis.text.y = element_text(size=12),  
        strip.text = element_text(hjust = 0, size=14), 
        panel.grid = element_blank())

# Combine plots using patchwork
# combined_plot <- p1 + p2 + 
#   plot_layout(widths = c(1, 2)) +
#   plot_annotation(
#     title = "FIGURE1. Comparing Metric vs. Ordinal View of Data",
#     subtitle = "Data: Light bars show observed (simulated) proportions\nMetric view: Dark line shows linear model-implied normal curve\nOrdinal view: Point-intervals show ordered probit model-implied predicted probabilities & 95% CI",
#     theme = theme(
#       plot.title = element_text(hjust = 0.5),
#       plot.subtitle = element_text(hjust = 0)
#     )
#   )


# Create spacers for above/below normal-like
top_spacer <- plot_spacer()
bottom_spacer <- plot_spacer()

# Combine plots w/custom layout
combined_plot <- top_spacer + p2 + 
                 p1 + bottom_spacer + 
  plot_layout(
    design = "
    AB
    CB
    DB
    ",
    heights = c(1, 2, 1), 
    widths = c(1, 2)
  ) +
  plot_annotation(
    title = "FIGURE1. Comparing Metric vs. Ordinal View of Data",
    subtitle = "Data: Light bars show observed (simulated) proportions\nMetric view: Dark line shows linear model-implied normal curve\nOrdinal view: Point-intervals show ordered probit model-implied predicted probabilities & 95% CI",
    theme = theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0)
    )
  )

3.1 FIG1

combined_plot

Clearly, these distributions are very different. However, a linear model would predict no mean differences between three of these groups (A, B, C). Meanwhile, although a linear model might detect differences between the two skewed distributions (D vs. E) with a large enough sample, the model-implied differences on the mean scale are hard to interpret, and predicted differences on response probability scale may depart drastically from observed differences in the data. For example, the linear model “view” only sees about ~1 response category difference between D & E, with means implying the typical observation in D is a 3 or 4 and the typical observation in E is a 2 or 3. Yet, the majority of observations in D were in categories 4 & 5 (not 3 & 4), whereas the majority of observations in E were in 1 & 2 (not 2 & 3)!

Before moving on, it might help to visualize the latent response thresholds.

# Define the threshold annotation function
add_probit_thresholds <- function(plot) {
  plot +
    # Add vertical lines at each response category boundary
    geom_vline(xintercept = 1:4 + 0.5, linetype = "dashed", color = "gray60", alpha = 0.5) +
    annotate("text", x = 2.5, y = 0.4, 
             label = "Cumulative\nThresholds", 
             color = "gray50", 
             hjust = 0.5)
}

# Modify both p1 and p2 to add threshold labels
p1 <- ggplot(df %>% dplyr::filter(group == "1. Symmetric Normal"), aes(x = response)) +
  geom_bar(aes(y = after_stat(prop)), fill = "#FEC488FF", color = "#FEC488FF") +
  geom_line(data = curves_df %>% dplyr::filter(group == "1. Symmetric Normal") %>% mutate(y = y/100), 
            aes(x = x, y = y), color = "#341069FF", linewidth = 1.5) +
  geom_pointrange(data = probit_df %>% dplyr::filter(group == "1. Symmetric Normal") %>% 
                   mutate(y = y/100, lower = lower/100, upper = upper/100), 
                 aes(x = x, y = y, ymin = lower, ymax = upper),
                 color = "#D3436EFF", size = 1.5, linewidth = 1, fatten = 2) +
  # Add vertical threshold lines with tau labels
  geom_vline(xintercept = 1:4 + 0.5, linetype = "dashed", color = "gray60", alpha = 0.5) +
  annotate("text", x=1.65, y=0.42, label = "tau[1]", color = "gray50", 
           parse = TRUE, size=6) +
  annotate("text", x=2.65, y=0.42, label = "tau[2]", color = "gray50", 
           parse = TRUE, size=6) +
  annotate("text", x=3.65, y=0.42, label = "tau[3]", color = "gray50", 
           parse = TRUE, , size=6) +
  annotate("text", x=4.65, y=0.42, label = "tau[4]", color = "gray50", 
           parse = TRUE, , size=6) +
  facet_wrap(~group, ncol = 1, labeller = as_labeller(setNames(stats_df$new_label, stats_df$group))) +
  scale_x_continuous(breaks = 1:5) +
  scale_y_continuous(limits = c(0, 0.45), 
                    expand = expansion(mult = c(0.02, 0.08))) +
  labs(x = "Response Category", y = "Proportion") +
  theme_minimal() +
  theme(
        text = element_text(family = "serif"),
        strip.text = element_text(hjust = 0, size=14), 
        panel.grid = element_blank(), 
        axis.text = element_text(size=14),
        axis.title = element_text(size=14)
        )

p2 <- ggplot(df %>% dplyr::filter(group != "1. Symmetric Normal"), aes(x = response)) +
  geom_bar(aes(y = after_stat(prop)), fill = "#FEC488FF", color = "#FEC488FF") +
  geom_line(data = curves_df %>% dplyr::filter(group != "1. Symmetric Normal") %>% mutate(y = y/100), 
            aes(x = x, y = y), color = "#341069FF", linewidth = 1.5) +
  geom_pointrange(data = probit_df %>% dplyr::filter(group != "1. Symmetric Normal") %>% 
                   mutate(y = y/100, lower = lower/100, upper = upper/100), 
                 aes(x = x, y = y, ymin = lower, ymax = upper),
                 color = "#D3436EFF", size = 1.5, linewidth = 1, fatten = 2) +
  # Add vertical threshold lines with tau labels
  geom_vline(xintercept = 1:4 + 0.5, linetype = "dashed", color = "gray60", alpha = 0.5) +
  annotate("text", x=1.65, y=0.42, label = "tau[1]", color = "gray50", 
           parse = TRUE, size=6) +
  annotate("text", x=2.65, y=0.42, label = "tau[2]", color = "gray50", 
           parse = TRUE, size=6) +
  annotate("text", x=3.65, y=0.42, label = "tau[3]", color = "gray50", 
           parse = TRUE, , size=6) +
  annotate("text", x=4.65, y=0.42, label = "tau[4]", color = "gray50", 
           parse = TRUE, , size=6) +
facet_wrap(~group, ncol = 2, labeller = as_labeller(setNames(stats_df$new_label, stats_df$group))) +
  scale_x_continuous(breaks = 1:5) +
  scale_y_continuous(limits = c(0, 0.45), 
                    expand = expansion(mult = c(0.02, 0.08))) +
  labs(x = "Response Category", y = "") +
  theme_minimal() +
  theme(
        text = element_text(family = "serif"),
        strip.text = element_text(hjust = 0.5, size=14), 
        panel.grid = element_blank(), 
        axis.text = element_text(size=14),
        axis.title = element_text(size=14)
)

combined_plot <- p1 + p2 + 
  plot_layout(widths = c(1, 2)) +
  plot_annotation(
    title = "Comparing Metric vs. Ordinal View of Data",
    subtitle = "Data: Light bars show observed (simulated) proportions\nMetric view: Dark line shows linear model-implied normal curve\nOrdinal view: Point-intervals show ordered probit model predicted probabilities & 95% CI\n                      Dashed lines represent category thresholds in latent space",
    theme = theme(
      plot.title = element_text(hjust = 0.5),
      plot.subtitle = element_text(hjust = 0))
  )

combined_plot

Following Kurz’s brms adaptation of Kruschke, we can also visualize the latent thresholds by fitting Bayesian threshold-only models and then plotting draws from the posterior distribution.

# Function to fit and save brms model for each distribution
fit_brms_model <- function(data, group_name) {
  # Prepare the data
  my_data <- data.frame(response = data)
  
  # Fit a Bayesian cumulative probit model
  fit <- brm(
    data = my_data,
    family = cumulative(probit),
    response ~ 1,
    prior = prior(normal(0, 4), class = Intercept),
    iter = 3000, warmup = 1000, chains = 4, cores = 4,
    seed = 123,
    file = paste0("Models/sim_fit_", group_name),
    file_refit = "on_change"
  )
  
  return(fit)
}

# Fit models for each distribution
fit_symmetric <- fit_brms_model(data_symmetric, "symmetric")
fit_uniform <- fit_brms_model(data_uniform, "uniform")
fit_bimodal <- fit_brms_model(data_bimodal, "bimodal")
fit_skew_left <- fit_brms_model(data_skew_left, "skew_left")
fit_skew_right <- fit_brms_model(data_skew_right, "skew_right")

# Create a list of fits for easy reference
fits_list <- list(
  symmetric = fit_symmetric,
  uniform = fit_uniform,
  bimodal = fit_bimodal,
  skew_left = fit_skew_left,
  skew_right = fit_skew_right
)

# Function to create threshold plot
create_threshold_plot <- function(fit, group_name) {
  # Extract draws 
  draws <- as_draws_df(fit)
  
  # Compute the posterior means for each threshold
  means <- draws %>% 
    summarise(across(starts_with("b_Intercept"), mean)) %>% 
    pivot_longer(everything(), names_to = "name", values_to = "mean")
  
  # Prepare the data with explicit threshold calculation
  plot_data <- draws %>% 
    pivot_longer(starts_with("b_Intercept"), names_to = "name", values_to = "threshold") %>% 
    group_by(.draw) %>% 
    # Calculate theta_bar as the mean of all thresholds for this draw
    mutate(theta_bar = mean(threshold))
  
  # Create the plot
  ggplot(plot_data, aes(x = threshold, y = theta_bar, color = name)) +
    # Vertical lines at mean thresholds
    geom_vline(data = means,
               aes(xintercept = mean, color = name),
               linetype = 2) +
    # Points with low alpha to show density
    geom_point(alpha = 1/10) +
    scale_color_viridis_d(option = "magma", begin = 0.2, end = 0.8) +
    # Dynamic x-axis limits based on actual threshold values
    scale_x_continuous(
      breaks = means$mean,
      labels = c("τ1", "τ2", "τ3", "τ4")
    ) +
    labs(title = group_name,
         y = "Mean of Thresholds per Draw",
         x = "Latent Thresholds") +
    theme_minimal() +
    theme(
      legend.position = "none", 
      axis.text.y = element_blank(), 
      axis.ticks.y = element_blank()
    )
}
# Create threshold plots for each model
threshold_plots <- list(
  create_threshold_plot(fit_symmetric, "A. Normal-like"),
  create_threshold_plot(fit_uniform, "B. Uniform"),
  create_threshold_plot(fit_bimodal, "C. Bimodal"),
  create_threshold_plot(fit_skew_left, "D. Skewed Left"),
  create_threshold_plot(fit_skew_right, "E. Skewed Right")
)

# Combine the threshold plots
combined_threshold_plot <- wrap_plots(threshold_plots, ncol = 3)
combined_threshold_plot

4 Supp. Sim. S1a: Distributional Change w/Ideological Polarization

This simulation demonstrates a scenario where linear and ordinal models yield different substantive conclusions despite zero overall average treatment effects. We model a situation where a highly publicized police incident causes polarized reactions that vary by political ideology, illustrating how distributional changes can be obscured when focusing solely on average effects.

The key insight: Even when ordinal data appears “continuous-like” with symmetric distributions across many categories, the underlying categorical social processes can create meaningful distributional shifts that the cumulative ordinal model captures more accurately compared to the linear model - even with explicit proportional odds assumption violation.

4.1 Simulation Design

set.seed(1138)
n_per_group <- 1000

# Create symmetric baseline distributions that look "normal enough" for OLS
# This mimics what researchers see when they think linear regression should be fine
generate_baseline <- function(n, ideology) {
  # Both groups start with nearly identical, symmetric 7-point distributions
  if(ideology == "Conservative") {
    probs <- c(0.06, 0.11, 0.21, 0.28, 0.19, 0.10, 0.05)
  } else {
    probs <- c(0.05, 0.10, 0.19, 0.28, 0.21, 0.11, 0.06)
  }
  
  responses <- sample(1:7, size = n, replace = TRUE, prob = probs)
  data.frame(
    id = 1:n,
    ideology = ideology,
    time = 0,
    response = responses
  )
}

baseline_data <- rbind(
  generate_baseline(n_per_group, "Conservative"),
  generate_baseline(n_per_group, "Liberal")
)

# Generate post-incident responses with BALANCED group-specific polarization
generate_post_incident <- function(baseline_df, shock_intensity = 0.75) {
  post_data <- baseline_df
  post_data$time <- 1
  
  for(i in 1:nrow(post_data)) {
    if(runif(1) < shock_intensity) {
      current_response <- post_data$response[i]
      ideology <- post_data$ideology[i]
      
      # BALANCED polarization: same proportions, opposite directions
      if(ideology == "Conservative") {
        # Conservatives: 80% move toward "unafraid", 20% toward "afraid"
        move_toward_unafraid <- runif(1) < 0.80
      } else {
        # Liberals: 20% move toward "unafraid", 80% toward "afraid"  
        move_toward_unafraid <- runif(1) < 0.20
      }
      
      if(move_toward_unafraid) {
        # Push toward extreme "unafraid" categories (1-2)
        if(current_response >= 5) {
          post_data$response[i] <- ifelse(runif(1) < 0.7, 1, 2)
        } else if(current_response == 4) {
          post_data$response[i] <- ifelse(runif(1) < 0.6, 1, 2)
        } else {
          post_data$response[i] <- 1
        }
      } else {
        # Push toward extreme "afraid" categories (6-7)
        if(current_response <= 3) {
          post_data$response[i] <- ifelse(runif(1) < 0.7, 7, 6)
        } else if(current_response == 4) {
          post_data$response[i] <- ifelse(runif(1) < 0.6, 7, 6)
        } else {
          post_data$response[i] <- 7
        }
      }
    }
  }
  return(post_data)
}

post_data <- generate_post_incident(baseline_data, shock_intensity = 0.75)

# Combine datasets
sim_data <- rbind(baseline_data, post_data) %>%
  mutate(
    ideology = factor(ideology, levels = c("Conservative", "Liberal")),
    response_ord = factor(response, levels = 1:7, ordered = TRUE)
  )

4.2 Descriptive Analysis

# Calculate group-specific means
means_summary <- sim_data %>%
  group_by(ideology, time) %>%
  summarise(mean_response = mean(response), .groups = "drop") %>%
  pivot_wider(names_from = time, values_from = mean_response, names_prefix = "time_") %>%
  mutate(change = time_1 - time_0)

print(means_summary)

# A tibble: 2 × 4
  ideology     time_0 time_1 change
  <fct>         <dbl>  <dbl>  <dbl>
1 Conservative   3.88   2.71  -1.17
2 Liberal        4.07   5.25   1.18

# Overall treatment effect
overall_ate <- sim_data %>%
  group_by(time) %>%
  summarise(mean_resp = mean(response), .groups = "drop") %>%
  summarise(ate = diff(mean_resp)) %>%
  pull(ate)

cat("Overall Average Treatment Effect:", round(overall_ate, 3), "\n")

Overall Average Treatment Effect: 0.005 

# Show distributional changes by group
dist_summary <- sim_data %>%
  group_by(ideology, time, response) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(ideology, time) %>%
  mutate(proportion = count / sum(count)) %>%
  ungroup()

# Create distribution comparison table
dist_table <- dist_summary %>%
  dplyr::select(ideology, time, response, proportion) %>%
  pivot_wider(names_from = c(ideology, time), values_from = proportion, values_fill = 0) %>%
  mutate(
    conservative_change = Conservative_1 - Conservative_0,
    liberal_change = Liberal_1 - Liberal_0
  )

print(round(dist_table, 3))

# A tibble: 7 × 7
  response Conservative_0 Conservative_1 Liberal_0 Liberal_1 conservative_change
     <dbl>          <dbl>          <dbl>     <dbl>     <dbl>               <dbl>
1        1          0.054          0.488     0.045     0.131               0.434
2        2          0.118          0.155     0.096     0.062               0.037
3        3          0.228          0.054     0.191     0.041              -0.174
4        4          0.278          0.08      0.294     0.071              -0.198
5        5          0.189          0.043     0.22      0.066              -0.146
6        6          0.079          0.045     0.095     0.143              -0.034
7        7          0.054          0.135     0.059     0.486               0.081
# ℹ 1 more variable: liberal_change <dbl>

4.3 Visualization

# Create visualization showing distributional changes by group
plot_data <- dist_summary %>%
  mutate(
    time_label = factor(time, levels = c(0, 1), labels = c("Pre-Incident", "Post-Incident")),
    response_label = factor(response, levels = 1:7, 
                           labels = c("Very\nUnafraid", "Unafraid", "Somewhat\nUnafraid", 
                                    "Neutral", "Somewhat\nAfraid", "Afraid", "Very\nAfraid"))
  )

polarization_viz <- ggplot(plot_data, aes(x = response_label, y = proportion, fill = time_label)) +
  geom_col(position = "dodge", alpha = 0.8, color = "white", linewidth = 0.3) +
  facet_wrap(~ideology, labeller = label_both) +
  scale_fill_manual(values = c("Pre-Incident" = "#FD9B6BFF", "Post-Incident" = "#782281FF"),
                   name = "Period") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                    expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "Ideological Polarization Following Police Incident",
    subtitle = paste0("Overall ATE = ", round(overall_ate, 3), 
                     " | Conservative Δ = ", round(means_summary$change[means_summary$ideology == "Conservative"], 3),
                     " | Liberal Δ = ", round(means_summary$change[means_summary$ideology == "Liberal"], 3)),
    x = "Fear of Police Response",
    y = "Proportion of Responses"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 11),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    strip.text = element_text(size = 11, face = "bold")
  )

polarization_viz

4.4 Model Comparison

# Check if cached models exist
linear_model_path <- here("Models", "polarization_linear_model.rds")
ordinal_model_path <- here("Models", "polarization_ordinal_model.rds")

if (file.exists(linear_model_path)) {
  linear_model <- readRDS(linear_model_path)
} else {
  # Fit linear model (using numeric version for proper LOO comparison)
  linear_model <- brm(
    response ~ time * ideology,
    data = sim_data,
    family = gaussian(),
    prior = c(
      prior(normal(4, 1), class = "Intercept"),
      prior(normal(0, 0.5), class = "b"),
      prior(exponential(1), class = "sigma")
    ),
    chains = 4, iter = 2000, warmup = 1000,
    cores = 4, seed = 1138, refresh = 0
  )
  saveRDS(linear_model, linear_model_path)
}

if (file.exists(ordinal_model_path)) {
  ordinal_model <- readRDS(ordinal_model_path)
} else {
  # Fit ordinal model (using same numeric outcome for LOO comparison)
  ordinal_model <- brm(
    response ~ time * ideology,
    data = sim_data,
    family = cumulative("probit"),
    prior = c(
      prior(normal(-1.28, 1), class = "Intercept", coef = 1),
      prior(normal(-0.67, 1), class = "Intercept", coef = 2),
      prior(normal(-0.25, 1), class = "Intercept", coef = 3),
      prior(normal(0.25, 1), class = "Intercept", coef = 4),
      prior(normal(0.67, 1), class = "Intercept", coef = 5),
      prior(normal(1.28, 1), class = "Intercept", coef = 6),
      prior(normal(0, 0.5), class = "b")
    ),
    chains = 4, iter = 2000, warmup = 1000,
    cores = 4, seed = 1138, refresh = 0
  )
  saveRDS(ordinal_model, ordinal_model_path)
}

# Extract key effects
linear_effects <- fixef(linear_model)
ordinal_effects <- fixef(ordinal_model)

# Compare main effects and interactions
effect_comparison <- tibble(
  Effect = c("Time (Main)", "Ideology (Liberal)", "Time × Ideology Interaction"),
  Linear_Est = round(c(linear_effects["time", "Estimate"],
                      linear_effects["ideologyLiberal", "Estimate"], 
                      linear_effects["time:ideologyLiberal", "Estimate"]), 3),
  Linear_CI = paste0("[", round(c(linear_effects["time", "Q2.5"],
                                 linear_effects["ideologyLiberal", "Q2.5"],
                                 linear_effects["time:ideologyLiberal", "Q2.5"]), 3), ", ",
                    round(c(linear_effects["time", "Q97.5"],
                           linear_effects["ideologyLiberal", "Q97.5"], 
                           linear_effects["time:ideologyLiberal", "Q97.5"]), 3), "]"),
  Ordinal_Est = round(c(ordinal_effects["time", "Estimate"],
                       ordinal_effects["ideologyLiberal", "Estimate"],
                       ordinal_effects["time:ideologyLiberal", "Estimate"]), 3),
  Ordinal_CI = paste0("[", round(c(ordinal_effects["time", "Q2.5"],
                                  ordinal_effects["ideologyLiberal", "Q2.5"],
                                  ordinal_effects["time:ideologyLiberal", "Q2.5"]), 3), ", ",
                     round(c(ordinal_effects["time", "Q97.5"],
                            ordinal_effects["ideologyLiberal", "Q97.5"],
                            ordinal_effects["time:ideologyLiberal", "Q97.5"]), 3), "]")
)

gt(effect_comparison) %>%
  tab_header(
    title = "Model Comparison: Linear vs. Ordinal Estimates",
    subtitle = "Coefficients and 95% credible intervals"
  ) %>%
  cols_label(
    Linear_Est = "Estimate",
    Linear_CI = "95% CrI",
    Ordinal_Est = "Estimate", 
    Ordinal_CI = "95% CrI"
  ) %>%
  tab_spanner(label = "Linear Model", columns = c(Linear_Est, Linear_CI)) %>%
  tab_spanner(label = "Ordinal Model", columns = c(Ordinal_Est, Ordinal_CI)) %>%
  tab_style(style = cell_text(weight = "bold"), locations = cells_column_labels()) %>%
  tab_style(style = cell_text(align = "center"), 
           locations = cells_body(columns = c(Linear_Est, Linear_CI, Ordinal_Est, Ordinal_CI)))

Model Comparison: Linear vs. Ordinal Estimates
Coefficients and 95% credible intervals
Effect	Linear Model		Ordinal Model
	Estimate	95% CrI	Estimate	95% CrI
Time (Main)	-1.08	[-1.238, -0.923]	-0.69	[-0.789, -0.602]
Ideology (Liberal)	0.26	[0.1, 0.419]	0.10	[0.01, 0.191]
Time × Ideology Interaction	2.20	[1.983, 2.413]	1.37	[1.235, 1.504]

4.5 Predicted Probabilities: Revealing the Polarization

# Generate predictions for all ideology × time combinations
newdata <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  time = c(0, 1)
)

# Get predicted probabilities from ordinal model
ordinal_preds <- posterior_epred(ordinal_model, newdata = newdata, re_formula = NA)

# Calculate summary statistics
pred_summary <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  time = c(0, 1),
  category = 1:7
) %>%
  mutate(
    row_idx = case_when(
      ideology == "Conservative" & time == 0 ~ 1,
      ideology == "Conservative" & time == 1 ~ 2,
      ideology == "Liberal" & time == 0 ~ 3,
      ideology == "Liberal" & time == 1 ~ 4
    ),
    prob = map2_dbl(row_idx, category, ~mean(ordinal_preds[, .x, .y])),
    lower = map2_dbl(row_idx, category, ~quantile(ordinal_preds[, .x, .y], 0.025)),
    upper = map2_dbl(row_idx, category, ~quantile(ordinal_preds[, .x, .y], 0.975)),
    time_label = factor(time, levels = c(0, 1), labels = c("Pre-Incident", "Post-Incident")),
    category_label = factor(category, levels = 1:7,
                           labels = c("Very\nUnafraid", "Unafraid", "Somewhat\nUnafraid", 
                                    "Neutral", "Somewhat\nAfraid", "Afraid", "Very\nAfraid"))
  )

# Plot predicted probabilities
prob_plot <- pred_summary %>%
  ggplot(aes(x = category_label, y = prob, color = time_label, shape = time_label)) +
  geom_point(size = 2.5, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.15, position = position_dodge(width = 0.4)) +
  facet_wrap(~ideology, labeller = label_both) +
  scale_color_manual(values = c("Pre-Incident" = "#FD9B6BFF", "Post-Incident" = "#782281FF"),
                    name = "Period") +
  scale_shape_manual(values = c("Pre-Incident" = 16, "Post-Incident" = 17),
                    name = "Period") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Ordinal Model Predicted Probabilities",
    subtitle = "Ideology-specific polarization patterns revealed through category probabilities",
    x = "Fear Response Category",
    y = "Predicted Probability"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 11),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    strip.text = element_text(size = 11, face = "bold")
  )

prob_plot

4.6 Model Fit Assessment

# Compare predictive performance
loo_linear <- loo(linear_model)
loo_ordinal <- loo(ordinal_model)

# Manual comparison since loo_compare might have issues
elpd_diff <- loo_ordinal$estimates["elpd_loo", "Estimate"] - loo_linear$estimates["elpd_loo", "Estimate"]
se_diff <- sqrt(loo_ordinal$estimates["elpd_loo", "SE"]^2 + loo_linear$estimates["elpd_loo", "SE"]^2)

cat("ELPD Linear:", round(loo_linear$estimates["elpd_loo", "Estimate"], 1), "\n")

ELPD Linear: -8190 

cat("ELPD Ordinal:", round(loo_ordinal$estimates["elpd_loo", "Estimate"], 1), "\n")

ELPD Ordinal: -7229 

cat("ELPD Difference (Ordinal - Linear):", round(elpd_diff, 1), "\n")

ELPD Difference (Ordinal - Linear): 962 

cat("Standard Error:", round(se_diff, 1), "\n")

Standard Error: 58 

cat("Difference in standard errors:", round(abs(elpd_diff)/se_diff, 1), "\n")

Difference in standard errors: 17 

if(abs(elpd_diff) > 2 * se_diff) {
  cat("\n*** SUBSTANTIAL EVIDENCE favoring", ifelse(elpd_diff > 0, "ORDINAL", "LINEAR"), "model ***\n")
} else {
  cat("\nNo strong evidence favoring either model\n")
}

*** SUBSTANTIAL EVIDENCE favoring ORDINAL model ***

# Additional diagnostics
cat("\nDetailed LOO comparison:\n")

Detailed LOO comparison:

loo_comp <- loo_compare(loo_linear, loo_ordinal)
print(loo_comp)

              elpd_diff se_diff
ordinal_model    0.0       0.0 
linear_model  -961.5      36.7 

# Check if models are actually different by comparing log-likelihoods
ll_linear <- log_lik(linear_model)
ll_ordinal <- log_lik(ordinal_model)

cat("\nMean log-likelihood (linear):", round(mean(rowSums(ll_linear)), 2), "\n")

Mean log-likelihood (linear): -8188 

cat("Mean log-likelihood (ordinal):", round(mean(rowSums(ll_ordinal)), 2), "\n")

Mean log-likelihood (ordinal): -7224 

# Alternative: Compare WAIC
waic_linear <- waic(linear_model)
waic_ordinal <- waic(ordinal_model)
waic_comp <- loo_compare(waic_linear, waic_ordinal)

cat("\nWAIC comparison:\n")

WAIC comparison:

print(waic_comp)

              elpd_diff se_diff
ordinal_model    0.0       0.0 
linear_model  -961.5      36.7 

# Check actual prediction accuracy using posterior predictive
# Generate predictions for held-out test set approach
set.seed(42)
test_indices <- sample(nrow(sim_data), size = 200)
test_data <- sim_data[test_indices, ]
train_data <- sim_data[-test_indices, ]

# Calculate prediction accuracy manually
linear_preds <- posterior_predict(linear_model, newdata = test_data)
ordinal_preds <- posterior_predict(ordinal_model, newdata = test_data)

# Round linear predictions to nearest integer for comparison
linear_preds_rounded <- round(linear_preds)
linear_preds_rounded[linear_preds_rounded < 1] <- 1
linear_preds_rounded[linear_preds_rounded > 7] <- 7

# Calculate accuracy
linear_accuracy <- mean(apply(linear_preds_rounded, 1, function(x) mean(x == test_data$response)))
ordinal_accuracy <- mean(apply(ordinal_preds, 1, function(x) mean(x == test_data$response)))

cat("\nPrediction accuracy on test set:\n")

Prediction accuracy on test set:

cat("Linear model:", round(linear_accuracy, 3), "\n")

Linear model: 0.17 

cat("Ordinal model:", round(ordinal_accuracy, 3), "\n")

Ordinal model: 0.19 

# Posterior predictive checks
pp_linear <- pp_check(linear_model, ndraws = 50) + 
  ggtitle("Linear Model: Posterior Predictive Check") +
  labs(subtitle = "Model assumes continuous normal distribution") +
  theme_minimal(base_family = "serif") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

pp_ordinal <- pp_check(ordinal_model, ndraws = 50, type = "bars") +
  ggtitle("Ordinal Model: Posterior Predictive Check") +
  labs(subtitle = "Model captures discrete categorical structure") +
  theme_minimal(base_family = "serif") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

# Combine pp checks
pp_combined <- pp_linear / pp_ordinal
pp_combined

4.7 The Critical Difference: Predicted Probabilities

While both models show similar coefficient patterns, they make fundamentally different predictions about response probabilities. This reveals why focusing solely on coefficient interpretation can be misleading.

# Generate predictions from both models for post-incident period
newdata_post <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  time = 1  # Post-incident only
)

# Get predictions from linear model (need to convert to probabilities)
linear_pred_draws <- posterior_epred(linear_model, newdata = newdata_post, re_formula = NA)

# Convert linear predictions to "probabilities" by assuming normal distribution
# and calculating probability of falling in each category range
convert_linear_to_probs <- function(linear_preds, categories = 1:7) {
  n_draws <- nrow(linear_preds)
  n_conditions <- ncol(linear_preds)
  
  # Initialize probability array
  prob_array <- array(dim = c(n_draws, n_conditions, length(categories)))
  
  for(draw in 1:n_draws) {
    for(cond in 1:n_conditions) {
      mean_pred <- linear_preds[draw, cond]
      # Get residual SD from model using modern syntax
      sigma <- as_draws_df(linear_model)$sigma[draw]
      
      # Calculate probability of each category
      for(cat in 1:length(categories)) {
        if(cat == 1) {
          prob_array[draw, cond, cat] <- pnorm(1.5, mean_pred, sigma)
        } else if(cat == length(categories)) {
          prob_array[draw, cond, cat] <- 1 - pnorm(cat - 0.5, mean_pred, sigma)
        } else {
          prob_array[draw, cond, cat] <- pnorm(cat + 0.5, mean_pred, sigma) - 
                                       pnorm(cat - 0.5, mean_pred, sigma)
        }
      }
    }
  }
  return(prob_array)
}

# Get ordinal model predictions
ordinal_pred_draws <- posterior_epred(ordinal_model, newdata = newdata_post, re_formula = NA)

# Convert linear predictions to probabilities
linear_probs <- convert_linear_to_probs(linear_pred_draws)

# Calculate observed proportions for post-incident period
observed_post <- sim_data %>%
  dplyr::filter(time == 1) %>%
  group_by(ideology, response) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(ideology) %>%
  mutate(proportion = count / sum(count))

# Create comparison dataset
prob_comparison <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  category = 1:7
) %>%
  mutate(
    condition_idx = case_when(
      ideology == "Conservative" ~ 1,
      ideology == "Liberal" ~ 2
    ),
    # Linear model predictions
    linear_prob = map2_dbl(condition_idx, category, 
                          ~mean(linear_probs[, .x, .y])),
    linear_lower = map2_dbl(condition_idx, category, 
                           ~quantile(linear_probs[, .x, .y], 0.025)),
    linear_upper = map2_dbl(condition_idx, category, 
                           ~quantile(linear_probs[, .x, .y], 0.975)),
    # Ordinal model predictions
    ordinal_prob = map2_dbl(condition_idx, category, 
                           ~mean(ordinal_pred_draws[, .x, .y])),
    ordinal_lower = map2_dbl(condition_idx, category, 
                            ~quantile(ordinal_pred_draws[, .x, .y], 0.025)),
    ordinal_upper = map2_dbl(condition_idx, category, 
                            ~quantile(ordinal_pred_draws[, .x, .y], 0.975)),
    category_label = factor(category, levels = 1:7,
                           labels = c("Very\nUnafraid", "Unafraid", "Somewhat\nUnafraid", 
                                    "Neutral", "Somewhat\nAfraid", "Afraid", "Very\nAfraid"))
  ) %>%
  # Add observed proportions
  left_join(
    observed_post %>% dplyr::select(ideology, response, proportion),
    by = c("ideology", "category" = "response")
  ) %>%
  # Reshape for plotting
  pivot_longer(
    cols = c(linear_prob, ordinal_prob, proportion),
    names_to = "type", 
    values_to = "probability"
  ) %>%
  mutate(
    lower = case_when(
      type == "linear_prob" ~ linear_lower,
      type == "ordinal_prob" ~ ordinal_lower,
      type == "proportion" ~ NA_real_
    ),
    upper = case_when(
      type == "linear_prob" ~ linear_upper,
      type == "ordinal_prob" ~ ordinal_upper,
      type == "proportion" ~ NA_real_
    ),
    type_label = factor(type, 
                       levels = c("proportion", "linear_prob", "ordinal_prob"),
                       labels = c("Observed", "Linear Model", "Ordinal Model"))
  )

# Plot model comparison with colorblind-friendly palette
model_comparison_plot <- prob_comparison %>%
  ggplot(aes(x = category_label, y = probability, color = type_label, shape = type_label)) +
  geom_point(size = 2.5, position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = lower, ymax = upper), 
                width = 0.15, position = position_dodge(width = 0.4),
                na.rm = TRUE) +
  geom_line(aes(group = type_label), position = position_dodge(width = 0.4), alpha = 0.7) +
  facet_wrap(~ideology, labeller = label_both) +
  scale_color_manual(values = c("Observed" = "#000000", "Linear Model" = "#D55E00", "Ordinal Model" = "#0072B2"),
                    name = "") +
  scale_shape_manual(values = c("Observed" = 17, "Linear Model" = 15, "Ordinal Model" = 16),
                    name = "") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "Model Accuracy: Observed vs. Predicted Probabilities",
    subtitle = "Ordinal model captures bimodal reality; linear model predicts impossible intermediate peaks",
    x = "Fear Response Category",
    y = "Probability"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 13, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 11),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 11),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    strip.text = element_text(size = 11, face = "bold")
  )

model_comparison_plot

4.8 Combined Visualization

# Panel A: Overall symmetric distribution of response variable
overall_dist <- sim_data %>%
  dplyr::filter(time == 0) %>%  # Pre-incident baseline
  group_by(response) %>%
  summarise(count = n(), .groups = "drop") %>%
  mutate(
    proportion = count / sum(count),
    response_label = factor(response, levels = 1:7,
                       labels = c("1\nVery\nUnafraid", "2\nUnafraid",
                                  "3\nSomewhat\nUnafraid", "4\nNeutral",
                                  "5\nSomewhat\nAfraid", "6\nAfraid",
                                  "7\nVery\nAfraid"))
  )

panel_a <- overall_dist %>%
  ggplot(aes(x = response_label, y = proportion)) +
  geom_col(fill = "#756bb1", alpha = 0.8, color = "white", linewidth = 0.3) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                    expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "A. Baseline Distribution",
    subtitle = "Symmetric 7-category response distribution",
    x = "Fear Response Category",
    y = "Proportion"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 10),
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank()
  )

# Panel B: Use the original plot_data that should have response_label defined
# If plot_data doesn't have proper labels, recreate it here
plot_data_fixed <- dist_summary %>%
  mutate(
    time_label = factor(time, levels = c(0, 1), labels = c("Pre-Incident", "Post-Incident")),
    response_label = factor(response, levels = 1:7, 
                       labels = c("1", "2", "3", "4", "5", "6", "7"))
  )

panel_b <- plot_data_fixed %>%
  ggplot(aes(x = response_label, y = proportion, fill = time_label)) +
  geom_col(position = "dodge", alpha = 0.8, color = "white", linewidth = 0.3) +
  facet_wrap(~ideology, labeller = label_both) +
  scale_fill_manual(values = c("Pre-Incident" = "#FD9B6BFF", "Post-Incident" = "#782281FF"),
                   name = "Period") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1),
                    expand = expansion(mult = c(0, 0.05))) +
  labs(
    title = "B. Post-Incident Polarization by Ideology",
    subtitle = paste0("Overall ATE = ", round(overall_ate, 3)),
    x = "Fear of Police Response",
    y = "Proportion of Responses"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    strip.text = element_text(size = 11, face = "bold")
  )

# Panel C1: Linear model pp check
panel_c1 <- pp_check(linear_model, ndraws = 50) + 
  ggtitle("C1. Linear Model PP Check") +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 9)
  )

# Panel C2: Ordinal model pp check
panel_c2 <- pp_check(ordinal_model, ndraws = 50, type = "bars") +
  ggtitle("C2. Ordinal Model PP Check") +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 9)
  )

# Panel D: Create comparison data with proper factor labels preserved
newdata_post <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  time = 1
)

linear_pred_draws <- posterior_epred(linear_model, newdata = newdata_post, re_formula = NA)

convert_linear_to_probs <- function(linear_preds, categories = 1:7) {
  n_draws <- nrow(linear_preds)
  n_conditions <- ncol(linear_preds)
  prob_array <- array(dim = c(n_draws, n_conditions, length(categories)))
  
  for(draw in 1:n_draws) {
    for(cond in 1:n_conditions) {
      mean_pred <- linear_preds[draw, cond]
      sigma <- as_draws_df(linear_model)$sigma[draw]
      
      for(cat in 1:length(categories)) {
        if(cat == 1) {
          prob_array[draw, cond, cat] <- pnorm(1.5, mean_pred, sigma)
        } else if(cat == length(categories)) {
          prob_array[draw, cond, cat] <- 1 - pnorm(cat - 0.5, mean_pred, sigma)
        } else {
          prob_array[draw, cond, cat] <- pnorm(cat + 0.5, mean_pred, sigma) - 
                                       pnorm(cat - 0.5, mean_pred, sigma)
        }
      }
    }
  }
  return(prob_array)
}

ordinal_pred_draws <- posterior_epred(ordinal_model, newdata = newdata_post, re_formula = NA)
linear_probs <- convert_linear_to_probs(linear_pred_draws)

observed_post <- sim_data %>%
  filter(time == 1) %>%
  group_by(ideology, response) %>%
  summarise(count = n(), .groups = "drop") %>%
  group_by(ideology) %>%
  mutate(proportion = count / sum(count))

# Create comparison dataset with explicit factor labels
prob_comparison_fixed <- expand_grid(
  ideology = c("Conservative", "Liberal"),
  category = 1:7
) %>%
  mutate(
    condition_idx = case_when(
      ideology == "Conservative" ~ 1,
      ideology == "Liberal" ~ 2
    ),
    linear_prob = map2_dbl(condition_idx, category, ~mean(linear_probs[, .x, .y])),
    ordinal_prob = map2_dbl(condition_idx, category, ~mean(ordinal_pred_draws[, .x, .y])),
    # Explicitly create factor with labels
    category_label = factor(category, levels = 1:7,
                       labels = c("1", "2", "3", "4", "5", "6", "7"))
    ) %>%
  left_join(
    observed_post %>% dplyr::select(ideology, response, proportion),
    by = c("ideology", "category" = "response")
  ) %>%
  pivot_longer(
    cols = c(linear_prob, ordinal_prob, proportion),
    names_to = "type", 
    values_to = "probability"
  ) %>%
  mutate(
    type_label = factor(type, 
                       levels = c("proportion", "linear_prob", "ordinal_prob"),
                       labels = c("Observed", "Linear Model", "Ordinal Model"))
  )

panel_d <- prob_comparison_fixed %>%
  ggplot(aes(x = category_label, y = probability, color = type_label, shape = type_label)) +
  geom_point(size = 2.5, position = position_dodge(width = 0.4)) +
  geom_line(aes(group = type_label), position = position_dodge(width = 0.4), alpha = 0.7) +
  facet_wrap(~ideology, labeller = label_both) +
  scale_color_manual(values = c("Observed" = "#000000", "Linear Model" = "#D55E00", "Ordinal Model" = "#0072B2"),
                    name = "") +
  scale_shape_manual(values = c("Observed" = 17, "Linear Model" = 15, "Ordinal Model" = 16),
                    name = "") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(
    title = "D. Model Predictions vs. Observed Data",
    subtitle = "Ordinal model captures bimodal reality",
    x = "Fear Response Category",
    y = "Probability"
  ) +
  theme_minimal(base_family = "serif") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 11, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, size = 10),
    axis.text = element_text(size = 9),
    axis.title = element_text(size = 10),
    legend.position = "bottom",
    panel.grid.minor = element_blank(),
    panel.grid.major.x = element_blank(),
    strip.text = element_text(size = 10, face = "bold")
  )

# Create the layout
design <- "
  AB
  AB
  CE
  DE
"

comprehensive_figure <- wrap_plots(
  A = panel_a, 
  B = panel_b,
  C = panel_c1, 
  D = panel_c2,
  E = panel_d,
  design = design,
  heights = c(1, 0.5, 1)
) + 
  plot_annotation(
    title = "Figure S1A. Ordinal Model Superiority Despite Symmetric Baseline Distributions and Assumption Violations", 
    theme = theme(
      plot.title = element_text(hjust = 0.5, size = 12, face = "bold")
      )
    )
  

comprehensive_figure