# Working with shrinkr in the Tidy Bayesian Ecosystem

#### Jacob M. Maronge

#### 2026-06-29

    knitr::opts_chunk$set(
      collapse = TRUE,
      comment = "#>",
      fig.width = 8,
      fig.height = 6,
      warning = FALSE,
      message = FALSE
    )

## Overview

shrinkr is designed to work seamlessly with modern R workflows. This vignette shows practical examples of using shrinkr with:

- tidyverse: dplyr, ggplot2, tidyr for data manipulation and visualization
- posterior: Working with MCMC draws, computing summaries and diagnostics
- bayesplot: MCMC diagnostic plots (trace plots, pairs plots, etc.)
- tidybayes: Tidy manipulation of Bayesian posteriors
- ggdist: Modern distribution visualizations

    library(shrinkr)
    library(distributional)
    library(MASS)

    # Bayesian ecosystem
    library(posterior)
    library(bayesplot)
    library(tidybayes)
    library(ggdist)

    # Tidyverse
    library(dplyr)
    library(tidyr)
    library(ggplot2)

    theme_set(theme_minimal(base_size = 12))

## Example: Multi-Region Clinical Trial

Imagine a clinical trial run across 5 regions testing a new treatment. We have Stage 1 posterior samples from region-specific analyses.

### Simulate Stage 1 Results

    set.seed(1104)

    # True effects (unknown in practice)
    true_effects <- c(0.45, 0.60, 0.38, -0.10, 0.65)
    region_names <- c("North", "South", "East", "West", "Central")

    # Simulate posterior samples from Stage 1
    samples_list <- lapply(1:5, function(i) {
      matrix(rnorm(2000, true_effects[i], 0.20), ncol = 1)
    })
    names(samples_list) <- region_names

