Programmatic workflow demonstration

In shinymrp: Interface for Multilevel Regression and Poststratification

This vignette provides an in-depth walk-through of the MRPWorkflow and MRPModel classes in the shinymrp package. You will find practical information about the purpose, arguments, and caveats of each method, enabling you to conduct MRP analyses programmatically.

For an introduction to the package and workflow concepts, begin with the Getting started with shinymrp vignette.

run_stan <- requireNamespace("cmdstanr", quietly = TRUE)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

knitr::opts_template$set(
  tall_plot = list(
    fig.width = 10,
    fig.height = 6,
    fig.align = "center",
    out.width = "90%"
  ),
  long_plot = list(
    fig.width = 13,
    fig.height = 6,
    fig.align = "center",
    out.width = "90%"
  ),
  map = list(
    fig.height = 4,
    out.width = "100%"
  )
)

sample_data <- readr::read_csv(
  "data/timevarying_binomial_prep.csv",
  show_col_types = FALSE
)
workflow <- qs2::qs_read("data/workflow/workflow.qs2")
model1 <- qs2::qs_read("data/workflow/model.qs2")
compare_df <- readr::read_csv(
  "data/workflow/loo.csv",
  show_col_types = FALSE
)

Initializing the workflow

MRP analysis begins by instantiating an MRPWorkflow object with mrp_workflow(). The object-oriented design efficiently manages your data and supports a reproducible workflow without requiring repeated user input. All workflow methods are accessed via the R6 $ operator.

library(shinymrp)
workflow <- mrp_workflow()

1. Data preparation

MRP requires two key data components:

• Sample data: Your analysis sample, which must include the outcome of interest (binary or continuous) and predictors (demographic, geographic, etc.), such as COVID-19 test records or survey responses.
• Poststratification data: A table containing population sizes of groups defined by demographic and geographic factors.

For details on accepted formats and preprocessing, see the Data Preprocessing vignette.

1.1 Preprocessing the sample data

Sample data should be a tidy data frame with your outcome measure and predictors. Time-varying data are also supported, allowing dates or indices (e.g., week, month, or year) to be incorporated for time-varying summaries and estimates.

The example_sample_data() function demonstrates accepted structures and common options:

• whether a time variable is included
• whether the data are aggregated
• whether the data follow a specific standard (e.g., medical records from the Epic system)
• whether the outcome measure is binary (modeled with binomial distributions) or continuous (modeled with Gaussian/normal distributions)

sample_data <- example_sample_data(
  is_timevar = TRUE,
  is_aggregated = TRUE,
  special_case = NULL,
  family = "binomial"
)
head(sample_data)

head(sample_data)

Data preprocessing is performed with the $preprocess() method, which standardizes formats (e.g., recodes age groups, creates time indices), removes duplicates, and prepares data for modeling.

Note: Only basic cleaning is performed. Conduct preliminary data checks beforehand.

workflow$preprocess(
  sample_data,
  is_timevar = TRUE,
  is_aggregated = TRUE,
  special_case = NULL,
  family = "binomial"
)

1.2 Assembling poststratification data

1.2.1 Linking to ACS

Use the $link_acs() method to retrieve poststratification data from the American Community Survey (ACS), specifying the desired geographic level and reference year (i.e., the last year of the five-year ACS data collection). Supported levels include ZIP code, county, and state. If not specified, poststratification data are aggregated for the U.S. population.

workflow$link_acs(link_geo = "zip", acs_year = 2021)

1.2.2 Importing custom poststratification data

You can also import your own poststratification data using $load_pstrat(). Requirements mirror those for sample data.

pstrat_data <- example_pstrat_data()
workflow$load_pstrat(pstrat_data, is_aggregated = TRUE)

2. Descriptive statistics

The MRPWorkflow class provides methods for creating visualizations of demographic distributions, geographic coverage, and outcome summaries.

Use $demo_bars() to generate bar plots comparing demographic distributions in the sample and poststratification data.

workflow$demo_bars(demo = "sex")

