Getting Started with shrinkr

What is shrinkr?

shrinkr lets you apply Bayesian hierarchical shrinkage to group-specific estimates in a modular, two-stage workflow:

Stage 1: Fit your model separately for each group (region, hospital, study, etc.) using flat priors
Stage 2: Borrow strength across groups through hierarchical shrinkage

Why use shrinkr?

You’ve already fit separate models and want to pool information afterward
You want modularity to separate model fitting from shrinkage estimation
You want to simplify shrinkage estimation in complex models
You want to try different shrinkage priors without refitting expensive Stage 1 models
You need transparency in how much your results depend on the hierarchical prior
You’re doing meta-analysis or federated learning with summary statistics

library(shrinkr)
library(distributional)
library(posterior)
library(tidyverse)

A Complete Example: Regional Clinical Trial

Imagine a clinical trial run independently in 5 regions. Each region estimated a treatment effect, and now you want to apply shrinkage analysis.

Stage 1: Fit Independent Models with Stan

First, let’s create synthetic trial data:

library(rstan)
set.seed(1104)

true_mu <- 0.5
true_tau <- 0.3
true_effects <- c(0.45, 0.72, 0.38, 0.55, 0.61)

regions <- c("North", "South", "East", "West", "Central")
n_per_region <- c(100, 80, 120, 90, 70)

trial_data <- lapply(seq_along(regions), function(i) {
  n <- n_per_region[i]
  data.frame(
    region = regions[i],
    treatment = rep(c(0, 1), each = n/2),
    outcome = c(
      rnorm(n/2, mean = 0, sd = 1),
      rnorm(n/2, mean = true_effects[i], sd = 1)
    )
  )
}) %>% bind_rows()

head(trial_data)
table(trial_data$region, trial_data$treatment)

#>   region treatment     outcome
#> 1  North         0 -0.04382439
#> 2  North         0  0.64111730
#> 3  North         0 -0.33395868
#> 4  North         0 -2.60279243
#> 5  North         0  1.17838867
#> 6  North         0 -0.20705477
#>          
#>            0  1
#>   Central 35 35
#>   East    60 60
#>   North   50 50
#>   South   40 40
#>   West    45 45

Now we write a Stan model with treatment-by-region interaction:

stan_code <- "
data {
  int<lower=0> N;
  int<lower=1> G;
  vector[N] y;
  vector[N] treatment;
  array[N] int<lower=1,upper=G> region;
}
parameters {
  vector[G] beta_region;
  real<lower=0> sigma;
}
model {
  // IMPORTANT: Flat prior on beta_region - critical for the two-stage approach!
  sigma ~ normal(0, 2);
  for (n in 1:N) {
    y[n] ~ normal(treatment[n] * beta_region[region[n]], sigma);
  }
}
"

regions <- c("North", "South", "East", "West", "Central")
region_indices <- as.integer(factor(trial_data$region, levels = regions))

fit_stan <- stan(
  model_code = stan_code,
  data = list(
    N = nrow(trial_data), G = length(regions),
    y = trial_data$outcome, treatment = trial_data$treatment,
    region = region_indices
  ),
  chains = 4, iter = 2000, warmup = 1000, refresh = 0, seed = 123
)

beta_draws <- rstan::extract(fit_stan, pars = "beta_region")$beta_region
samples_list <- lapply(seq_along(regions), function(i) beta_draws[, i])
names(samples_list) <- regions
samples <- lapply(samples_list, function(x) matrix(x, ncol = 1))

Let’s examine what we got from Stage 1:

stage1_summary <- data.frame(
  region = regions,
  mean   = sapply(samples, mean),
  sd     = sapply(samples, sd),
  lower  = sapply(samples, function(x) quantile(x, 0.025)),
  upper  = sapply(samples, function(x) quantile(x, 0.975))
)

print(stage1_summary)