### Fit shrinkr Model

    # Fit mixture approximation
    mix <- fit_mixture(samples_list, K_max = 3, verbose = FALSE)

    # Specify hierarchical priors
    priors <- list(
      mu = dist_normal(0, 5),
      tau = dist_truncated(dist_student_t(3, 0, 1), lower = 0)
    )

    # Run hierarchical shrinkage
    fit <- shrink(
      mixture = mix,
      hierarchical_priors = priors,
      chains = 4,
      iter = 2000,
      warmup = 1000,
      cores = 1,
      seed = 2024,
      refresh = 0
    )
    #>
    #> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 1).
    #> Chain 1:
    #> Chain 1: Gradient evaluation took 9e-06 seconds
    #> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
    #> Chain 1: Adjust your expectations accordingly!
    #> Chain 1:
    #> Chain 1:
    #> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
    #> Chain 1: Iteration:  100 / 2000 [  5%]  (Warmup)
    #> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
    #> Chain 1: Iteration:  300 / 2000 [ 15%]  (Warmup)
    #> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
    #> Chain 1: Iteration:  500 / 2000 [ 25%]  (Warmup)
    #> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
    #> Chain 1: Iteration:  700 / 2000 [ 35%]  (Warmup)
    #> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
    #> Chain 1: Iteration:  900 / 2000 [ 45%]  (Warmup)
    #> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
    #> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
    #> Chain 1: Iteration: 1100 / 2000 [ 55%]  (Sampling)
    #> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
    #> Chain 1: Iteration: 1300 / 2000 [ 65%]  (Sampling)
    #> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
    #> Chain 1: Iteration: 1500 / 2000 [ 75%]  (Sampling)
    #> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
    #> Chain 1: Iteration: 1700 / 2000 [ 85%]  (Sampling)
    #> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
    #> Chain 1: Iteration: 1900 / 2000 [ 95%]  (Sampling)
    #> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
    #> Chain 1:
    #> Chain 1:  Elapsed Time: 0.144 seconds (Warm-up)
    #> Chain 1:                0.078 seconds (Sampling)
    #> Chain 1:                0.222 seconds (Total)
    #> Chain 1:
    #>
    #> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 2).
    #> Chain 2:
    #> Chain 2: Gradient evaluation took 6e-06 seconds
    #> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
    #> Chain 2: Adjust your expectations accordingly!
    #> Chain 2:
    #> Chain 2:
    #> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
    #> Chain 2: Iteration:  100 / 2000 [  5%]  (Warmup)
    #> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
    #> Chain 2: Iteration:  300 / 2000 [ 15%]  (Warmup)
    #> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
    #> Chain 2: Iteration:  500 / 2000 [ 25%]  (Warmup)
    #> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
    #> Chain 2: Iteration:  700 / 2000 [ 35%]  (Warmup)
    #> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
    #> Chain 2: Iteration:  900 / 2000 [ 45%]  (Warmup)
    #> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
    #> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
    #> Chain 2: Iteration: 1100 / 2000 [ 55%]  (Sampling)
    #> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
    #> Chain 2: Iteration: 1300 / 2000 [ 65%]  (Sampling)
    #> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
    #> Chain 2: Iteration: 1500 / 2000 [ 75%]  (Sampling)
    #> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
    #> Chain 2: Iteration: 1700 / 2000 [ 85%]  (Sampling)
    #> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
    #> Chain 2: Iteration: 1900 / 2000 [ 95%]  (Sampling)
    #> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
    #> Chain 2:
    #> Chain 2:  Elapsed Time: 0.109 seconds (Warm-up)
    #> Chain 2:                0.085 seconds (Sampling)
    #> Chain 2:                0.194 seconds (Total)
    #> Chain 2:
    #>
    #> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 3).
    #> Chain 3:
    #> Chain 3: Gradient evaluation took 6e-06 seconds
    #> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
    #> Chain 3: Adjust your expectations accordingly!
    #> Chain 3:
    #> Chain 3:
    #> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
    #> Chain 3: Iteration:  100 / 2000 [  5%]  (Warmup)
    #> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
    #> Chain 3: Iteration:  300 / 2000 [ 15%]  (Warmup)
    #> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
    #> Chain 3: Iteration:  500 / 2000 [ 25%]  (Warmup)
    #> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
    #> Chain 3: Iteration:  700 / 2000 [ 35%]  (Warmup)
    #> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
    #> Chain 3: Iteration:  900 / 2000 [ 45%]  (Warmup)
    #> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
    #> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
    #> Chain 3: Iteration: 1100 / 2000 [ 55%]  (Sampling)
    #> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
    #> Chain 3: Iteration: 1300 / 2000 [ 65%]  (Sampling)
    #> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
    #> Chain 3: Iteration: 1500 / 2000 [ 75%]  (Sampling)
    #> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
    #> Chain 3: Iteration: 1700 / 2000 [ 85%]  (Sampling)
    #> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
    #> Chain 3: Iteration: 1900 / 2000 [ 95%]  (Sampling)
    #> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
    #> Chain 3:
    #> Chain 3:  Elapsed Time: 0.103 seconds (Warm-up)
    #> Chain 3:                0.082 seconds (Sampling)
    #> Chain 3:                0.185 seconds (Total)
    #> Chain 3:
    #>
    #> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 4).
    #> Chain 4:
    #> Chain 4: Gradient evaluation took 6e-06 seconds
    #> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
    #> Chain 4: Adjust your expectations accordingly!
    #> Chain 4:
    #> Chain 4:
    #> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
    #> Chain 4: Iteration:  100 / 2000 [  5%]  (Warmup)
    #> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
    #> Chain 4: Iteration:  300 / 2000 [ 15%]  (Warmup)
    #> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
    #> Chain 4: Iteration:  500 / 2000 [ 25%]  (Warmup)
    #> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
    #> Chain 4: Iteration:  700 / 2000 [ 35%]  (Warmup)
    #> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
    #> Chain 4: Iteration:  900 / 2000 [ 45%]  (Warmup)
    #> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
    #> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
    #> Chain 4: Iteration: 1100 / 2000 [ 55%]  (Sampling)
    #> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
    #> Chain 4: Iteration: 1300 / 2000 [ 65%]  (Sampling)
    #> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
    #> Chain 4: Iteration: 1500 / 2000 [ 75%]  (Sampling)
    #> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
    #> Chain 4: Iteration: 1700 / 2000 [ 85%]  (Sampling)
    #> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
    #> Chain 4: Iteration: 1900 / 2000 [ 95%]  (Sampling)
    #> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
    #> Chain 4:
    #> Chain 4:  Elapsed Time: 0.101 seconds (Warm-up)
    #> Chain 4:                0.081 seconds (Sampling)
    #> Chain 4:                0.182 seconds (Total)
    #> Chain 4:

## Working with posterior Package

The posterior package provides the foundation for working with MCMC draws.

### Extract Draws

    # Extract all parameters as draws_df
    draws <- as_draws_df(fit)

    # See what's available
    variables(draws)
    #> [1] "mu"          "tau"         "theta[1]"    "theta[2]"    "theta[3]"
    #> [6] "theta[4]"    "theta[5]"    "tau_squared" "lp__"

    # Extract specific parameters
    mu_tau_draws <- extract_mu_tau(fit)
    theta_draws <- extract_theta(fit)

