fuzzy-stacked-ologit-frequentist.md

library(MASS)
library(tidyverse)
library(marginaleffects)
library(palmerpenguins)

# Make a categorical weight column
penguins <- penguins |> 
  drop_na(sex) |> 
  mutate(weight_cat = cut(
    body_mass_g, breaks = 3, ordered_result = TRUE, 
    labels = c("Light", "Medium", "Heavy"))
  )

# Predict weight based on sex
mod <- polr(weight_cat ~ sex, data = penguins, Hess = TRUE)

# Predictions
preds_exact <- predictions(mod, by = "sex")
preds_exact
#> 
#>   Group    sex Estimate Std. Error     z Pr(>|z|)     S  2.5 % 97.5 %
#>  Light  male     0.2348     0.0314  7.47   <0.001  43.5 0.1732  0.296
#>  Light  female   0.6517     0.0360 18.11   <0.001 241.1 0.5812  0.722
#>  Medium male     0.4868     0.0327 14.89   <0.001 164.2 0.4227  0.551
#>  Medium female   0.2888     0.0296  9.77   <0.001  72.4 0.2308  0.347
#>  Heavy  male     0.2784     0.0337  8.26   <0.001  52.6 0.2124  0.344
#>  Heavy  female   0.0595     0.0125  4.77   <0.001  19.0 0.0351  0.084
#> 
#> Columns: group, sex, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high 
#> Type:  probs

# Predictions through bootstrapping
preds_bootstrapped <- predictions(mod, by = "sex") |> 
  inferences(method = "rsample", R = 500)

preds_bootstrapped
#> 
#>   Group    sex Estimate  2.5 % 97.5 %
#>  Light  male     0.2348 0.1717 0.2995
#>  Light  female   0.6517 0.5781 0.7135
#>  Medium male     0.4868 0.4279 0.5502
#>  Medium female   0.2888 0.2351 0.3482
#>  Heavy  male     0.2784 0.2135 0.3558
#>  Heavy  female   0.0595 0.0397 0.0824
#> 
#> Columns: group, sex, estimate, conf.low, conf.high 
#> Type:  probs

# These aren't actually posterior draws, but posteriordraws() can extract the bootstrapped samples 
preds_long <- posteriordraws(preds_bootstrapped)

# Extract labels from exact predictions
preds_labs <- preds_exact |>
  arrange(desc(group)) |> 
  group_by(sex) |> 
  mutate(
    y_end = cumsum(estimate),
    y_start = lag(y_end, default = 0),
    y = y_end - ((y_end - y_start) / 2),
    y_lab = glue::glue("{round(estimate, 2)} ({round(std.error, 2)})")
  )

# Plot
ggplot(preds_long, aes(x = as.numeric(drawid), y = draw)) + 
  geom_area(aes(fill = group), position = position_stack()) + 
  geom_label(
    data = preds_labs, aes(x = 250, y = y, label = y_lab), 
    fill = scales::alpha("white", 0.4), label.size = 0, fontface = "bold"
  ) +
  labs(y = "Cumulative predicted probability", fill = "Weight category") +
  facet_wrap(vars(sex))

^{Created on 2024-06-21 with reprex v2.1.0}