ARD MCMC Sampler for Multiply-Imputed Regression

Description

Implements Bayesian variable selection using the Automatic Relevance Determination (ARD) prior across multiply-imputed datasets. The ARD prior imposes feature-specific shrinkage by placing a prior proportional to inverse of precision of each coefficient.

Usage

ARD_mcmc(
  X,
  Y,
  intercept = TRUE,
  nburn = 4000,
  npost = 4000,
  seed = NULL,
  verbose = TRUE,
  printevery = 1000,
  chain_index = 1
)

Arguments

Value

A named list with components:

post_beta

Array npost × D × p of sampled regression coefficients.

post_alpha

Matrix npost × D of sampled intercepts (if used).

post_sigma2

Numeric vector length npost, sampled residual variances.

post_psi2

Matrix npost × p of sampled precision parameters for each coefficient.

post_fitted_Y

Array npost × D × n of posterior predictive draws (with noise).

post_pool_beta

Matrix (npost * D) × p of pooled coefficient draws.

post_pool_fitted_Y

Matrix (npost * D) × n of pooled predictive draws (with noise).

hat_matrix_proj

Matrix D × n × n of averaged projection hat-matrices. To avoid recalculate for estimating degree of freedom.

Examples

sim <- sim_B(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5,
low_missing = TRUE, n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- ARD_mcmc(X, Y, nburn = 100, npost = 100)

Bayesian MI-LASSO for Multiply-Imputed Regression

Description

Fit a Bayesian multiple-imputation LASSO (BMI-LASSO) model across multiply-imputed datasets, using one of four priors: Multi-Laplace, Horseshoe, ARD, or Spike-Laplace. Automatically standardizes data, runs MCMC in parallel, performs variable selection via three-step projection predictive variable selection, and selects a final submodel by BIC.

Usage

BMI_LASSO(
  X,
  Y,
  model,
  standardize = TRUE,
  SNC = TRUE,
  grid = seq(0, 1, 0.01),
  orthogonal = FALSE,
  nburn = 4000,
  npost = 4000,
  seed = NULL,
  nchain = 1,
  ncores = 1,
  verbose = TRUE,
  printevery = 1000,
  ...
)

Arguments

Value

A named list with elements:

posterior

List of length nchain of MCMC outputs (posterior draws).

select

List of length nchain of logical matrices showing which variables are selected at each grid value.

best_select

List of length nchain of the single best selection (by BIC) for each chain.

posterior_best_models

List of length nchain of projected posterior draws for the best submodel.

bic_models

List of length nchain of BIC values and degrees-of-freedom for each candidate submodel.

summary_table_full

A data frame summarizing rank-normalized split-Rhat and other diagnostics for the full model.

summary_table_selected

A data frame summarizing diagnostics for the selected submodel after projection.

Examples

sim <- sim_A(n = 100, p = 20, type = "MAR", SNP = 1.5, low_missing = TRUE, n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- BMI_LASSO(X, Y, model = "Horseshoe",
                 nburn = 100, npost = 100,
                 nchain = 1, ncores = 1)
str(fit$best_select)

Multiple-Imputation LASSO (MI-LASSO)

Description

Fit a LASSO-like penalty across D multiply-imputed datasets by iteratively reweighted ridge regressions (Equation (4) of the manuscript). For each tuning parameter in lamvec, it returns the pooled coefficient estimates, the BIC, and the selected variables.

Usage

MI_LASSO(
  X,
  Y,
  lamvec = (2^(seq(-1, 4, by = 0.05)))^2/2,
  maxiter = 200,
  eps = 1e-20,
  ncores = 1
)

Arguments

Value

If length(lamvec) > 1, a list with elements:

best

List for the lambda with minimal BIC containing: coefficients ((p+1)×D intercept + slopes), bic (BIC scalar), varsel (logical length-p vector of selected predictors), lambda (the chosen penalty).

lambda_path

length(lamvec)×2 matrix of each lambda and its corresponding BIC.

If length(lamvec) == 1, returns a single list (as above) for that penalty.

Examples

sim <- sim_A(n = 100, p = 20, type = "MAR", SNP = 1.5, low_missing = TRUE, n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- MI_LASSO(X, Y, lamvec = c(0.1))

Horseshoe MCMC Sampler for Multiply-Imputed Regression

Description

Implements Bayesian variable selection using the hierarchical Horseshoe prior across multiply-imputed datasets. This model applies global–local shrinkage to regression coefficients via a global scale (tau2), local scales (lambda2), and auxiliary hyperpriors (kappa, eta).

Usage

horseshoe_mcmc(
  X,
  Y,
  intercept = TRUE,
  nburn = 4000,
  npost = 4000,
  seed = NULL,
  verbose = TRUE,
  printevery = 1000,
  chain_index = 1
)

Arguments

Value

A named list with components:

post_beta

Array npost × D × p of sampled regression coefficients.

post_alpha

Matrix npost × D of sampled intercepts (if used).

post_sigma2

Numeric vector of length npost, sampled residual variances.

post_lambda2

Matrix npost × p of local shrinkage parameters \lambda_j^2.

post_kappa

Matrix npost × p of auxiliary local hyperparameters \kappa_j.

post_tau2

Numeric vector of length npost, sampled global scale \tau^2.

post_eta

Numeric vector of length npost, sampled auxiliary global hyperparameter \eta.

post_fitted_Y

Array npost × D × n of posterior predictive draws (with noise).

post_pool_beta

Matrix (npost * D) × p of pooled coefficient draws.

post_pool_fitted_Y

Matrix (npost * D) × n of pooled predictive draws (with noise).

hat_matrix_proj

Matrix D × n × n of averaged projection hat-matrices. To avoid recalculate for estimating degree of freedom.

Examples

sim <- sim_B(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5,
low_missing = TRUE, n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- horseshoe_mcmc(X, Y, nburn = 100, npost = 100)

Multi-Laplace MCMC Sampler for Multiply-Imputed Regression

Description

Implements Bayesian variable selection under the Multi-Laplace prior on regression coefficients across multiply-imputed datasets. The prior shares local shrinkage parameters (lambda2) across imputations and places a Gamma(h, v) hyperprior on the global parameter rho.

Usage

multi_laplace_mcmc(
  X,
  Y,
  intercept = TRUE,
  h = 2,
  v = NULL,
  nburn = 4000,
  npost = 4000,
  seed = NULL,
  verbose = TRUE,
  printevery = 1000,
  chain_index = 1
)

Arguments

Value

A named list with elements:

post_beta

Array npost × D × p of sampled regression coefficients.

post_alpha

Matrix npost × D of sampled intercepts (if used).

post_sigma2

Numeric vector of length npost, sampled residual variances.

post_lambda2

Matrix npost × p of sampled local shrinkage parameters.

post_rho

Numeric vector of length npost, sampled global parameters.

post_fitted_Y

Array npost × D × n of posterior predictive draws (with noise).

post_pool_beta

Matrix (npost * D) × p of pooled coefficient draws.

post_pool_fitted_Y

Matrix (npost * D) × n of pooled predictive draws (with noise).

hat_matrix_proj

Matrix D × n × n of averaged projection hat-matrices. To avoid recalculate for estimating degree of freedom.

h, v

Numeric; the shape and scale hyperparameters used.

Examples

sim <- sim_B(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5, low_missing = TRUE,
n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- multi_laplace_mcmc(X, Y, intercept = TRUE, nburn = 100, npost = 100)

Projecting Posterior Means of Full-Model Coefficients onto a Reduced Subset Model

Description

Given posterior means of beta1_mat (and optional intercepts alpha1_vec) from a full model fitted on D imputed datasets, compute the predictive projection onto the submodel defined by xs_vec. Returns the projected coefficients (and intercepts, if requested).

Usage

projection_mean(X_arr, beta1_mat, xs_vec, sigma2, alpha1_vec = NULL)

Arguments

Value

A list with components:

beta2_mat

A D * p matrix of projected submodel coefficients.

alpha2_vec

(If alpha1_vec provided) numeric vector length D of projected intercepts.

Examples

# Simulate a single imputation with n=50, p=5:
D <- 3; n <- 50; p <- 5
X_arr <- array(rnorm(D * n * p), c(D, n, p))
beta1_mat <- matrix(rnorm(D * p), nrow = D)
# Suppose full-model sigma2 pooled is 1.2
sigma2 <- 1.2
# Project onto predictors 1 and 4 only:
xs_vec <- c(TRUE, FALSE, FALSE, TRUE, FALSE)
proj <- projection_mean(X_arr, beta1_mat, xs_vec, sigma2)
str(proj)

# With intercept:
alpha1_vec <- rnorm(D)
proj2 <- projection_mean(X_arr, beta1_mat, xs_vec, sigma2, alpha1_vec)
str(proj2)

Projection of Full-Posterior Draws onto a Reduced-Subset Model

Description

Given posterior draws beta1_arr (and optional intercepts alpha1_arr) from a full model fitted on D imputed datasets, compute the predictive projection of each draw onto the submodel defined by xs_vec. Returns the projected coefficients (and intercepts, if requested) plus the projected residual variance for each posterior draw.

Usage

projection_posterior(X_arr, beta1_arr, sigma1_vec, xs_vec, alpha1_arr = NULL)

Arguments

Value

A list with components:

beta2_arr

Array npost * D * p of projected submodel coefficients.

alpha2_arr

(If alpha1_arr provided) matrix npost * D of projected intercepts.

sigma2_opt

Numeric vector length npost of projected residual variances.

Examples

D <- 3; n <- 50; p <- 5; npost <- 100
X_arr      <- array(rnorm(D*n*p), c(D, n, p))
beta1_arr  <- array(rnorm(npost*D*p), c(npost, D, p))
sigma1_vec <- runif(npost, 0.5, 2)
xs_vec     <- c(TRUE, FALSE, TRUE, FALSE, TRUE)
# Without intercept
proj <- projection_posterior(X_arr, beta1_arr, sigma1_vec, xs_vec)
str(proj)
# With intercept draws
alpha1_arr <- matrix(rnorm(npost*D), nrow = npost, ncol = D)
proj2 <- projection_posterior(X_arr, beta1_arr, sigma1_vec, xs_vec, alpha1_arr)
str(proj2)

Simulate dataset A: Independent continuous covariates with MCAR/MAR missingness

Description

Generates a dataset for Scenario A used in Bayesian MI-LASSO benchmarking. Covariates are iid standard normal, with a fixed true coefficient vector, linear outcome, missingness imposed on specified columns under MCAR or MAR, and multiple imputations via predictive mean matching.

Usage

sim_A(
  n = 100,
  p = 20,
  type = "MAR",
  SNP = 1.5,
  low_missing = TRUE,
  n_imp = 5,
  seed = NULL
)

Arguments

Value

A list with components:

data_O

A list of complete covariate matrix and outcomes before missingness.

data_mis

A list of covariate matrix and outcomes with missing values.

data_MI

A list of array of imputed covariates (n_imp × n × p) and a matrix of imputed outcomes (n_imp × n).

data_CC

A list of complete-case covariate matrix and outcomes.

important

Logical vector of true nonzero coefficient indices.

covmat

True covariance matrix used for X.

beta

True coefficient vector.

Examples

sim <- sim_A(n = 100, p = 20, type = "MAR", SNP = 1.5,
             low_missing = TRUE, n_imp = 5, seed = 123)
str(sim)

Simulate dataset B: AR(1)-correlated continuous covariates with MCAR/MAR missingness

Description

Generates a dataset for Scenario B used in Bayesian MI-LASSO benchmarking. Covariates are multivariate normal with AR(1) covariance, with a fixed true coefficient vector, linear outcome, missingness imposed on specified columns under MCAR or MAR, and multiple imputations via predictive mean matching.

Usage

sim_B(
  n = 100,
  p = 20,
  low_missing = TRUE,
  type = "MAR",
  SNP = 1.5,
  corr = 0.5,
  n_imp = 5,
  seed = NULL
)

Arguments

Value

A list with components:

data_O

A list of complete covariate matrix and outcomes before missingness.

data_mis

A list of covariate matrix and outcomes with missing values.

data_MI

A list of array of imputed covariates (n_imp × n × p) and a matrix of imputed outcomes (n_imp × n).

data_CC

A list of complete-case covariate matrix and outcomes.

important

Logical vector of true nonzero coefficient indices.

covmat

True covariance matrix used for X.

beta

True coefficient vector.

Examples

sim <- sim_B(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5,
             low_missing = TRUE, n_imp = 5, seed = 123)
str(sim)

Simulate dataset C: AR(1)-latent Gaussian dichotomized to binary covariates with MCAR/MAR missingness

Description

Generates binary covariates by thresholding an AR(1) latent Gaussian, then proceeds as in sim_B.

Usage

sim_C(
  n = 100,
  p = 20,
  low_missing = TRUE,
  type = "MAR",
  SNP = 1.5,
  corr = 0.5,
  n_imp = 5,
  seed = NULL
)

Arguments

Value

A list with components:

data_O

A list of complete covariate matrix and outcomes before missingness.

data_mis

A list of covariate matrix and outcomes with missing values.

data_MI

A list of array of imputed covariates (n_imp × n × p) and a matrix of imputed outcomes (n_imp × n).

data_CC

A list of complete-case covariate matrix and outcomes.

important

Logical vector of true nonzero coefficient indices.

covmat

True covariance matrix used for X.

beta

True coefficient vector.

Examples

sim <- sim_C(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5,
             low_missing = TRUE, n_imp = 5, seed = 123)
str(sim)

Spike-and-Laplace MCMC Sampler for Multiply-Imputed Regression

Description

Implements Bayesian variable selection using a spike-and-slab prior with a Laplace (double-exponential) slab on nonzero coefficients. Latent inclusion indicators gamma follow Bernoulli(theta), and their probabilities follow independent Beta(a, b) priors.

Usage

spike_laplace_partially_mcmc(
  X,
  Y,
  intercept = TRUE,
  a = 2,
  b = NULL,
  nburn = 4000,
  npost = 4000,
  seed = NULL,
  verbose = TRUE,
  printevery = 1000,
  chain_index = 1
)

Arguments

Value

A named list with components:

post_rho

Numeric vector length npost, sampled global scale \rho.

post_gamma

Matrix npost * p of sampled inclusion indicators.

post_theta

Matrix npost * p of sampled Beta parameters \theta_j.

post_alpha

Matrix npost * D of sampled intercepts (if used).

post_lambda2

Matrix npost * p of sampled local scale parameters \lambda_j^2.

post_sigma2

Numeric vector length npost, sampled residual variances.

post_beta

Array npost * D * p of sampled regression coefficients.

post_fitted_Y

Array npost * D * n of posterior predictive draws (including noise).

post_pool_beta

Matrix (npost * D) * p of pooled coefficient draws.

post_pool_fitted_Y

Matrix (npost * D) * n of pooled predictive draws (with noise).

hat_matrix_proj

Matrix D * n * n of averaged projection hat-matrices. To avoid recalculate for estimating degree of freedom.

a, b

Numeric values of the rho hyperparameters used.

Examples

sim <- sim_B(n = 100, p = 20, type = "MAR", SNP = 1.5, corr = 0.5,
low_missing = TRUE, n_imp = 5, seed = 123)
X <- sim$data_MI$X
Y <- sim$data_MI$Y
fit <- spike_laplace_partially_mcmc(X, Y, nburn = 10, npost = 10)