Package {shrinkr}

Type:	Package
Title:	Modular Bayesian Hierarchical Shrinkage Models
Version:	0.4.5
Date:	2026-06-29
Description:	Implements a two-stage Bayesian hierarchical modeling framework for applying shrinkage to subgroup-specific effects. The package separates model fitting (Stage 1) from hierarchical shrinkage (Stage 2), enabling modular sensitivity analyses without refitting expensive Markov chain Monte Carlo (MCMC) chains. Supports flexible prior specifications through the 'distributional' package, mixture approximations via 'mclust', and efficient 'Stan'-based inference.
License:	GPL (≥ 3)
Encoding:	UTF-8
Depends:	R (≥ 4.1)
Imports:	rstan (≥ 2.32.7), mclust, posterior, distributional, cli, tibble, dplyr, purrr, stats, utils, methods, Rcpp (≥ 0.12.0), RcppParallel (≥ 5.0.1), rstantools (≥ 2.5.0)
Suggests:	testthat (≥ 3.0.0), ggplot2, knitr, rmarkdown, patchwork, tidybayes, bayesplot, tidyverse, ggdist, ggridges, brms, tidyr, survival, MASS, Matrix, beastt
LinkingTo:	BH (≥ 1.66.0), Rcpp (≥ 0.12.0), RcppEigen (≥ 0.3.3.3.0), RcppParallel (≥ 5.0.1), rstan (≥ 2.32.7), StanHeaders (≥ 2.32.10)
SystemRequirements:	GNU make, C++17 compiler
Biarch:	true
URL:	https://gsk-biostatistics.github.io/shrinkr/
BugReports:	https://github.com/GSK-Biostatistics/shrinkr/issues
VignetteBuilder:	knitr
NeedsCompilation:	yes
Config/roxygen2/version:	8.0.0
Packaged:	2026-06-29 20:46:16 UTC; jmaro
Author:	Jacob M. Maronge [aut, cre], GlaxoSmithKline Research & Development Limited [cph, fnd], Trustees of Columbia University [cph] (R/stanmodels.R, configure, configure.win)
Maintainer:	Jacob M. Maronge <jacob.m.maronge@gsk.com>
Repository:	CRAN
Date/Publication:	2026-07-06 13:40:02 UTC

shrinkr: Modular Bayesian Hierarchical Shrinkage Models

Description

The shrinkr package provides a flexible framework for two-stage Bayesian hierarchical modeling. It enables post-hoc shrinkage of subgroup-specific posterior estimates from any Bayesian model, with support for diverse prior specifications and diagnostic tools.

Key Features

Two-Stage Workflow:

• Stage 1: Fit any Bayesian model without shrinkage

• Stage 2: Apply hierarchical shrinkage with flexible priors

Flexible Priors:

• Standard families (Normal, Student-t, Cauchy, Lognormal)

• Heavy-tailed (Inverse-Gamma, Half-Cauchy, Half-t)

• Bounded (Uniform)

• Mixture priors (spike-and-slab)

• Truncated distributions

Input Methods:

• Full posterior samples (via mixture approximation)

• Point estimates + variance/covariance

Main Functions

Core Workflow:

• fit_mixture(): Approximate Stage 1 posteriors with Gaussian mixture

• shrink(): Main user interface for hierarchical shrinkage

Prior Specification:

• prior_spike_slab(): Create spike-and-slab mixture prior

• prior_mixture(): Create custom mixture prior

• sample_prior_predictive(): Generate prior predictive samples for checking

Extraction & Visualization:

• extract_mu_tau(): Extract hyperparameter draws

• extract_theta(): Extract group-level draws

• summarise_mu_tau(): Summarize hyperparameters

• summarise_theta(): Summarize group-level estimates

• theta_contrasts(): Compute linear combinations of theta

• plot(): Visualize shrinkage effect and mixture approximation quality

Getting Started

See vignette("getting_started", package = "shrinkr") for a basic workflow, or vignette("brms_integration", package = "shrinkr") for a survival analysis example.

Use Cases

• Meta-analysis: Shrink study-specific effects