### Basic Summaries

    # Quick summary of all parameters
    summarize_draws(draws)
    #> # A tibble: 9 × 10
    #>   variable     mean  median    sd    mad       q5    q95  rhat ess_bulk ess_tail
    #>   <chr>       <dbl>   <dbl> <dbl>  <dbl>    <dbl>  <dbl> <dbl>    <dbl>    <dbl>
    #> 1 mu         0.398   0.395  0.182 0.155   1.11e-1  0.700  1.00    1597.    1394.
    #> 2 tau        0.306   0.269  0.203 0.173   4.73e-2  0.687  1.00    1175.    1258.
    #> 3 theta[1]   0.427   0.424  0.162 0.154   1.68e-1  0.704  1.00    5061.    3283.
    #> 4 theta[2]   0.517   0.509  0.171 0.168   2.50e-1  0.811  1.00    4120.    3799.
    #> 5 theta[3]   0.378   0.383  0.162 0.154   1.09e-1  0.639  1.00    5106.    3403.
    #> 6 theta[4]   0.0997  0.108  0.211 0.219  -2.63e-1  0.430  1.00    2176.    2757.
    #> 7 theta[5]   0.549   0.539  0.179 0.183   2.80e-1  0.856  1.00    3853.    3141.
    #> 8 tau_squa…  0.135   0.0723 0.196 0.0816  2.23e-3  0.472  1.00    1175.    1258.
    #> 9 lp__      -6.44   -6.11   3.01  2.95   -1.17e+1 -2.05   1.00    1240.    2011.

    # Focus on theta parameters
    summarize_draws(theta_draws, mean, sd, median, mad, ~quantile(.x, c(0.025, 0.975)))
    #> # A tibble: 19 × 7
    #>    variable         mean    sd   median    mad     `2.5%` `97.5%`
    #>    <chr>           <dbl> <dbl>    <dbl>  <dbl>      <dbl>   <dbl>
    #>  1 mu           0.398    0.182  0.395   0.155    0.0220     0.800
    #>  2 tau          0.306    0.203  0.269   0.173    0.0255     0.806
    #>  3 theta_c[1]   0.0282   0.987  0.0260  0.993   -1.89       1.99
    #>  4 theta_c[2]   0.000676 1.02  -0.00862 1.02    -1.93       2.03
    #>  5 theta_c[3]  -0.0272   0.990 -0.0175  0.998   -1.96       1.90
    #>  6 theta_c[4]  -0.00466  0.980 -0.0140  0.971   -1.93       1.88
    #>  7 theta_c[5]   0.0108   1.04  -0.00245 1.06    -1.97       2.09
    #>  8 z[1]         0.108    0.748  0.103   0.711   -1.37       1.63
    #>  9 z[2]         0.421    0.727  0.404   0.701   -1.01       1.88
    #> 10 z[3]        -0.0531   0.711 -0.0742  0.673   -1.46       1.41
    #> 11 z[4]        -1.02     0.783 -1.02    0.770   -2.55       0.518
    #> 12 z[5]         0.524    0.753  0.508   0.730   -0.962      2.02
    #> 13 theta[1]     0.427    0.162  0.424   0.154    0.118      0.771
    #> 14 theta[2]     0.517    0.171  0.509   0.168    0.204      0.872
    #> 15 theta[3]     0.378    0.162  0.383   0.154    0.0491     0.696
    #> 16 theta[4]     0.0997   0.211  0.108   0.219   -0.334      0.473
    #> 17 theta[5]     0.549    0.179  0.539   0.183    0.230      0.918
    #> 18 tau_squared  0.135    0.196  0.0723  0.0816   0.000649   0.649
    #> 19 lp__        -6.44     3.01  -6.11    2.95   -13.1       -1.45

    # Convergence diagnostics
    summarize_draws(draws, default_convergence_measures())
    #> # A tibble: 9 × 4
    #>   variable     rhat ess_bulk ess_tail
    #>   <chr>       <dbl>    <dbl>    <dbl>
    #> 1 mu           1.00    1597.    1394.
    #> 2 tau          1.00    1175.    1258.
    #> 3 theta[1]     1.00    5061.    3283.
    #> 4 theta[2]     1.00    4120.    3799.
    #> 5 theta[3]     1.00    5106.    3403.
    #> 6 theta[4]     1.00    2176.    2757.
    #> 7 theta[5]     1.00    3853.    3141.
    #> 8 tau_squared  1.00    1175.    1258.
    #> 9 lp__         1.00    1240.    2011.

    # Custom summaries
    summarise_draws(
      theta_draws,
      mean,
      sd,
      prob_positive = ~mean(.x > 0),
      prob_large = ~mean(.x > 0.5)
    )
    #> # A tibble: 19 × 5
    #>    variable         mean    sd prob_positive prob_large
    #>    <chr>           <dbl> <dbl>         <dbl>      <dbl>
    #>  1 mu           0.398    0.182        0.980     0.247
    #>  2 tau          0.306    0.203        1         0.145
    #>  3 theta_c[1]   0.0282   0.987        0.513     0.313
    #>  4 theta_c[2]   0.000676 1.02         0.494     0.309
    #>  5 theta_c[3]  -0.0272   0.990        0.491     0.298
    #>  6 theta_c[4]  -0.00466  0.980        0.494     0.298
    #>  7 theta_c[5]   0.0108   1.04         0.498     0.323
    #>  8 z[1]         0.108    0.748        0.560     0.286
    #>  9 z[2]         0.421    0.727        0.724     0.448
    #> 10 z[3]        -0.0531   0.711        0.457     0.208
    #> 11 z[4]        -1.02     0.783        0.0855    0.026
    #> 12 z[5]         0.524    0.753        0.764     0.503
    #> 13 theta[1]     0.427    0.162        0.995     0.306
    #> 14 theta[2]     0.517    0.171        1.000     0.518
    #> 15 theta[3]     0.378    0.162        0.988     0.212
    #> 16 theta[4]     0.0997   0.211        0.688     0.0138
    #> 17 theta[5]     0.549    0.179        1.000     0.580
    #> 18 tau_squared  0.135    0.196        1         0.0435
    #> 19 lp__        -6.44     3.01         0.001     0.00025

