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