#>          region      mean        sd      lower     upper
#> North     North 0.6040241 0.1405588 0.33617456 0.8799941
#> South     South 0.6000233 0.1540522 0.29866430 0.8924705
#> East       East 0.6152066 0.1265935 0.36515463 0.8666375
#> West       West 0.7049122 0.1494923 0.41466597 0.9940618
#> Central Central 0.3626256 0.1676080 0.02935275 0.6966270

Stage 1 visualization:

ggplot(stage1_summary, aes(x = region, y = mean)) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(ymin = lower, ymax = upper),
                  size = 0.8, color = "steelblue") +
  labs(
    title    = "Stage 1: Independent Regional Estimates",
    subtitle = "Each region analyzed separately with flat priors",
    x = "Region", y = "Treatment Effect",
    caption = "Points show posterior means; bars show 95% credible intervals"
  ) +
  theme_minimal(base_size = 12)

Notice: Central has the widest interval (n=70) and East the narrowest (n=120). Estimates vary considerably — hierarchical shrinkage will borrow strength across regions.

Stage 2: Apply Hierarchical Shrinkage

Step 1: Fit Mixture Approximation

mix <- fit_mixture(samples = samples, K_max = 3, verbose = TRUE)

print(mix)

Check the approximation quality:

# Blue line should overlay the histogram well
plot(mix, draws = samples, type = "density")


# Points should fall near the diagonal
plot(mix, draws = samples, type = "qq")

Step 2: Specify Hierarchical Priors

hierarchical_priors <- list(
  mu  = dist_normal(0, 1),
  tau = dist_truncated(dist_student_t(3, 0, 0.5), lower = 0)
)

Check prior implications before fitting:

prior_pred <- sample_prior_predictive(
  hierarchical_priors = hierarchical_priors,
  n_groups = 5,
  n_draws  = 1000
)

cat("Prior on tau (between-region SD):\n")
cat("  Median:", round(median(prior_pred$tau), 2), "\n")
cat("  95% interval:", round(quantile(prior_pred$tau, c(0.025, 0.975)), 2), "\n\n")

cat("Implied variation in regional effects:\n")
cat("  Typical range:", round(median(prior_pred$implied_range), 2), "\n")
cat("  95% interval:", round(quantile(prior_pred$implied_range, c(0.025, 0.975)), 2), "\n")

#> Prior on tau (between-region SD):
#>   Median: 0.38
#>   95% interval: 0.02 1.95
#> Implied variation in regional effects:
#>   Typical range: 0.81
#>   95% interval: 0.03 4.83

plot(prior_pred)

Check what the prior implies about pairwise subgroup differences:

The implied_range above measures max(theta) - min(theta) across all groups for each draw. For a more detailed view, prior_pairwise_differences() computes the distribution of |theta_i - theta_j| for every pair of groups. This is particularly useful for calibrating whether your prior places reasonable probability on clinically meaningful differences.

pw <- prior_pairwise_differences(prior_pred)
print(pw)
#> == Prior Predictive: Pairwise |theta_i - theta_j| ==
#> 
#> Groups:  5 
#> Pairs:   10 
#> Draws:   1000 
#> 
#> Overall quantiles of |theta_i - theta_j|:
#>    q2.5 = 0.005, q25 = 0.094, q50 = 0.302, q75 = 0.752, q97.5 = 2.958 
#> 
#> Per-pair summary:
#> # A tibble: 10 × 6
#>    pair             median    q2.5 q97.5 prob_gt_0.5 prob_gt_1
#>    <chr>             <dbl>   <dbl> <dbl>       <dbl>     <dbl>
#>  1 group1 vs group2  0.314 0.00517  2.94       0.364     0.185
#>  2 group1 vs group3  0.320 0.00507  2.93       0.37      0.182
#>  3 group1 vs group4  0.318 0.00522  3.02       0.362     0.193
#>  4 group1 vs group5  0.298 0.00463  2.95       0.376     0.177
#>  5 group2 vs group3  0.316 0.00555  2.81       0.361     0.185
#>  6 group2 vs group4  0.297 0.00444  2.71       0.35      0.172
#>  7 group2 vs group5  0.297 0.00569  3.07       0.36      0.189
#>  8 group3 vs group4  0.306 0.00333  2.92       0.368     0.167
#>  9 group3 vs group5  0.294 0.00443  3.04       0.336     0.159
#> 10 group4 vs group5  0.277 0.00573  2.90       0.341     0.164
#> 
#> -----------------------------------------------------
#> Use plot() to visualize