### Check Convergence

    # Check Rhat for all parameters
    all_rhats <- summarise_draws(draws, "rhat")
    max(all_rhats$rhat, na.rm = TRUE)
    #> [1] 1.001981

    # Check effective sample size
    summarise_draws(draws, "ess_bulk", "ess_tail") %>%
      filter(ess_bulk < 400 | ess_tail < 400)
    #> # A tibble: 0 × 3
    #> # ℹ 3 variables: variable <chr>, ess_bulk <dbl>, ess_tail <dbl>

    # Detailed diagnostics for specific parameters
    summarise_draws(
      subset_draws(draws, variable = c("mu", "tau")),
      default_convergence_measures()
    )
    #> # A tibble: 2 × 4
    #>   variable  rhat ess_bulk ess_tail
    #>   <chr>    <dbl>    <dbl>    <dbl>
    #> 1 mu        1.00    1597.    1394.
    #> 2 tau       1.00    1175.    1258.

## Diagnostic Plots with bayesplot

bayesplot provides essential MCMC diagnostic visualizations.

### Trace Plots

Check for mixing and stationarity:

    # Check hyperparameters
    mcmc_trace(draws, pars = c("mu", "tau", "tau_squared"))

    # Check first few thetas
    mcmc_trace(draws, regex_pars = "theta\\[[1-3]\\]")

    # All thetas at once (if not too many)
    mcmc_trace(draws, regex_pars = "theta")

### Density Plots

Compare chains and check for multimodality:

    # Overlay densities from different chains
    mcmc_dens_overlay(draws, pars = c("mu", "tau"))

    # Individual densities
    mcmc_dens(draws, pars = c("mu", "tau", "tau_squared"))

    # Compare all thetas
    mcmc_dens_overlay(draws, regex_pars = "theta")

### Interval Plots

Visualize posterior uncertainties:

    # All thetas with 50% and 95% intervals
    mcmc_intervals(draws, regex_pars = "theta", prob = 0.5, prob_outer = 0.95)

    # With point estimates
    mcmc_intervals_data(draws, regex_pars = "theta") %>%
      ggplot(aes(y = parameter)) +
      geom_pointrange(aes(x = m, xmin = ll, xmax = hh)) +
      geom_point(aes(x = m), size = 3) +
      labs(title = "Posterior Intervals for Regional Effects", x = "Effect Size", y = NULL)

### Area Plots

Density plots with shaded intervals:

    # Hyperparameters
    mcmc_areas(draws, pars = c("mu", "tau"), prob = 0.95, prob_outer = 0.99)

    # All thetas
    mcmc_areas(draws, regex_pars = "theta", prob = 0.8)

## Tidy Analysis with tidybayes

tidybayes makes it easy to manipulate and visualize posteriors using tidy principles.

### Spread and Gather Draws

    # Gather theta parameters into long format
    theta_tidy <- draws %>%
      gather_draws(theta[region]) %>%
      mutate(region = region_names[region])

    head(theta_tidy)
    #> # A tibble: 6 × 6
    #> # Groups:   region, .variable [1]
    #>   region .chain .iteration .draw .variable .value
    #>   <chr>   <int>      <int> <int> <chr>      <dbl>
    #> 1 North       1          1     1 theta      0.596
    #> 2 North       1          2     2 theta      0.527
    #> 3 North       1          3     3 theta      0.526
    #> 4 North       1          4     4 theta      0.355
    #> 5 North       1          5     5 theta      0.467
    #> 6 North       1          6     6 theta      0.386

    # Spread into wide format
    theta_wide <- draws %>%
      spread_draws(theta[region]) %>%
      mutate(region = region_names[region])

    head(theta_wide)
    #> # A tibble: 6 × 5
    #> # Groups:   region [1]
    #>   region theta .chain .iteration .draw
    #>   <chr>  <dbl>  <int>      <int> <int>
    #> 1 North  0.596      1          1     1
    #> 2 North  0.527      1          2     2
    #> 3 North  0.526      1          3     3
    #> 4 North  0.355      1          4     4
    #> 5 North  0.467      1          5     5
    #> 6 North  0.386      1          6     6