• Clinical trials: Borrow information across subgroups or historical controls

• Genomics: Regularize gene-specific effects

• Simulation studies: Compare shrinkage methods systematically

Package Options

• shrinkr.refresh: Controls Stan sampling progress output (default: 100)

Author(s)

Maintainer: Jacob M. Maronge jacob.m.maronge@gsk.com (ORCID)

Authors:

• Jacob M. Maronge jacob.m.maronge@gsk.com (ORCID)

Other contributors:

• GlaxoSmithKline Research & Development Limited [copyright holder, funder]

• Trustees of Columbia University (R/stanmodels.R, configure, configure.win) [copyright holder]

References

Maronge, J. M. (2026). shrinkr: Modular Bayesian Hierarchical Shrinkage Models. R package version 0.4.3.

See Also

• Stan: https://mc-stan.org/

• distributional package: https://pkg.mitchelloharawild.com/distributional/

Examples

## Not run: 
# This example fits a Stan model, so it is not run during package checks.
library(shrinkr)
priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_student_t(3, 0, 1), lower = 0)
)
fit <- shrink(
  mle = c(0.0, 0.5, 1.0),
  var_matrix = c(0.25, 0.25, 0.25),
  hierarchical_priors = priors,
  iter = 1000, chains = 2, seed = 1
)
summary(fit)

## End(Not run)

Internal: convert covariance to Cholesky factor for Stan

Description

Uses a jitter-and-retry strategy when the input matrix is not quite positive-definite. Falls back to Matrix::nearPD() if available, and errors rather than passing a non-lower-triangular factor to Stan.

Usage

.as_chol_factor(Sigma)

Convert distributional priors to Stan format

Description

Handles Normal, Student-t, Cauchy, Gamma, Exponential, Lognormal, Inverse-Gamma, Uniform, truncated distributions, and mixture priors (including spike-and-slab) for both mu and tau.

Usage

.coerce_priors_to_stan(prior_mu, prior_tau)

Prepare Stan data from mixture

Description

Prepare Stan data from mixture

Usage

.prep_bhm_data_from_mixture(mixture, pri, centered)

Prepare Stan data from mle and cov

Description

Prepare Stan data from mle and cov

Usage

.prep_bhm_data_from_mle(mle, var_matrix, pri, centered)

Internal: stop if fewer than 2 groups

Description

Internal: stop if fewer than 2 groups

Usage

.stop_if_lt_two_groups(groups)

Convert shrinkr_fit to data.frame

Description

Extracts posterior draws as a regular data frame. This is a convenience wrapper around as_draws_df() that returns a plain data.frame.

Usage

## S3 method for class 'shrinkr_fit'
as.data.frame(
  x,
  row.names = NULL,
  optional = FALSE,
  variables = NULL,
  include_internals = FALSE,
  ...
)

Arguments

Value

A data.frame with columns for chain, iteration, draw, and requested parameters.

See Also

as_draws_df.shrinkr_fit() for posterior package format, extract_mu_tau() for hyperparameters only

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
draws_df <- as.data.frame(fit)
head(draws_df)
mu_tau_df <- as.data.frame(fit, variables = c("mu", "tau"))

Convert mixture fit to data frame

Description

Converts a fitted mixture model to a tidy long-format data frame. This is essentially an accessor for the components element.

Usage

## S3 method for class 'shrinkr_mixture'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A data frame (or tibble if available) containing the component specifications with columns:

See Also

fit_mixture for fitting mixture models

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
df <- as.data.frame(mix)
head(df)

Convert prior predictive samples to data frame

Description

Converts prior predictive samples to a tidy long-format data frame suitable for analysis and visualization with tidyverse tools.

Usage

## S3 method for class 'shrinkr_prior_pred'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A data frame (or tibble if tibble package is available) with columns:

See Also

sample_prior_predictive for generating prior predictive samples

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
df <- as.data.frame(prior_pred)
head(df)

Convert shrinkr_fit to draws_df

Description

Extracts posterior draws in tidy format using the posterior package. By default returns user-facing parameters (mu, tau, theta, etc.) and excludes internal parameterization details. Set include_internals = TRUE to access all parameters including theta_c and z.

