factorverse

Tools for Bayesian Factor Modeling via Rcpp

---

Installation

# Requires the devtools package
devtools::install_github("peterdunson/factorverse")

---

Overview

factorverse fits Bayesian factor models with a range of shrinkage and sparsity priors on the loading matrix, using Gibbs samplers implemented in C++ via RcppArmadillo. All methods share a single, consistent interface, so switching priors is a one-word change. The package also ships simulation and evaluation tools for studying covariance recovery, calibration, and computational efficiency against a known truth, and diagnostics for real-data fits.

Available methods: "mgps" (multiplicative gamma process shrinkage), "ssl" (spike-and-slab LASSO), "dl" (Dirichlet–Laplace), "hs" (horseshoe), "bass" (Bayesian structured sparsity), and "mnl" (mass-nonlocal score).

---

Quickstart

library(factorverse)

# 1. Build a block-sparse loading matrix and simulate data from it
Lambda <- sim_lambda(P = 12, K = 3, val = 1.2)
Sigma  <- diag(rep(0.3, 12))
Omega  <- tcrossprod(Lambda) + Sigma

dat <- sim_data(list(Lambda = Lambda, Sigma = Sigma, Omega = Omega),
                N = 300, seed = 42)

# 2. Fit a model (Dirichlet–Laplace prior)
fit <- fit_bfm(
  method = "dl", dat = dat, k0 = 3,
  chains = 1, nrun = 10000, burn = 5000,
  save_rds = FALSE, verbose = FALSE
)

# 3. Evaluate covariance recovery against the truth
omega_mse(fit, Omega)        # mean squared error of Omega-hat
omega_coverage(fit, Omega)   # 95% credible-interval coverage and width

To use your own data, replace dat with an n × P numeric matrix (rows are observations, columns variables) and drop the truth-based evaluation in step 3, relying instead on the posterior summaries from the fit.

---

Fitting Models

• fit_bfm(method, dat, k0, ...) Fit a Bayesian factor model with the chosen prior. The main user-facing interface; forwards to run_rcpp_model.

  fit <- fit_bfm(
    method = "ssl", dat = Y, k0 = 5,
    chains = 4, nrun = 10000, burn = 5000
  )

• run_rcpp_model(method, dat, ..., chains = 4, nrun = 16000, burn = 8000) The underlying sampler driver. Equivalent to fit_bfm; exposes the full set of arguments (factor count, adaptation, post-processing thresholds, output options).

The returned object contains the posterior draws in $posterior (the loadings Lambda, scores b, residual variances Sigma, and per-iteration factor count) and the posterior-mean loadings in $Lambda_hat.

---

Simulation & Evaluation

Simulation designs use block-sparse loadings: each factor loads on a disjoint block of P/K variables (K = 1 reduces to a single dense loading vector). The design and evaluation tools are built around this structure, letting you reproduce the package's own simulation studies at any (K, P, N), or build a custom design.

• sim_lambda(P, K, val) Build a P x K block-sparse loading matrix with loading value val.

  Lambda <- sim_lambda(P = 10, K = 2, val = 5)

• sim_params(P, K, lambdas = c(1, 5), sigma_equal, sigma_unequal) Build a list of parameter sets crossing loading magnitudes with residual variance structures (defaults reproduce the package's four-set design), each containing Lambda, Sigma, and the implied Omega.

  param_sets <- sim_params(P = 10, K = 2)
  param_sets[[1]]$Omega

• sim_data(params, N, seed = NULL) Generate one dataset from a parameter set (an element of sim_params(), or any list with Lambda, Sigma, Omega).

  dat <- sim_data(param_sets[[1]], N = 200, seed = 1)

• sim_fit(dat, method, k0, ...) Fit one dataset with one method, returning point estimates, credible intervals for Omega, wall-clock fit time, and (optionally) effective sample size (ESS). Automatically switches to a memory-light streaming path at high dimension (P > 50 by default) so large runs don't require holding a full P x P x ndraws array in memory.

  fit_info <- sim_fit(dat, method = "dl", k0 = 2, compute_ess = TRUE)

• sim_run(P, K, N = 200, B = 1, methods, ...) The full replicate-loop driver: builds the parameter sets, generates B datasets per set, fits every method to every dataset, and returns per-replicate recovery, coverage, and efficiency metrics. Reproduces any scenario from the package's simulation studies by setting K, P, N, and B; saves one .rds file per method as it completes if out_dir is supplied, so progress survives an interruption. Timing is recorded by default (record_time = TRUE); ESS is computed on the first ess_reps replicates of each parameter set (ess_reps = 20 by default, set to 0 to disable).

  study <- sim_run(
    P = 10, K = 2, N = 200, B = 100,
    methods = c("dl", "mgps"),
    nrun = 10000, burn = 5000
  )

• sim_summarize(results) Collapse the per-replicate output of sim_run() into the three summary tables used throughout the package's simulation results: mse (mean/median/SD), coverage (empirical coverage and interval width), and efficiency (mean fit time, mean ESS, and ESS per second).

  tables <- sim_summarize(study)
  tables$mse
  tables$coverage
  tables$efficiency

• omega_mse(fits, Omega_true, method = c("posterior_mean", "plugin")) Mean squared error of the estimated covariance against the truth, for one already-fitted model or a list of fits. Reports error over all entries and split by diagonal/off-diagonal.

  omega_mse(fit, Omega)
  omega_mse(list(DL = fit1, MGPS = fit2), Omega)

• omega_coverage(fits, Omega_true, level = 0.95) Empirical coverage and mean width of posterior credible intervals for the entries of Omega, split by diagonal (variances) and off-diagonal (covariances), for one already-fitted model or a list of fits.

  omega_coverage(fit, Omega)

omega_mse()/omega_coverage() evaluate a model you've already fit, useful for a one-off check against a validation truth. sim_run()/sim_summarize() handle the full replicated-study workflow (many datasets, many methods, one call) and are what the package's own simulation studies use.

---

Real-Data Diagnostics

• analyze_factor_fits(Y, fits, fit_labels = NULL, n_perm = 1000, plot = TRUE, side_by_side = FALSE) Compare residual correlation structure across multiple fits against a permutation null.

  analyze_factor_fits(Y, fits = list(SSL = fit1, HS = fit2))

• score_correlation(Y, fit, label = "Model") Summarize residual–residual, factor–residual, and factor–factor correlations for one fit.

  score_correlation(Y, fit, label = "SSL Model")

• bfm_convergence(fit) Compute R̂ and bulk/tail effective sample sizes for the loadings, residual variances, and linear predictors.

  conv_df <- bfm_convergence(fit)

• lambda_by_sample(Y, fit) Assess loading sparsity and covariance-reconstruction error across posterior draws.

  lambda_stats <- lambda_by_sample(Y, fit)

• residuals_by_sample(Y, fit, n_perm = 1000) Per-draw residual diagnostics: Shapiro–Wilk p-values, Wasserstein distances to the null, skewness, and kurtosis.

  res_stats <- residuals_by_sample(Y, fit, n_perm = 500)

• summarize_fit_diagnostics(data_matrix, fits, n_perm = 1000, prob_seq = seq(0.1, 0.9, 0.1), null_sd = 0.02) Aggregate quantile-MSE and normality metrics across fits.

  diag_df <- summarize_fit_diagnostics(Y, fits = list(SSL = fit1, HS = fit2))

• summarize_normality(mat) Column-wise Shapiro–Wilk tests; returns min and median p-values.

  norm_stats <- summarize_normality(Y)

---

Documentation & Vignettes

help(package = "factorverse")
browseVignettes("factorverse")

---

License

This project is licensed under the MIT License.