# Pooled histogram of |theta_i - theta_j| across all pairs
plot(pw)

The prob_gt_0.5 and prob_gt_1 columns in the summary show the prior probability of observing pairwise differences exceeding those thresholds — useful for assessing whether your prior is consistent with your clinical expectations about subgroup heterogeneity.

Step 3: Fit the Hierarchical Model

fit <- shrink(
  mixture             = mix,
  hierarchical_priors = hierarchical_priors,
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 456,
  refresh = 0
)
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 8e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 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.093 seconds (Warm-up)
#> Chain 1:                0.087 seconds (Sampling)
#> Chain 1:                0.18 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 9e-06 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 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.097 seconds (Warm-up)
#> Chain 2:                0.059 seconds (Sampling)
#> Chain 2:                0.156 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.096 seconds (Warm-up)
#> Chain 3:                0.072 seconds (Sampling)
#> Chain 3:                0.168 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 7e-06 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.07 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.104 seconds (Warm-up)
#> Chain 4:                0.078 seconds (Sampling)
#> Chain 4:                0.182 seconds (Total)
#> Chain 4:

print(fit)
#> # A tibble: 3 × 7
#>   variable      mean     sd      q2.5     q50 q97.5  rhat
#>   <chr>        <dbl>  <dbl>     <dbl>   <dbl> <dbl> <dbl>
#> 1 mu          0.583  0.0890 0.402     0.586   0.763  1.00
#> 2 tau         0.109  0.0973 0.00380   0.0827  0.368  1.00
#> 3 tau_squared 0.0214 0.0408 0.0000144 0.00684 0.136  1.00

Step 4: Examine Results

mu_tau <- extract_mu_tau(fit)

cat("Overall treatment effect (mu):\n")
#> Overall treatment effect (mu):
cat("  Mean:", round(mean(mu_tau$mu), 3), "\n")
#>   Mean: 0.583
cat("  95% CI:", round(quantile(mu_tau$mu, c(0.025, 0.975)), 3), "\n\n")
#>   95% CI: 0.402 0.763

cat("Between-region heterogeneity (tau):\n")
#> Between-region heterogeneity (tau):
cat("  Mean:", round(mean(mu_tau$tau), 3), "\n")
#>   Mean: 0.109
cat("  95% CI:", round(quantile(mu_tau$tau, c(0.025, 0.975)), 3), "\n")
#>   95% CI: 0.004 0.368

theta_summary <- summarize_theta(fit)
print(theta_summary)
#> # A tibble: 5 × 9
#>   group    mean     sd  q2.5   q50 q97.5  rhat ess_bulk ess_tail
#>   <chr>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 North   0.590 0.0980 0.400 0.590 0.790  1.00    4037.    2233.
#> 2 South   0.587 0.101  0.384 0.587 0.792  1.00    3329.    1901.
#> 3 East    0.596 0.0879 0.427 0.594 0.778  1.00    4630.    3305.
#> 4 West    0.621 0.102  0.434 0.612 0.851  1.00    3943.    3258.
#> 5 Central 0.523 0.121  0.244 0.538 0.721  1.00    3389.    2517.

Step 5: Visualize Shrinkage

plot(fit)

Key observations:

Central (most uncertain) shrinks most toward the mean
East (most precise) shrinks least
This is adaptive shrinkage: uncertain estimates borrow more

plot(fit, type = "diagnostics")

Alternative Input: Using Summary Statistics Only

If you only have published means and standard errors (no full posteriors), you can use the MLE approach:

mle_estimates <- sapply(samples, mean)
mle_variances <- sapply(samples, var)