Usage

## S3 method for class 'shrinkr_fit'
as_draws_df(x, variables = NULL, include_internals = FALSE, ...)

Arguments

Value

A posterior::draws_df with columns for chain, iteration, draw, and requested parameters.

See Also

shrink() for fitting models, extract_mu_tau() for hyperparameters only

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
all_draws <- posterior::as_draws_df(fit)
posterior::variables(all_draws)
theta_draws <- posterior::as_draws_df(fit, variables = c("theta[1]", "theta[2]"))

Extract mu and tau parameters

Description

Extracts posterior draws for the hyperparameters mu (global mean) and tau (heterogeneity standard deviation) from a fitted shrinkage model.

Usage

extract_mu_tau(x, ...)

Arguments

Value

A posterior::draws_df with columns:

.chain

Chain index

.iteration

Iteration within chain

.draw

Overall draw index

mu

Global mean parameter

tau

Heterogeneity parameter

tau_squared

Variance (tau^2)

See Also

shrink() for fitting models, summarise_mu_tau() for summary statistics, as_draws_df.shrinkr_fit() for all parameters

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
mu_tau <- extract_mu_tau(fit)
summarise_mu_tau(fit)
posterior::summarise_draws(mu_tau)

Extract theta (group-level effect) parameters

Description

Extracts posterior draws for the group-level effects (theta parameters) from a fitted shrinkage model. This is the hierarchically shrunk version of the subgroup effects.

Usage

extract_theta(x, ...)

Arguments

Value

A posterior::draws_df with columns:

.chain

Chain index

.iteration

Iteration within chain

.draw

Overall draw index

theta[1], theta[2], ...

Group-level effects

See Also

shrink() for fitting models, extract_mu_tau() for hyperparameters, summarise_theta() for summary statistics, theta_contrasts() for pairwise comparisons

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
theta_draws <- extract_theta(fit)
posterior::summarise_draws(theta_draws)
summarise_theta(fit)

Fit Gaussian mixture models to posterior samples

Description

Fits a multivariate Gaussian mixture model (GMM) jointly across all supplied variables using mclust. The function is intended for approximating posterior draws from Bayesian models and produces a tidy component table suitable for visualization and shrinkage modeling.

Usage

fit_mixture(samples, K_max = 5L, verbose = FALSE, model_names = NULL, ...)

Arguments

Details

• Requires at least two variables; shrinkage across a single variable is not meaningful.

• Rows with any missing values are removed before fitting.

• Component weights are normalized to sum to one.

• The sd column reports marginal standard deviations (square roots of diagonal entries) from each component covariance matrix.

Value

A list with class "shrinkr_mixture" containing:

• components — data frame with columns: group, component, variable, weight, mean, and sd (marginal standard deviations).

• K — number of mixture components selected.

• vars — character vector of variable names.

• weights — vector of component weights (mixing proportions).

• covs — list of component covariance matrices (p × p each).

• model_name — selected mclust covariance structure (e.g., "VVV").

• bic — Bayesian Information Criterion for the fitted model.

• n_samples — number of samples used in fitting.

• n_vars — number of variables.

• mclust_fit — the complete mclust model object (for advanced use).

• diagnostics — list with sample size details, removed rows, and quality warnings.

Covariance structures in mclust

mclust parameterizes component covariances via eigen-decomposition and offers a set of model families controlling volume (V), shape (S), and orientation (O). Common codes include (non-exhaustive):

Spherical

"EII": equal volume, spherical
"VII": variable volume, spherical

Diagonal

"EEI": equal volume & shape (axis-aligned)
"VEI": variable volume, equal shape
"EVI": equal volume, variable shape
"VVI": variable volume & shape

Ellipsoidal (full covariance)

"EEE": equal volume, shape, orientation
"EEV": equal volume & shape, variable orientation
"VEV": variable volume, equal shape, variable orientation
"VVV": variable volume, shape, and orientation (most flexible)