For analyzing COVID-19 test data with linked geographic predictors, the $covar_hist() method creates frequency histograms of the geographic predictors. This method will give the following error if it is called with incompatible data that do not include geographic predictors.

workflow$covar_hist(covar = "income")

$sample_size_map() visualizes sample sizes across geographic areas in an interactive choropleth map, powered by highcharter.

workflow$sample_size_map()

$outcome_plot() creates plots of average outcome values, including time variation if present.

workflow$outcome_plot()

$outcome_map() visualizes outcome summaries by geography, with options for the maximum or minimum value across time.

workflow$outcome_map(summary_type = "max")

3. Model building

The MRPModel class is a wrapper around CmdStanR objects, providing a high-level interface for obtaining posterior summaries and diagnostics.

3.1 Model specification

Use $create_model() to specify the model, including the mean structure (fixed and varying effects, two-way interactions) and priors. Varying effects assume normal distributions with unknown standard deviations (with specified priors). Stan code is generated for compilation via CmdStanR. Supported prior types include:

• normal(mu, sigma)
• student_t(nu, mu, sigma)
• structured (see Si et al., 2020)
• icar (see Besag et al., 1991)
• bym2 (see Riebler et al., 2016)

Default priors are used if "" (empty string) is supplied:

• Overall intercept: normal(0, 5)
• Coefficient: normal(0, 3)
• Standard deviation (main effect): normal₊(0, 3)
• Standard deviation (interaction): normal₊(0, 1)

model1 <- workflow$create_model(
  intercept_prior = "normal(0, 4)",
  fixed = list(
    sex = "normal(0, 2)",
    race = "normal(0, 2)"
  ),
  varying = list(
    age = "",
    time = ""
  )
)

3.2 Model fitting

The $fit() method runs Stan's main Markov chain Monte Carlo algorithm via CmdStanR.

model1$fit(
  n_iter = 500,
  n_chains = 2,
  seed = 123,
  show_messages = FALSE,
  show_exceptions = FALSE
)

3.3 Posterior summaries and diagnostics

The fitted MRPModel object extracts draws from the internal CmdStanMCMC object for summaries and diagnostics.

• Posterior summary and convergence diagnostics: Use $summary() for posterior statistics, effective sample sizes, and R-hat diagnostics.

model1$summary()

• Additional diagnostics: Use $diagnostics() for detailed diagnostic information about the MCMC process.

model1$diagnostics()

• Posterior predictive check: Use $pp_check() to compare posterior predictive samples against observed data.

workflow$pp_check(model1)

• Model comparison: Compare models via leave-one-out cross-validation with $compare_models().

model2 <- workflow$create_model(
  intercept_prior = "normal(0, 4)",
  fixed = list(
    sex = "normal(0, 2)",
    race = "normal(0, 2)"
  ),
  varying = list(
    age = "",
    time = ""
  ),
  interaction = list(
    `age:time` = "normal(0, 1)",
    `race:time` = "normal(0, 1)",
    `sex:time` = "normal(0, 1)"
  )
)

model2$fit(
  n_iter = 500,
  n_chains = 2,
  seed = 123,
  show_messages = FALSE,
  show_exceptions = FALSE
)

workflow$compare_models(model1, model2)

compare_df

3.4 Saving and loading models

You can easily save and restore models using the qs2 package. The $save() gathers all posterior draws and diagnostics, then calls qs2::qs_save() internally.

model1$save(file = "model1.qs2")

# load back into workspace
model1 <- qs2::qs_read("model1.qs2")

4. Result visualization

Use $estimate_plot() to visualize the model estimates, either overall or by group, with a 95% credible interval by default. Users can specify alternative intervals.

workflow$estimate_plot(model1, group = "sex", interval = 0.95)

For time-varying data, $estimate_map() creates snapshots of geography-based estimates at specific time points. time_index is irrelevant for cross-sectional data.

workflow$estimate_map(model1, geo = "county", interval = 0.95, time_index = 1)

Further reading

• Getting Started with shinymrp
• Data Preprocessing
• Stan Model Diagnostics

shinymrp documentation built on Dec. 4, 2025, 5:07 p.m.