fit_mle <- shrink(
  mle                 = mle_estimates,
  var_matrix          = mle_variances,
  hierarchical_priors = hierarchical_priors,
  chains  = 4,
  iter    = 2000,
  warmup  = 1000,
  seed    = 456,
  refresh = 0
)
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 8e-06 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.08 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.087 seconds (Warm-up)
#> Chain 1:                0.076 seconds (Sampling)
#> Chain 1:                0.163 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 8e-06 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.08 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.098 seconds (Warm-up)
#> Chain 2:                0.096 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 7e-06 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 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.084 seconds (Warm-up)
#> Chain 3:                0.082 seconds (Sampling)
#> Chain 3:                0.166 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'stage2_shrinkage' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 1.1e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.11 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.059 seconds (Sampling)
#> Chain 4:                0.16 seconds (Total)
#> Chain 4:

mu_tau_mle <- extract_mu_tau(fit_mle)

cat("Mixture approach:\n")
#> Mixture approach:
cat("  mu =",  round(mean(mu_tau$mu), 3), "\n")
#>   mu = 0.583
cat("  tau =", round(mean(mu_tau$tau), 3), "\n\n")
#>   tau = 0.109

cat("MLE approach:\n")
#> MLE approach:
cat("  mu =",  round(mean(mu_tau_mle$mu), 3), "\n")
#>   mu = 0.579
cat("  tau =", round(mean(mu_tau_mle$tau), 3), "\n")
#>   tau = 0.112

Method	Use When	Pros	Cons
Mixture	You have full posteriors	Captures non-normality; More accurate	Requires posterior samples
MLE	Only means + SEs available	Simpler; Works with published data	Assumes normality

Complete Stan → shrinkr Workflow Summary

# 1. Write Stan model with FLAT PRIORS on parameters of interest
stan_code <- "
data {
  int<lower=0> N;
  int<lower=1> G;
  vector[N] y;
  vector[N] treatment;
  array[N] int<lower=1,upper=G> group;
}
parameters {
  vector[G] theta;  // Group-specific effects - NO PRIOR SPECIFIED
  real<lower=0> sigma;
}
model {
  sigma ~ normal(0, 2);
  for (n in 1:N) {
    y[n] ~ normal(treatment[n] * theta[group[n]], sigma);
  }
}
"

# 2. Fit model once to get all group effects, extract posteriors
group_indices <- as.integer(factor(data$group, levels = groups))
fit_stan <- stan(
  model_code = stan_code,
  data = list(N = nrow(data), G = length(groups),
              y = data$y, treatment = data$treatment,
              group = group_indices),
  chains = 4, iter = 2000, warmup = 1000, refresh = 0
)

theta_draws <- extract(fit_stan)$theta
samples <- lapply(seq_along(groups), function(i) matrix(theta_draws[, i], ncol = 1))
names(samples) <- groups

# 3. Fit mixture approximation and check quality
mix <- fit_mixture(samples, K_max = 3)
plot(mix, draws = samples)

# 4. Specify and check hierarchical priors
priors <- list(
  mu  = dist_normal(0, 5),
  tau = dist_truncated(dist_student_t(3, 0, 1), lower = 0)
)
plot(sample_prior_predictive(priors, n_groups = length(groups)))

# 5. Fit, extract, and visualize
fit <- shrink(mixture = mix, hierarchical_priors = priors)
plot(fit)
summarize_theta(fit)
extract_mu_tau(fit)

Key Concepts Checklist

Stage 1 uses flat priors - This is critical! Don’t specify a prior on the parameter you’ll shrink
Mixture approximation - Allows shrinkr to work with any posterior shape; check quality with plot()
Hierarchical priors - Applied only in Stage 2; use sample_prior_predictive() to understand implications
Adaptive shrinkage - Uncertain estimates borrow more; precise estimates borrow less
Diagnostics matter - Always check mixture fit, prior-data agreement, and MCMC convergence

Common Pitfalls to Avoid