### Point and Interval Summaries

    # Median and 95% quantile intervals
    theta_tidy %>%
      group_by(region) %>%
      median_qi(.value, .width = 0.95)
    #> # A tibble: 5 × 7
    #>   region  .value  .lower .upper .width .point .interval
    #>   <chr>    <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>
    #> 1 Central  0.539  0.230   0.918   0.95 median qi
    #> 2 East     0.383  0.0491  0.696   0.95 median qi
    #> 3 North    0.424  0.118   0.771   0.95 median qi
    #> 4 South    0.509  0.204   0.872   0.95 median qi
    #> 5 West     0.108 -0.334   0.473   0.95 median qi

    # Multiple interval widths
    theta_tidy %>%
      group_by(region) %>%
      median_qi(.value, .width = c(0.5, 0.8, 0.95))
    #> # A tibble: 15 × 7
    #>    region  .value  .lower .upper .width .point .interval
    #>    <chr>    <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>
    #>  1 Central  0.539  0.419   0.668   0.5  median qi
    #>  2 East     0.383  0.273   0.482   0.5  median qi
    #>  3 North    0.424  0.320   0.527   0.5  median qi
    #>  4 South    0.509  0.399   0.626   0.5  median qi
    #>  5 West     0.108 -0.0406  0.254   0.5  median qi
    #>  6 Central  0.539  0.332   0.788   0.8  median qi
    #>  7 East     0.383  0.174   0.574   0.8  median qi
    #>  8 North    0.424  0.227   0.630   0.8  median qi
    #>  9 South    0.509  0.304   0.744   0.8  median qi
    #> 10 West     0.108 -0.179   0.374   0.8  median qi
    #> 11 Central  0.539  0.230   0.918   0.95 median qi
    #> 12 East     0.383  0.0491  0.696   0.95 median qi
    #> 13 North    0.424  0.118   0.771   0.95 median qi
    #> 14 South    0.509  0.204   0.872   0.95 median qi
    #> 15 West     0.108 -0.334   0.473   0.95 median qi

    # Mean and HDI (highest density interval)
    theta_tidy %>%
      group_by(region) %>%
      mean_hdi(.value, .width = 0.95)
    #> # A tibble: 5 × 7
    #>   region  .value  .lower .upper .width .point .interval
    #>   <chr>    <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>
    #> 1 Central 0.549   0.198   0.881   0.95 mean   hdi
    #> 2 East    0.378   0.0620  0.706   0.95 mean   hdi
    #> 3 North   0.427   0.118   0.771   0.95 mean   hdi
    #> 4 South   0.517   0.194   0.856   0.95 mean   hdi
    #> 5 West    0.0997 -0.295   0.502   0.95 mean   hdi

### Custom Summaries with dplyr

    # Probability of positive effect
    theta_tidy %>%
      group_by(region) %>%
      summarise(
        mean_effect = mean(.value),
        sd_effect = sd(.value),
        prob_positive = mean(.value > 0),
        prob_clinically_meaningful = mean(.value > 0.3),
        .groups = "drop"
      ) %>%
      arrange(desc(prob_positive))
    #> # A tibble: 5 × 5
    #>   region  mean_effect sd_effect prob_positive prob_clinically_meaningful
    #>   <chr>         <dbl>     <dbl>         <dbl>                      <dbl>
    #> 1 Central      0.549      0.179         1.000                      0.934
    #> 2 South        0.517      0.171         1.000                      0.904
    #> 3 North        0.427      0.162         0.995                      0.787
    #> 4 East         0.378      0.162         0.988                      0.699
    #> 5 West         0.0997     0.211         0.688                      0.182

### Computing Contrasts

    # Method 1: Using shrinkr's built-in function
    L <- rbind(
      "South - North" = c(-1, 1, 0, 0, 0),
      "Central - North" = c(-1, 0, 0, 0, 1),
      "South - West" = c(0, 1, 0, -1, 0)
    )
    contrasts <- theta_contrasts(fit, L, labels = rownames(L))
    summarise_draws(contrasts)
    #> # A tibble: 3 × 10
    #>   variable        mean median    sd   mad       q5   q95  rhat ess_bulk ess_tail
    #>   <chr>          <dbl>  <dbl> <dbl> <dbl>    <dbl> <dbl> <dbl>    <dbl>    <dbl>
    #> 1 South - North 0.0904 0.0715 0.219 0.198 -0.251   0.460  1.00    5281.    3574.
    #> 2 Central - No… 0.123  0.102  0.226 0.216 -0.226   0.505  1.00    5179.    3583.
    #> 3 South - West  0.417  0.404  0.289 0.312 -0.00425 0.923  1.00    2129.    2363.

    # Method 2: Using tidybayes compare_levels
    theta_wide %>%
      compare_levels(theta, by = region, comparison = "pairwise") %>%
      group_by(region) %>%
      median_qi(theta) %>%
      arrange(desc(theta))
    #> # A tibble: 10 × 7
    #>    region            theta .lower .upper .width .point .interval
    #>    <chr>             <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
    #>  1 South - East     0.119  -0.283 0.615    0.95 median qi
    #>  2 South - North    0.0715 -0.334 0.550    0.95 median qi
    #>  3 North - East     0.0336 -0.366 0.501    0.95 median qi
    #>  4 South - Central -0.0211 -0.481 0.398    0.95 median qi
    #>  5 North - Central -0.102  -0.593 0.302    0.95 median qi
    #>  6 East - Central  -0.151  -0.670 0.241    0.95 median qi
    #>  7 West - East     -0.257  -0.797 0.135    0.95 median qi
    #>  8 West - North    -0.306  -0.907 0.0952   0.95 median qi
    #>  9 West - South    -0.404  -1.03  0.0373   0.95 median qi
    #> 10 West - Central  -0.437  -1.07  0.0277   0.95 median qi

## Modern Visualizations with ggdist

ggdist provides publication-ready distribution visualizations.

### Halfeye Plots

Eye + interval visualization:

    theta_tidy %>%
      ggplot(aes(y = region, x = .value)) +
      stat_halfeye(
        .width = c(0.66, 0.95),
        fill = "steelblue"
      ) +
      geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
      labs(
        title = "Regional Treatment Effects",
        subtitle = "Posterior distributions with median and 66%/95% intervals",
        x = "Treatment Effect",
        y = NULL
      )

### Slab + Interval

Density with separate interval layer:

    theta_tidy %>%
      ggplot(aes(y = region, x = .value)) +
      stat_slab(aes(fill_ramp = after_stat(level)), fill = "steelblue", alpha = 0.8) +
      stat_pointinterval(.width = c(0.66, 0.95), position = position_nudge(y = -0.15)) +
      scale_fill_ramp_discrete(range = c(1, 0.2), guide = "none") +
      labs(
        title = "Posterior Densities with Quantile Intervals",
        x = "Treatment Effect",
        y = NULL
      )

### Quantile Dotplots

Each dot = quantile of the distribution:

    theta_tidy %>%
      ggplot(aes(y = region, x = .value)) +
      stat_dots(quantiles = 100) +
      geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
      labs(
        title = "Quantile Dotplots",
        subtitle = "Each dot represents 1% of the posterior",
        x = "Treatment Effect",
        y = NULL
      )

### Gradient Intervals

Continuous representation of uncertainty:

    theta_tidy %>%
      ggplot(aes(y = region, x = .value)) +
      stat_gradientinterval(.width = ppoints(50)) +
      scale_color_brewer(palette = "Blues", guide = "none") +
      labs(
        title = "Gradient Interval Representation",
        x = "Treatment Effect",
        y = NULL
      )

## Comparing Pre- and Post-Shrinkage

### Extract Both Estimates

    # Get pre-shrunk estimates from mixture
    pre_shrunk <- summarise_theta(fit) %>%
      mutate(type = "Pre-shrunk")

    # Get post-shrunk estimates
    post_shrunk <- summarise_theta(fit) %>%
      mutate(type = "Post-shrunk")

    # Or use shrinkr's built-in plot
    plot(fit, group_names = region_names)

### Custom Comparison Plot

    # Get the hierarchical mean (mu)
    mu_draws <- draws %>% spread_draws(mu)
    mu_mean <- mean(mu_draws$mu)

    # Combine with Stage 1 samples
    stage1_draws <- lapply(seq_along(samples_list), function(i) {
      data.frame(
        region = region_names[i],
        .value = samples_list[[i]][,1],
        type = "Stage 1"
      )
    }) %>% bind_rows()

    stage2_draws <- theta_tidy %>%
      mutate(type = "Stage 2 (Shrunk)")

    # Plot side by side
    bind_rows(stage1_draws, stage2_draws) %>%
      ggplot(aes(y = region, x = .value, fill = type)) +
      stat_halfeye(alpha = 0.7, position = position_dodge(width = 0.4)) +
      geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5, color = "gray50") +
      geom_vline(xintercept = mu_mean, linetype = "solid", alpha = 0.8,
                 color = "darkred", linewidth = 1) +
      annotate("text", x = mu_mean, y = 0.5,
               label = sprintf("Global mean (μ) = %.2f", mu_mean),
               hjust = -0.1, color = "darkred", size = 3.5) +
      scale_fill_manual(values = c("Stage 1" = "gray70", "Stage 2 (Shrunk)" = "steelblue")) +
      labs(
        title = "Stage 1 vs Stage 2: Effect of Hierarchical Shrinkage",
        subtitle = "Stage 2 estimates are pulled toward the global mean",
        x = "Treatment Effect",
        y = NULL,
        fill = NULL
      ) +
      theme(legend.position = "bottom")

## Complete Workflow Example

