Last active
November 4, 2025 18:10
-
-
Save carlislerainey/a8bc0c725e9b60dc1a7a6a7a24cce0eb to your computer and use it in GitHub Desktop.
MRP Example
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # from: https://marginaleffects.com/chapters/mrp.html | |
| # load packages | |
| library(tidyverse) | |
| library(brms) | |
| library(marginaleffects) | |
| library(ggdist) | |
| library(maps) # for map_data("state") | |
| library(viridisLite) # palettes (used via scale_fill_viridis_c) | |
| # set seed for reproducibility | |
| set.seed(123) | |
| # ces survey subset used throughout | |
| survey <- get_dataset("ces_survey") |> | |
| glimpse() | |
| # priors | |
| priors_vague <- c( | |
| prior(normal(0, 1e6), class = "b"), | |
| prior(normal(0, 1e6), class = "Intercept") | |
| ) | |
| # fit model | |
| model <- brm( | |
| defund ~ age + education + military + gender + (1 + gender | state), | |
| family = bernoulli, | |
| prior = priors_vague, | |
| data = survey, | |
| chains = 4, | |
| cores = 4 | |
| ) | |
| # poststratification frame | |
| demographics <- get_dataset("ces_demographics") |> | |
| glimpse() | |
| # predictions for each row in the poststratification frame | |
| p_ps <- predictions(model, newdata = demographics) |> | |
| glimpse() | |
| # weighted state-level average predictions (MRP estimates) | |
| p_state <- avg_predictions( | |
| model, | |
| newdata = demographics, | |
| wts = "percent", | |
| by = "state" | |
| ) |> | |
| glimpse() | |
| # plot posterior density of state estimates | |
| p_state_draws <- p_state |> | |
| get_draws(shape = "rvar") |> | |
| sort_by(~ estimate) |> | |
| glimpse() | |
| p_state_draws$state <- factor(p_state_draws$state, levels = p_state_draws$state) | |
| ggplot(p_state_draws, aes(y = state, xdist = rvar)) + | |
| stat_slab(height = 2, color = "white") + | |
| labs(x = "Posterior density of Support for Budget Cuts to Police", y = NULL) + | |
| xlim(c(.3, .5)) | |
| # plot map | |
| us_map <- map_data("state") |> | |
| as_tibble() |> | |
| rename(region = region) |> | |
| glimpse() | |
| state_lookup <- tibble( | |
| state = state.abb, | |
| region = tolower(state.name) | |
| ) |> | |
| glimpse() | |
| plot_df <- p_state_draws %>% | |
| inner_join(state_lookup, by = "state") %>% | |
| inner_join(us_map, by = "region") |> | |
| glimpse() | |
| ggplot(plot_df, aes(x = long, y = lat, group = group, fill = estimate)) + | |
| geom_polygon(color = "white", linewidth = 0.2) + | |
| coord_map("albers", lat0 = 39, lat1 = 45) + | |
| scale_fill_viridis_c(name = "Support for Budget Cuts to Police") + | |
| theme_void(base_size = 12) + | |
| theme( | |
| legend.position = "right", | |
| plot.title = element_text(face = "bold") | |
| ) | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment