Central Coordinator: Stage 2 Shrinkage
The coordinator, for example a coordinating center or trusted third party, now applies hierarchical shrinkage.
Path A: Using Full Posteriors
# Fit mixture approximation
mix <- fit_mixture(
samples = site_posteriors,
K_max = 2, # Age effects should be fairly normal
verbose = TRUE
)
# Check quality
plot(mix, draws = site_posteriors, type = "density")
Specify Network-Level Priors
# Based on clinical knowledge:
# - Age effect should be positive but moderate
# - Some heterogeneity expected across hospital types
hierarchical_priors <- list(
mu = dist_normal(0.05, 0.025), # Centered on 5% increase per year
tau = dist_truncated(dist_student_t(3, 0, 0.01), lower = 0) # Modest heterogeneity
)
# Visualize prior implications
prior_pred <- sample_prior_predictive(
hierarchical_priors = hierarchical_priors,
n_groups = 6,
n_draws = 1000
)
plot(prior_pred, type = "both")
Fit Hierarchical Model
fit_full_post <- shrink(
mixture = mix,
hierarchical_priors = hierarchical_priors,
chains = 4,
iter = 2000,
warmup = 1000,
seed = 123
)print(fit_full_post)
#> # A tibble: 3 × 7
#> variable mean sd q2.5 q50 q97.5 rhat
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 mu 0.0497 0.00750 0.0349 0.0497 0.0643 1.00
#> 2 tau 0.00624 0.00536 0.000255 0.00498 0.0198 1.00
#> 3 tau_squared 0.0000677 0.000132 0.0000000650 0.0000248 0.000390 1.00
# Network-level estimates
mu_tau_full <- extract_mu_tau(fit_full_post)
cat("\nNetwork-level age effect (mu):\n")
#>
#> Network-level age effect (mu):
cat(" Posterior mean:", round(mean(mu_tau_full$mu), 4), "\n")
#> Posterior mean: 0.0497
cat(" 95% CI: [", round(quantile(mu_tau_full$mu, 0.025), 4), ",",
round(quantile(mu_tau_full$mu, 0.975), 4), "]\n")
#> 95% CI: [ 0.0349 , 0.0643 ]
cat("\nBetween-site heterogeneity (tau):\n")
#>
#> Between-site heterogeneity (tau):
cat(" Posterior mean:", round(mean(mu_tau_full$tau), 4), "\n")
#> Posterior mean: 0.0062
cat(" 95% CI: [", round(quantile(mu_tau_full$tau, 0.025), 4), ",",
round(quantile(mu_tau_full$tau, 0.975), 4), "]\n")
#> 95% CI: [ 3e-04 , 0.0198 ]Path B: Using Only Summary Statistics
Before using summary statistics, we must verify posteriors are approximately normal.
Step 1: Check Normality Assumption
# Visual checks for approximate normality
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(3, 2))
for(i in 1:6) {
site_name <- names(site_posteriors)[i]
samples_i <- as.vector(site_posteriors[[i]])
# QQ plot against normal
qqnorm(samples_i, main = paste("QQ Plot:", site_name))
qqline(samples_i, col = "red", lwd = 2)
}
# Quantitative checks
normality_checks <- data.frame(
site = names(site_posteriors),
skewness = sapply(site_posteriors, function(x) {
m3 <- mean((x - mean(x))^3)
s3 <- sd(x)^3
m3 / s3
}),
kurtosis = sapply(site_posteriors, function(x) {
m4 <- mean((x - mean(x))^4)
s4 <- sd(x)^4
m4 / s4 - 3 # Excess kurtosis
})
)
print(normality_checks)
#> site skewness kurtosis
#> Metro General Metro General 0.006067498 0.009199907
#> County Regional County Regional -0.050758390 -0.116123944
#> University Medical University Medical 0.035950881 -0.005600149
#> Community Hospital Community Hospital 0.008651589 0.118037106
#> Veterans Affairs Veterans Affairs -0.041185160 -0.015861386
#> Children's Specialty Children's Specialty 0.010783310 -0.090076276
cat("\nNormality assessment:\n")
#>
#> Normality assessment:
cat(" Skewness close to 0? (|skew| < 0.5 is good)\n")
#> Skewness close to 0? (|skew| < 0.5 is good)
cat(" Kurtosis close to 0? (|kurt| < 1.0 is good)\n")
#> Kurtosis close to 0? (|kurt| < 1.0 is good)
cat(" All sites pass:",
all(abs(normality_checks$skewness) < 0.5 & abs(normality_checks$kurtosis) < 1.0),
"\n")
#> All sites pass: TRUEDecision rule:
- If posteriors look approximately normal, Path B is likely valid
- If posteriors are skewed/heavy-tailed, use Path A
- If unsure, use Path A
When CLT Fails: A Counter-Example
To illustrate why normality matters, consider a scenario where we are estimating a variance parameter:
# Simulated: posterior for a variance parameter (boundary at 0)
set.seed(999)
variance_posterior <- matrix(rchisq(4000, df = 5) / 5, ncol = 1)
# Compute summary statistics
var_mean <- mean(variance_posterior)
var_se <- sd(variance_posterior)
# Check normality
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1, 2))
hist(variance_posterior, breaks = 30, main = "Variance Posterior",
xlab = "sigma^2", col = "lightblue")
qqnorm(variance_posterior, main = "QQ Plot")
qqline(variance_posterior, col = "red", lwd = 2)
par(oldpar)
# Skewness
skew <- mean((variance_posterior - var_mean)^3) / var_se^3
cat("Skewness:", round(skew, 2), "(should be about 0 for normal)\n")
#> Skewness: 1.24 (should be about 0 for normal)
cat("This posterior is right-skewed - CLT approximation would be poor!\n")
#> This posterior is right-skewed - CLT approximation would be poor!For this scenario:
- Path B summaries would give biased results
- Path A mixture handles the skewness correctly
Step 2: Fit Using CLT Approximation
# Extract means and variances
mle_estimates <- site_summaries$beta_age_mean
names(mle_estimates) <- site_summaries$site
mle_variances <- site_summaries$beta_age_se^2
names(mle_variances) <- site_summaries$site
# Fit using MLE path (CLT approximation)
fit_summaries <- shrink(
mle = mle_estimates,
var_matrix = mle_variances,
hierarchical_priors = hierarchical_priors,
chains = 4,
iter = 2000,
warmup = 1000,
seed = 123,
verbose = FALSE,
refresh = 0
)
print(fit_summaries)
#> # A tibble: 3 × 7
#> variable mean sd q2.5 q50 q97.5 rhat
#> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 mu 0.0496 0.00732 0.0349 0.0497 0.0637 1.00
#> 2 tau 0.00627 0.00545 0.000209 0.00496 0.0201 1.00
#> 3 tau_squared 0.0000690 0.000136 0.0000000438 0.0000246 0.000402 1.00Compare Paths
mu_tau_summaries <- extract_mu_tau(fit_summaries)
comparison <- data.frame(
parameter = c("mu", "tau"),
full_posteriors = c(
mean(mu_tau_full$mu),
mean(mu_tau_full$tau)
),
summaries_only = c(
mean(mu_tau_summaries$mu),
mean(mu_tau_summaries$tau)
)
) %>%
mutate(difference = abs(full_posteriors - summaries_only))
print(comparison)
#> parameter full_posteriors summaries_only difference
#> 1 mu 0.049654685 0.049563051 9.163467e-05
#> 2 tau 0.006243329 0.006270951 2.762236e-05
cat("\nMaximum difference:", round(max(comparison$difference), 5), "\n")
#>
#> Maximum difference: 9e-05Conclusion: Both paths give nearly identical results because posteriors are approximately normal in this case. This will not always be true.
When each path is appropriate:
| Situation | Recommended Path | Reason |
|---|---|---|
| Posteriors are normal (verified) | Path B acceptable | CLT holds; minimal sharing |
| Posteriors are skewed/multimodal | Path A required | CLT fails; mixture needed |
| Small sample sizes per site | Path A safer | CLT may not hold yet |
| Boundary constraints | Path A required | CLT assumes interior parameters |
| Unknown posterior shape | Path A safer | Conservative choice |
| Maximum privacy needed and normal posteriors | Path B acceptable | But verify normality |