Here’s a typical analysis workflow using tidy principles:

    # 1. Extract and prepare data
    analysis_data <- draws %>%
      spread_draws(mu, tau, theta[i]) %>%
      mutate(region = region_names[i])

    # 2. Compute summaries
    summary_table <- analysis_data %>%
      group_by(region) %>%
      summarise(
        mean = mean(theta),
        median = median(theta),
        sd = sd(theta),
        q025 = quantile(theta, 0.025),
        q975 = quantile(theta, 0.975),
        prob_positive = mean(theta > 0),
        prob_clinically_important = mean(theta > 0.3),
        .groups = "drop"
      ) %>%
      arrange(desc(median))

    print(summary_table)
    #> # A tibble: 5 × 8
    #>   region    mean median    sd    q025  q975 prob_positive prob_clinically_impo…¹
    #>   <chr>    <dbl>  <dbl> <dbl>   <dbl> <dbl>         <dbl>                  <dbl>
    #> 1 Central 0.549   0.539 0.179  0.230  0.918         1.000                  0.934
    #> 2 South   0.517   0.509 0.171  0.204  0.872         1.000                  0.904
    #> 3 North   0.427   0.424 0.162  0.118  0.771         0.995                  0.787
    #> 4 East    0.378   0.383 0.162  0.0491 0.696         0.988                  0.699
    #> 5 West    0.0997  0.108 0.211 -0.334  0.473         0.688                  0.182
    #> # ℹ abbreviated name: ¹​prob_clinically_important

    # 3. Create advanced figure
    library(patchwork)

    p1 <- analysis_data %>%
      ggplot(aes(y = reorder(region, theta), x = theta)) +
      stat_halfeye(.width = c(0.66, 0.95), fill = "steelblue") +
      geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
      labs(
        title = "A. Regional Treatment Effects",
        x = "Effect Size",
        y = NULL
      )

    p2 <- analysis_data %>%
      dplyr::ungroup() %>%
      dplyr::select(mu, tau, .draw) %>%
      dplyr::distinct() %>%
      tidyr::pivot_longer(cols = c(mu, tau), names_to = "name", values_to = "value") %>%
      ggplot(aes(x = value, fill = name)) +
      stat_halfeye(alpha = 0.7) +
      facet_wrap(~name, scales = "free", labeller = label_both) +
      scale_fill_brewer(palette = "Set2") +
      labs(
        title = "B. Hyperparameters",
        x = "Value",
        y = "Density"
      ) +
      theme(legend.position = "none")

    p3 <- analysis_data %>%
      dplyr::ungroup() %>%
      dplyr::select(.draw, region, theta) %>%
      compare_levels(theta, by = region) %>%
      ggplot(aes(y = region, x = theta)) +
      stat_halfeye(fill = "coral", alpha = 0.7) +
      geom_vline(xintercept = 0, linetype = "dashed", color = "red", alpha = 0.5) +
      labs(
        title = "C. Pairwise Regional Comparisons",
        x = "Difference in Effect Size",
        y = NULL
      )

    p4 <- analysis_data %>%
      dplyr::ungroup() %>%
      dplyr::select(.draw, mu, tau) %>%
      dplyr::distinct() %>%
      ggplot(aes(x = mu, y = tau)) +
      geom_hex(bins = 30) +
      stat_ellipse(level = 0.95, color = "red", linewidth = 1) +
      scale_fill_viridis_c() +
      labs(
        title = "D. Hyperparameter Correlation",
        x = expression(mu~"(global mean)"),
        y = expression(tau~"(heterogeneity)")
      )

    (p1 + p2) / (p3 + p4) +
      plot_annotation(
        title = "Complete Bayesian Shrinkage Analysis",
        subtitle = sprintf(
          "Global effect: %.2f [%.2f, %.2f] | Heterogeneity (tau): %.2f",
          median(analysis_data$mu),
          quantile(analysis_data$mu, 0.025),
          quantile(analysis_data$mu, 0.975),
          median(analysis_data$tau)
        )
      )

## Advanced: Custom Analyses

### Probability Statements

    # Which region is best?
    analysis_data %>%
      group_by(.draw) %>%
      slice_max(theta, n = 1) %>%
      ungroup() %>%
      count(region) %>%
      mutate(probability = n / sum(n)) %>%
      arrange(desc(probability))
    #> # A tibble: 5 × 3
    #>   region      n probability
    #>   <chr>   <int>       <dbl>
    #> 1 Central  1679      0.420
    #> 2 South    1260      0.315
    #> 3 North     625      0.156
    #> 4 East      395      0.0988
    #> 5 West       41      0.0102

    # Alternative: probability each region is best
    analysis_data %>%
      group_by(.draw) %>%
      mutate(rank = rank(-theta)) %>%
      ungroup() %>%
      group_by(region) %>%
      summarise(
        prob_best = mean(rank == 1),
        prob_top2 = mean(rank <= 2),
        mean_rank = mean(rank),
        .groups = "drop"
      ) %>%
      arrange(mean_rank)
    #> # A tibble: 5 × 4
    #>   region  prob_best prob_top2 mean_rank
    #>   <chr>       <dbl>     <dbl>     <dbl>
    #> 1 Central    0.420     0.709       2.02
    #> 2 South      0.315     0.630       2.23
    #> 3 North      0.156     0.372       2.87
    #> 4 East       0.0988    0.254       3.20
    #> 5 West       0.0102    0.0357      4.67

    # Pairwise comparisons: Probability that South > North
    # Create wide format for comparisons
    theta_wide_for_contrasts <- analysis_data %>%
      ungroup() %>%
      dplyr::select(.draw, region, theta) %>%
      tidyr::pivot_wider(names_from = region, values_from = theta)

    theta_wide_for_contrasts %>%
      summarise(
        prob_south_beats_north = mean(South > North),
        prob_south_beats_north_by_02 = mean((South - North) > 0.2),
        prob_central_beats_all = mean(
          Central > North & Central > South &
          Central > East & Central > West
        )
      )
    #> # A tibble: 1 × 3
    #>   prob_south_beats_north prob_south_beats_north_by_02 prob_central_beats_all
    #>                    <dbl>                        <dbl>                  <dbl>
    #> 1                  0.657                        0.289                  0.420