Using informative priors in Stage 1 - This creates “double priors” and breaks the equivalence
Skipping mixture quality checks - Poor approximation → biased shrinkage, run: plot(mix, draws = samples)
Ignoring prior implications - Your prior might be too strong or weak, run: sample_prior_predictive() and plot()
Not checking convergence - MCMC diagnostics apply to Stage 2 too, check: fit$diagnostics, Rhat, ESS

Next Steps

Advanced integration: See vignette("tidy_bayesian_workflow") for working with tidyverse tools
Real applications: See vignette("brms_integration") for survival analysis and brms integration
Sensitivity analysis: Learn to explore different prior specifications efficiently (this is done in the brms vignette)
Federated learning: See vignette("federated_learning") for distributed data analysis

Session Info

sessionInfo()
#> R version 4.4.2 (2024-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 10 x64 (build 19045)
#> 
#> Matrix products: default
#> 
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/Chicago
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] patchwork_1.3.2      posterior_1.7.0      survival_3.7-0      
#>  [4] lubridate_1.9.5      forcats_1.0.1        stringr_1.6.0       
#>  [7] dplyr_1.2.1          purrr_1.2.2          readr_2.2.0         
#> [10] tidyr_1.3.2          tibble_3.3.1         ggplot2_4.0.3       
#> [13] tidyverse_2.0.0      distributional_0.7.1 tidybayes_3.0.7     
#> [16] brms_2.23.0          Rcpp_1.1.1           shrinkr_0.4.5       
#> 
#> loaded via a namespace (and not attached):
#>  [1] tidyselect_1.2.1      svUnit_1.0.8          farver_2.1.2         
#>  [4] loo_2.9.0             S7_0.2.2              fastmap_1.2.0        
#>  [7] tensorA_0.36.2.1      digest_0.6.39         estimability_1.5.1   
#> [10] timechange_0.4.0      lifecycle_1.0.5       StanHeaders_2.32.10  
#> [13] magrittr_2.0.5        compiler_4.4.2        rlang_1.2.0          
#> [16] sass_0.4.10           tools_4.4.2           utf8_1.2.6           
#> [19] yaml_2.3.12           knitr_1.51            labeling_0.4.3       
#> [22] bridgesampling_1.2-1  pkgbuild_1.4.8        mclust_6.1.2         
#> [25] curl_7.1.0            RColorBrewer_1.1-3    abind_1.4-8          
#> [28] withr_3.0.2           grid_4.4.2            stats4_4.4.2         
#> [31] xtable_1.8-8          inline_0.3.21         emmeans_2.0.3        
#> [34] scales_1.4.0          cli_3.6.6             mvtnorm_1.4-1        
#> [37] rmarkdown_2.31        generics_0.1.4        otel_0.2.0           
#> [40] RcppParallel_5.1.11-2 rstudioapi_0.19.0     tzdb_0.5.0           
#> [43] cachem_1.1.0          rstan_2.32.7          splines_4.4.2        
#> [46] bayesplot_1.15.0      parallel_4.4.2        matrixStats_1.5.0    
#> [49] vctrs_0.7.3           V8_8.2.0              Matrix_1.7-1         
#> [52] jsonlite_2.0.0        hms_1.1.4             arrayhelpers_1.1-0   
#> [55] ggdist_3.3.3          jquerylib_0.1.4       glue_1.8.1           
#> [58] codetools_0.2-20      stringi_1.8.7         gtable_0.3.6         
#> [61] QuickJSR_1.9.2        pillar_1.11.1         htmltools_0.5.9      
#> [64] Brobdingnag_1.2-9     R6_2.6.1              evaluate_1.0.5       
#> [67] lattice_0.22-6        backports_1.5.1       bslib_0.11.0         
#> [70] rstantools_2.6.0      coda_0.19-4.1         gridExtra_2.3        
#> [73] nlme_3.1-166          checkmate_2.3.4       xfun_0.57            
#> [76] pkgconfig_2.0.3