If model_names = NULL (default), mclust::Mclust() selects among the models appropriate for the data dimension via BIC. In practice, this lets the data decide between parsimonious structures (e.g., "EII", "VVI") and fully flexible ones (e.g., "VVV"). You can restrict or expand the search space by supplying model_names.

See Also

plot.shrinkr_mixture for visualizing marginal fits, as.data.frame.shrinkr_mixture for extracting component data

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1),
  group3 = matrix(rnorm(100, 1.0, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
summary(mix)

Internal null coalescing operator

Description

Returns y if x is NULL, otherwise returns x. Used internally to simplify default argument handling.

Usage

x %||% y

Plot shrinkage fit

Description

Visualizes the hierarchical shrinkage model fit. Creates either:

• "shrinkage" — Shows pre/post-shrunk estimates with arrows

• "diagnostics" — Multi-panel view with hyperparameters and shrinkage factor

The shrinkage plot displays:

• Pre-shrunk estimates (hollow circles) — from stage 1 mixture or MLEs

• Post-shrunk estimates (filled circles) — posterior means of theta

• Global mean (dashed line) — posterior mean of mu

• Credible intervals (optional) — for shrunken estimates

• Arrows (optional) — showing direction and magnitude of shrinkage

Usage

## S3 method for class 'shrinkr_fit'
plot(
  x,
  type = c("shrinkage", "diagnostics"),
  group_names = NULL,
  show_arrows = FALSE,
  show_intervals = TRUE,
  interval_prob = 0.95,
  point_size = 3,
  arrow_alpha = 0.6,
  dodge_width = 0.3,
  title = NULL,
  subtitle = NULL,
  ...
)

Arguments

Value

A ggplot2 object (for type = "shrinkage"), or a patchwork object/list (for type = "diagnostics").

See Also

shrink() for fitting models, extract_mu_tau() for hyperparameter draws

Examples

# Plotting requires a fitted shrinkr_fit object from shrink().
# The full example is not run because it fits a Stan model.
## Not run: 
fit <- shrink(mixture = mix, hierarchical_priors = priors)
plot(fit)
plot(fit, show_arrows = TRUE, interval_prob = 0.95)

## End(Not run)

Plot fitted marginal densities or QQ plots for mixture models

Description

Overlays fitted marginal mixture densities from a shrinkr_mixture (returned by fit_mixture()) on top of the observed samples for selected variables, OR creates QQ plots comparing empirical vs fitted quantiles. The function uses the same coercion logic as fit_mixture() (via an internal .coerce_draws_df() helper), ensuring that variable names line up even when users pass a list of matrices.

Important: For multivariate joint fits, this produces marginal overlays (one panel per variable when faceting). Each marginal density is computed by summing the weighted component densities for that variable.

Usage

## S3 method for class 'shrinkr_mixture'
plot(
  x,
  draws = NULL,
  variables = NULL,
  type = c("density", "qq"),
  overlay = c("hist", "kde", "both", "none"),
  bins = 50,
  kde_bw = NULL,
  show_components = TRUE,
  facet = TRUE,
  n_points = 501,
  verbose = FALSE,
  ...
)

Arguments

Details

Density plots

The total marginal density for each variable j is computed as

f_j(x)=\sum_{k=1}^{K} w_k \,\phi\!\left(x \mid \mu_{jk},\, \sigma_{jk}\right),

using per-component marginal SDs (sd) already stored in x$components.

The plotting range per variable is taken from the sample range if available (with 5% padding), otherwise from ⁠mean +/- 4*sd⁠ across that variable's components–avoiding non-finite seq() errors when samples are absent.

QQ plots

When type = "qq", the function creates quantile-quantile plots by:

1. Computing empirical quantiles from the observed data

2. Computing theoretical quantiles from the fitted mixture CDF via numerical inversion

3. Plotting empirical vs theoretical quantiles with a 45-degree reference line

Points falling on the reference line indicate good agreement between the fitted mixture and the data. Systematic deviations suggest model misfit.

Value

A ggplot2 object.

See Also

fit_mixture() for fitting mixture models

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
plot(mix, draws = samples, type = "density", variables = c("group1", "group2"))

Plot prior predictive pairwise differences

Description

Creates a density plot of |\theta_i - \theta_j| from prior predictive samples. Useful for calibrating hierarchical priors. Styling matches plot.shrinkr_prior_pred() for visual consistency.

Usage

## S3 method for class 'shrinkr_prior_contrasts'
plot(x, by_pair = FALSE, ...)

Arguments

Value

A ggplot2 object.

See Also

prior_pairwise_differences

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
pw <- prior_pairwise_differences(prior_pred)
plot(pw)

Plot prior predictive samples

Description

Visualizes the prior predictive distribution for hyperparameters (mu and tau) and subgroup effects (theta). Requires the ggplot2 package.

Usage

## S3 method for class 'shrinkr_prior_pred'
plot(x, type = c("both", "hyperparameters", "theta"), ...)

Arguments

Value

A ggplot2 object, or a list of two ggplot2 objects if type = "both" and patchwork package is not available. If patchwork is available and type = "both", returns a combined plot.

See Also

sample_prior_predictive for generating samples, as.data.frame.shrinkr_prior_pred for extracting data

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_student_t(3, 0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
plot(prior_pred, type = "hyperparameters")

Print method for shrinkr_fit

Description

Displays a compact summary of the fitted model including dimensions, hyperparameter estimates, and diagnostics.

Usage

## S3 method for class 'shrinkr_fit'
print(x, digits = 3, ...)

Arguments

Value

Invisibly returns the input object x.

See Also

shrink() for fitting models, summarise_theta() for detailed theta estimates

Print method for mixture fits

Description

Displays summary information about a fitted mixture model in a readable format.

Usage

## S3 method for class 'shrinkr_mixture'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

See Also

fit_mixture for fitting mixture models, summary.shrinkr_mixture for detailed summaries

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
print(mix)

Print method for prior pairwise contrasts

Description

Print method for prior pairwise contrasts

Usage

## S3 method for class 'shrinkr_prior_contrasts'
print(x, digits = 3, ...)

Arguments

Value

Invisibly returns x.

Print method for prior predictive samples

Description

Displays summary information about prior predictive samples in a readable format.

Usage

## S3 method for class 'shrinkr_prior_pred'
print(x, ...)

Arguments

Value

Invisibly returns the input object x.

See Also

sample_prior_predictive for generating samples, summary.shrinkr_prior_pred for detailed summaries

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
print(prior_pred)

Print summary of mixture fit

Description

Displays formatted summary statistics for a fitted mixture model.

Usage

## S3 method for class 'summary.shrinkr_mixture'
print(x, digits = 3, ...)

Arguments

Value

Invisibly returns the input object x.

See Also

summary.shrinkr_mixture for creating summaries

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
summ <- summary(mix)
print(summ)

Print summary of prior predictive samples

Description

Displays formatted summary statistics for prior predictive samples.

Usage

## S3 method for class 'summary.shrinkr_prior_pred'
print(x, digits = 3, ...)

Arguments

Value

Invisibly returns the input object x.

See Also

summary.shrinkr_prior_pred for creating summaries

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
print(summary(prior_pred))

Create a mixture prior

Description

Creates a mixture distribution using distributional::dist_mixture(). All standard distributional operations (sampling, density, quantiles, formatting) work automatically.

Usage

prior_mixture(..., weights = NULL)

Arguments

Value

A distributional mixture distribution object

Examples

mix <- prior_mixture(
  distributional::dist_normal(0, 0.1),
  distributional::dist_normal(0, 1),
  weights = c(0.7, 0.3)
)

Compute prior predictive pairwise differences |theta_i - theta_j|

Description

Computes the prior-implied distribution of absolute pairwise differences between subgroup effects. This is useful for calibrating priors: if your prior implies that subgroup differences of 5 units are common but clinical relevance starts at 0.5, your prior may be too diffuse.

This implements the recommendation from the SDSIH Vignettes Library: inspect the prior distribution of |\theta_i - \theta_j| when choosing priors for Bayesian hierarchical models.

Usage

prior_pairwise_differences(prior_pred)

Arguments

Value

A list with class "shrinkr_prior_contrasts" containing:

See Also

sample_prior_predictive for generating prior predictive samples, plot.shrinkr_prior_contrasts for visualizing the result

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
pw <- prior_pairwise_differences(prior_pred)
print(pw)
pw$overall_summary

Spike-and-slab prior for testing homogeneity

Description

Creates a mixture of two Normal distributions. Since tau is a scale parameter, this must be wrapped in distributional::dist_truncated() with lower = 0 before passing to shrink():

Usage

prior_spike_slab(
  spike_location = 0,
  spike_scale = 0.01,
  slab_scale = 1,
  spike_prob = 0.5
)

Arguments

Details

tau_prior <- dist_truncated(prior_spike_slab(), lower = 0)

Value

A spike-and-slab mixture distribution

Examples

tau_prior <- distributional::dist_truncated(
  prior_spike_slab(spike_prob = 0.5, spike_scale = 0.01, slab_scale = 1),
  lower = 0
)

Sample from prior predictive distribution

Description

Generates samples from the prior predictive distribution for the hierarchical shrinkage model. Useful for prior elicitation and sensitivity analysis.

The generative process is:

1. Sample mu from p(mu)

2. Sample tau from p(tau)

3. Sample theta_i ~ N(mu, tau) for each group i

Usage

sample_prior_predictive(
  hierarchical_priors,
  n_groups,
  n_draws = 1000,
  group_names = NULL
)

Arguments

Value

A list with class "shrinkr_prior_pred" containing:

See Also

shrink for fitting the hierarchical model, plot.shrinkr_prior_pred for visualizing prior predictive samples

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
median(prior_pred$implied_range)
head(as.data.frame(prior_pred))

Bayesian Hierarchical Shrinkage Model

Description

Applies hierarchical shrinkage to group-specific estimates using a two-stage Bayesian approach. Takes either a Gaussian mixture approximation of Stage 1 posteriors or point estimates with variance, and applies a Normal hierarchical model with flexible hyperpriors.

Usage

shrink(
  mixture = NULL,
  mle = NULL,
  var_matrix = NULL,
  hierarchical_priors = list(mu = distributional::dist_normal(0, 5), tau =
    distributional::dist_truncated(distributional::dist_normal(0, 2.5), lower = 0)),
  centered = FALSE,
  verbose = TRUE,
  ...
)

Arguments

Details

Model Specification

Hierarchical model (Stage 2):

\theta_g \mid \mu, \tau \sim \text{Normal}(\mu, \tau^2), \quad g = 1, \ldots, G

\mu \sim \pi(\mu)

\tau \sim \pi(\tau)

Stage 1 likelihood (Gaussian mixture approximation):

\theta_g \mid D_g \sim q_g(\theta_g) \approx \sum_{k=1}^K w_k \, \text{MVN}(\mu_k, \Sigma_k)

Full posterior:

\pi(\theta, \mu, \tau \mid D) \propto \left[\prod_{g=1}^G q_g(\theta_g)\right] \left[\prod_{g=1}^G \text{Normal}(\theta_g \mid \mu, \tau^2)\right] \pi(\mu) \pi(\tau)

where q_g(\theta_g) approximates the Stage 1 posterior for group g.

What's Fixed vs. Flexible

Fixed:

• Hierarchical distribution: \theta_g \mid \mu, \tau \sim \text{Normal}(\mu, \tau^2)

Flexible:

• Hyperpriors \pi(\mu) and \pi(\tau): Normal, Student-t, Cauchy, Lognormal, Gamma, Inverse-Gamma, Exponential, Uniform, mixtures (including spike-and-slab), and truncated versions

• Stage 1 posteriors: Can be non-Normal (handled by mixture approximation)

Critical Requirements

1. Stage 1 must use flat/uninformative priors on \theta_g

   • Ensures two-stage = one-stage hierarchical model

   • Stan: Don't specify prior (defaults to flat)

   • JAGS/NIMBLE: Use very wide priors

2. Verify mixture quality: plot(mixture, draws = samples)

   • Check density overlays and QQ plots

   • Poor approximation → biased shrinkage

3. Check prior implications: sample_prior_predictive(hierarchical_priors)

   • Understand what priors imply before fitting

   • Avoid prior-data conflicts

4. Minimum 2 groups required for heterogeneity estimation

Common Prior Choices for \tau

• Half-Normal: dist_truncated(dist_normal(0, s), lower = 0) - Weakly informative

• Half-t: dist_truncated(dist_student_t(df, 0, s), lower = 0) - Heavier tails

• Half-Cauchy: dist_truncated(dist_cauchy(0, s), lower = 0) - Very diffuse

• Uniform: dist_uniform(0, U) - Bounded heterogeneity

• Inverse-Gamma: dist_inverse_gamma(a, b) - Traditional choice

See vignette("getting_started") for complete workflow, vignette("brms_integration") for real examples, and package README for mathematical justification.

Value

A shrinkr_fit object (list) containing:

See Also

Workflow functions: fit_mixture(), sample_prior_predictive()

Extract results: extract_mu_tau(), extract_theta(), summarize_mu_tau(), summarize_theta(), theta_contrasts()

Visualization: plot.shrinkr_fit(), plot.shrinkr_mixture()

Vignettes:

• vignette("getting_started") - Complete workflow with Stan example

• vignette("brms_integration") - Survival analysis example

Examples

## Not run: 
# This example fits a Stan model, so it is not run during package checks.
priors <- list(
  mu = distributional::dist_normal(0, 10),
  tau = distributional::dist_truncated(distributional::dist_student_t(3, 0, 2.5), lower = 0)
)
fit <- shrink(
  mle = c(0.5, 1.2, -0.3),
  var_matrix = c(0.1, 0.15, 0.12),
  hierarchical_priors = priors,
  iter = 1000, chains = 2, seed = 1
)
summarise_theta(fit)

## End(Not run)

Imports from stats

Description

Imports from stats

Summarize mu and tau hyperparameters

Description

Computes posterior summaries for the hierarchical hyperparameters (mu, tau, and tau_squared). Returns a data frame with one row per parameter containing posterior means, standard deviations, quantiles, and convergence diagnostics.

This is a focused alternative to summary(fit), which returns summaries for all parameters including theta.

Usage

summarise_mu_tau(fit, probs = c(0.025, 0.5, 0.975), measures = NULL)

summarize_mu_tau(fit, probs = c(0.025, 0.5, 0.975), measures = NULL)

Arguments

Value

A data frame (tibble if available) with one row per parameter and columns:

parameter

Parameter name (mu, tau, or tau_squared)

mean

Posterior mean

sd

Posterior standard deviation

q2.5, q50, q97.5

Quantiles (or custom quantiles from probs)

rhat

R-hat convergence diagnostic

ess_bulk

Effective sample size (bulk)

ess_tail

Effective sample size (tail)

See Also

shrink() for fitting models, extract_mu_tau() for raw hyperparameter draws, summarise_theta() for group-level summaries

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
summarise_mu_tau(fit)
summarise_mu_tau(fit, probs = c(0.05, 0.5, 0.95))

Summarize theta parameters by group

Description

Computes posterior summaries for subgroup effects (theta parameters). Returns a data frame with one row per group containing posterior means, standard deviations, quantiles, and convergence diagnostics.

This is a focused alternative to summary(fit), which returns summaries for all parameters including mu and tau.

Usage

summarise_theta(
  fit,
  probs = c(0.025, 0.5, 0.975),
  group_names = NULL,
  measures = NULL
)

summarize_theta(
  fit,
  probs = c(0.025, 0.5, 0.975),
  group_names = NULL,
  measures = NULL
)

Arguments

Value

A data frame (tibble if available) with one row per group and columns:

group

Group identifier

mean

Posterior mean

sd

Posterior standard deviation

q2.5, q50, q97.5

Quantiles (or custom quantiles from probs)

rhat

R-hat convergence diagnostic

ess_bulk

Effective sample size (bulk)

ess_tail

Effective sample size (tail)

See Also

shrink() for fitting models, summarise_mu_tau() for hyperparameter summaries, theta_contrasts() for computing contrasts

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
summarise_theta(fit)
summarise_theta(fit, probs = c(0.05, 0.5, 0.95))
summarise_theta(fit, group_names = c("Control", "A", "B"))

Summary method for shrinkr_fit

Description

Comprehensive posterior summary including hyperparameters and all subgroup effects with convergence diagnostics.

Usage

## S3 method for class 'shrinkr_fit'
summary(
  object,
  probs = c(0.025, 0.5, 0.975),
  group_names = NULL,
  digits = 3,
  ...
)

Arguments

Value

Invisibly returns a list with mu_tau and theta summary tables. Prints formatted output.

See Also

shrink() for fitting models, summarise_theta() for theta-only summaries

Summary statistics for mixture fits

Description

Computes comprehensive summary statistics for fitted mixture components, including component-wise and variable-wise summaries.

Usage

## S3 method for class 'shrinkr_mixture'
summary(object, ...)

Arguments

Value

A list with class "summary.shrinkr_mixture" containing:

See Also

fit_mixture for fitting mixture models, print.shrinkr_mixture for quick overview

Examples

set.seed(1)
samples <- list(
  group1 = matrix(rnorm(100, 0.0, 0.5), ncol = 1),
  group2 = matrix(rnorm(100, 0.5, 0.5), ncol = 1)
)
mix <- fit_mixture(samples, K_max = 2, verbose = FALSE)
summ <- summary(mix)
summ$components

Summary statistics for prior predictive samples

Description

Computes comprehensive summary statistics for hyperparameters (mu and tau) and theta parameters from prior predictive samples.

Usage

## S3 method for class 'shrinkr_prior_pred'
summary(object, probs = c(0.025, 0.5, 0.975), ...)

Arguments

Value

A list with class "summary.shrinkr_prior_pred" containing:

Each data frame includes columns for mean, sd, and the requested quantiles.

See Also

sample_prior_predictive for generating samples, print.shrinkr_prior_pred for quick overview

Examples

priors <- list(
  mu = distributional::dist_normal(0, 5),
  tau = distributional::dist_truncated(distributional::dist_normal(0, 1), lower = 0)
)
prior_pred <- sample_prior_predictive(priors, n_groups = 3, n_draws = 50)
summ <- summary(prior_pred)
summ$theta

Linear combinations of theta

Description

Computes posterior draws for linear combinations of subgroup effects. Useful for pairwise contrasts (e.g., treatment vs control), weighted averages, or any custom linear estimand involving theta parameters.

Usage

theta_contrasts(fit, contrast_matrix, labels = NULL)

Arguments

Value

A posterior::draws_df with columns .chain, .iteration, .draw, and one column per contrast.

See Also

shrink() for fitting models, summarise_theta() for basic theta summaries

Examples

set.seed(1)
draws <- data.frame(
  mu = rnorm(20, 0.2, 0.05),
  tau = abs(rnorm(20, 0.3, 0.03)),
  `theta[1]` = rnorm(20, 0.0, 0.1),
  `theta[2]` = rnorm(20, 0.3, 0.1),
  `theta[3]` = rnorm(20, 0.5, 0.1),
  check.names = FALSE
)
draws$tau_squared <- draws$tau^2
fit <- list(
  fit = posterior::as_draws_df(draws),
  data = list(
    G = 3, K = 1, centered = FALSE,
    vars = c("group1", "group2", "group3"),
    quantiles = data.frame(
      q2.5 = c(-0.20, 0.10, 0.30),
      q50 = c(0.00, 0.30, 0.50),
      q97.5 = c(0.20, 0.50, 0.70)
    )
  ),
  summary = posterior::summarise_draws(
    posterior::as_draws_df(draws),
    "mean", "sd",
    ~posterior::quantile2(., probs = c(0.025, 0.5, 0.975))
  ),
  diagnostics = list(n_divergent = 0, max_treedepth = 0, n_leapfrog = 0)
)
class(fit) <- "shrinkr_fit"
L <- matrix(c(-1, 1, 0), nrow = 1)
contrast <- theta_contrasts(fit, L, labels = "group2_vs_group1")
posterior::summarise_draws(contrast)