### Tail Probabilities

    # Classify effects into categories
    theta_tidy %>%
      group_by(region) %>%
      summarise(
        prob_harm = mean(.value < -0.1),
        prob_null = mean(abs(.value) < 0.1),
        prob_small_benefit = mean(.value > 0.1 & .value < 0.3),
        prob_large_benefit = mean(.value > 0.3),
        .groups = "drop"
      ) %>%
      arrange(desc(prob_large_benefit))
    #> # A tibble: 5 × 5
    #>   region  prob_harm prob_null prob_small_benefit prob_large_benefit
    #>   <chr>       <dbl>     <dbl>              <dbl>              <dbl>
    #> 1 Central   0         0.00325             0.0632              0.934
    #> 2 South     0         0.003               0.0932              0.904
    #> 3 North     0.001     0.019               0.193               0.787
    #> 4 East      0.00325   0.043               0.254               0.699
    #> 5 West      0.170     0.316               0.331               0.182

    # Visualize classification
    theta_tidy %>%
      mutate(
        category = case_when(
          .value < -0.1 ~ "Harm",
          abs(.value) < 0.1 ~ "Null",
          .value > 0.1 & .value < 0.3 ~ "Small Benefit",
          .value > 0.3 ~ "Large Benefit"
        )
      ) %>%
      count(region, category) %>%
      group_by(region) %>%
      mutate(probability = n / sum(n)) %>%
      ggplot(aes(x = probability, y = region, fill = category)) +
      geom_col(position = "stack") +
      scale_fill_manual(
        values = c(
          "Harm" = "red",
          "Null" = "gray",
          "Small Benefit" = "lightblue",
          "Large Benefit" = "darkblue"
        )
      ) +
      labs(
        title = "Classification of Treatment Effects",
        x = "Probability",
        y = NULL,
        fill = "Effect Category"
      ) +
      theme(legend.position = "bottom")

### Ranking Analysis

    # Compute ranks for each draw
    rank_data <- analysis_data %>%
      group_by(.draw) %>%
      mutate(rank = rank(-theta)) %>%
      ungroup()

    # Summary statistics
    rank_summary <- rank_data %>%
      group_by(region) %>%
      summarise(
        mean_rank = mean(rank),
        median_rank = median(rank),
        prob_rank1 = mean(rank == 1),
        prob_rank2 = mean(rank == 2),
        prob_top3 = mean(rank <= 3),
        .groups = "drop"
      ) %>%
      arrange(mean_rank)

    print(rank_summary)
    #> # A tibble: 5 × 6
    #>   region  mean_rank median_rank prob_rank1 prob_rank2 prob_top3
    #>   <chr>       <dbl>       <dbl>      <dbl>      <dbl>     <dbl>
    #> 1 Central      2.02           2     0.420      0.289     0.874
    #> 2 South        2.23           2     0.315      0.314     0.854
    #> 3 North        2.87           3     0.156      0.216     0.662
    #> 4 East         3.20           3     0.0988     0.155     0.53
    #> 5 West         4.67           5     0.0102     0.0255    0.0788

    # Visualize ranking distribution
    rank_data %>%
      ggplot(aes(x = rank, y = reorder(region, -theta))) +
      stat_dots(quantiles = 100) +
      scale_x_continuous(breaks = 1:5) +
      labs(
        title = "Ranking Distribution",
        subtitle = "Each dot represents 1% of posterior draws",
        x = "Rank (1 = best, 5 = worst)",
        y = NULL
      )

    # Alternative: bar chart of ranking probabilities
    rank_data %>%
      count(region, rank) %>%
      group_by(region) %>%
      mutate(probability = n / sum(n)) %>%
      ggplot(aes(x = rank, y = probability, fill = region)) +
      geom_col() +
      facet_wrap(~region, ncol = 1) +
      scale_x_continuous(breaks = 1:5) +
      scale_fill_brewer(palette = "Set2") +
      labs(
        title = "Probability of Each Rank by Region",
        x = "Rank (1 = best)",
        y = "Probability"
      ) +
      theme(legend.position = "none")

## Further Reading

- posterior package: https://mc-stan.org/posterior/
- bayesplot package: https://mc-stan.org/bayesplot/
- tidybayes package: http://mjskay.github.io/tidybayes/
- ggdist package: https://mjskay.github.io/ggdist/