# 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)
)
)
