From c9d78d9fb038e04428003aa47f9e4d1e2e1e607d Mon Sep 17 00:00:00 2001 From: Gao Wang Date: Fri, 24 Apr 2026 17:30:37 -0400 Subject: [PATCH 1/2] Bring SuSiE/mvSuSiE interface up to date with susieR and mvsusieR susieR: rename sketch LD parameter from stochastic_ld_sample to sketch_samples in susie_rss_wrapper. mvsusieR: drop hard-coded option overrides that no longer exist on the mvsusieR refactor (precompute_covariances -> precompute_cache with default TRUE; compute_objective, approximate, verbosity removed). Rely on mvsusie defaults for estimate_residual_variance / estimate_prior_variance / estimate_prior_mixture_weights / estimate_prior_method. In the post-processor, read the renamed outcome_names field, compute coef via coef.mvsusie(), and drop b1_rescaled (trivially derivable from alpha, mu, and X_column_scale_factors; documented inline). Co-Authored-By: Claude Opus 4.7 (1M context) --- R/multivariate_pipeline.R | 5 ++--- R/regularized_regression.R | 10 +--------- R/susie_wrapper.R | 21 ++++++++++++--------- 3 files changed, 15 insertions(+), 21 deletions(-) diff --git a/R/multivariate_pipeline.R b/R/multivariate_pipeline.R index 67f49715..2a664441 100644 --- a/R/multivariate_pipeline.R +++ b/R/multivariate_pipeline.R @@ -286,9 +286,8 @@ multivariate_analysis_pipeline <- function( message("Fitting mvSuSiE model on input data ...") res$mvsusie_fitted <- mvsusieR::mvsusie(X, Y, L = max_L, prior_variance = mvsusie_reweighted_mixture_prior$data_driven_prior_matrices, - residual_variance = resid_Y, precompute_covariances = FALSE, compute_objective = TRUE, - estimate_residual_variance = FALSE, estimate_prior_variance = TRUE, estimate_prior_method = "EM", - max_iter = mvsusie_max_iter, n_thread = 1, approximate = FALSE, verbosity = verbose, coverage = coverage[1] + residual_variance = resid_Y, + max_iter = mvsusie_max_iter, coverage = coverage[1] ) # Process mvSuSiE results diff --git a/R/regularized_regression.R b/R/regularized_regression.R index 27d4b980..a003d4e6 100644 --- a/R/regularized_regression.R +++ b/R/regularized_regression.R @@ -324,18 +324,10 @@ mvsusie_weights <- function(mvsusie_fit = NULL, X = NULL, Y = NULL, prior_varian stop("Both X and Y must be provided if mvsusie_fit is NULL.") } if (is.null(prior_variance)) prior_variance <- mvsusieR::create_mixture_prior(R = ncol(Y)) - if (is.null(residual_variance)) { - if (!requireNamespace("mr.mashr", quietly = TRUE)) { - stop("Package 'mr.mashr' is required for residual variance estimation. Install with: devtools::install_github('stephenslab/mr.mashr')") - } - residual_variance <- mr.mashr:::compute_cov_flash(Y) - } mvsusie_fit <- mvsusieR::mvsusie( X = X, Y = Y, L = L, prior_variance = prior_variance, - residual_variance = residual_variance, precompute_covariances = FALSE, - compute_objective = TRUE, estimate_residual_variance = FALSE, estimate_prior_variance = TRUE, - estimate_prior_method = "EM", approximate = FALSE, ... + residual_variance = residual_variance, ... ) } return(mvsusieR::coef.mvsusie(mvsusie_fit)[-1, ]) diff --git a/R/susie_wrapper.R b/R/susie_wrapper.R index 0bd5b4ca..8af0fa7b 100644 --- a/R/susie_wrapper.R +++ b/R/susie_wrapper.R @@ -180,7 +180,7 @@ susie_rss_wrapper <- function(z, R = NULL, X = NULL, n = NULL, # Build argument list for susie_rss base_args <- list(z = z, n = n, coverage = coverage, - stochastic_ld_sample = sketch_samples, ...) + sketch_samples = sketch_samples, ...) if (!is.null(X)) base_args$X <- X else base_args$R <- R run_with_L <- function(L_val) do.call(susie_rss, c(base_args, list(L = L_val))) @@ -354,7 +354,7 @@ susie_post_processor <- function(susie_output, data_x, data_y, X_scalar, y_scala if (analysis_script != "") res$analysis_script <- analysis_script if (!is.null(other_quantities)) res$other_quantities <- other_quantities if (mode == "mvsusie") { - res$context_names <- susie_output$condition_names + res$context_names <- susie_output$outcome_names } if (!is.null(data_y)) { # Mode-specific processing @@ -453,14 +453,17 @@ susie_post_processor <- function(susie_output, data_x, data_y, X_scalar, y_scala res$susie_result_trimmed$mu2 <- susie_output$mu2[eff_idx, , drop = FALSE] } if (mode == "mvsusie") { - # res$susie_result_trimmed$b1 = susie_output$b1[eff_idx, , , drop = FALSE] - # res$susie_result_trimmed$b2 = susie_output$b2[eff_idx, , , drop = FALSE] - res$susie_result_trimmed$b1_rescaled <- susie_output$b1_rescaled[eff_idx, , , drop = FALSE] - res$susie_result_trimmed$coef <- susie_output$coef + res$susie_result_trimmed$coef <- mvsusieR::coef.mvsusie(susie_output) res$susie_result_trimmed$clfsr <- susie_output$conditional_lfsr[eff_idx, , , drop = FALSE] - # other lfsr can be computed: - # se_lfsr <- mvsusie_single_effect_lfsr(clfsr, alpha) - # lfsr <- mvsusie_get_lfsr(clfsr, alpha) + # b1_rescaled (per-effect posterior mean on the original X scale, L x J x R or L x J) + # is not stored: it is trivially derivable from alpha, mu, and X_column_scale_factors. + # To reconstruct from a fit `s`: + # csd <- s$X_column_scale_factors + # if (length(dim(s$mu)) == 3) { + # b1_rescaled <- sweep(sweep(s$mu, c(1,2), s$alpha, "*"), 2, csd, "/") + # } else { + # b1_rescaled <- sweep(s$alpha * s$mu, 2, csd, "/") + # } } class(res$susie_result_trimmed) <- "susie" } From 76486a691db531e72e623eb98ce1a0c1a27dae7a Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Fri, 24 Apr 2026 17:43:15 -0700 Subject: [PATCH 2/2] update tests --- tests/testthat/test_multivariate_pipeline.R | 62 +++++++++++++++++++++ tests/testthat/test_rr_mrmash_mvsusie.R | 44 +++++++++++++++ tests/testthat/test_susie_wrapper.R | 50 +++++++++++++++++ 3 files changed, 156 insertions(+) diff --git a/tests/testthat/test_multivariate_pipeline.R b/tests/testthat/test_multivariate_pipeline.R index f08e711a..fc630197 100644 --- a/tests/testthat/test_multivariate_pipeline.R +++ b/tests/testthat/test_multivariate_pipeline.R @@ -470,3 +470,65 @@ test_that("initialize_mvsusie_prior runs with provided data_driven_prior_matrice expect_true("reweighted_mixture_prior_cv" %in% names(result)) expect_true("mvsusie_fitted" %in% names(result)) }) + +test_that("pipeline propagates outcome_names from mvsusie through susie_post_processor", { + skip_if(!requireNamespace("mvsusieR", quietly = TRUE), + "mvsusieR not installed") + skip_if(!requireNamespace("susieR", quietly = TRUE), + "susieR not installed") + d <- make_mv_data() + r <- ncol(d$Y) + p <- ncol(d$X) + L <- 5 + vnames <- colnames(d$X) + cnames <- colnames(d$Y) + fake_coef <- matrix(rnorm((p + 1) * r), nrow = p + 1, ncol = r) + + prior_U <- list( + udd_1 = matrix(0.5, r, r) + 0.5 * diag(r), + udd_2 = diag(r) + ) + for (k in seq_along(prior_U)) { + rownames(prior_U[[k]]) <- colnames(prior_U[[k]]) <- cnames + } + prior_mats <- list(U = prior_U, w = c(udd_1 = 0.5, udd_2 = 0.5)) + fake_w0 <- c("udd_1_a" = 0.4, "udd_2_a" = 0.4, "null" = 0.2) + + local_mocked_bindings( + mrmash_wrapper = function(X, Y, ...) { + list(V = diag(ncol(Y)), w0 = fake_w0, + w1 = matrix(0.1, nrow = ncol(X), ncol = 1)) + }, + ) + # Mock mvsusie to return outcome_names (the new field name) + local_mocked_bindings( + mvsusie = function(...) { + list( + pip = setNames(rep(0.1, p), vnames), + alpha = matrix(1 / p, nrow = L, ncol = p), + lbf_variable = matrix(0, nrow = L, ncol = p), + V = rep(1, L), + sets = list(cs = NULL, requested_coverage = 0.95), + niter = 10, + outcome_names = cnames, + conditional_lfsr = array(0.5, dim = c(L, p, r)) + ) + }, + create_mixture_prior = function(...) list(matrices = prior_U, weights = c(0.5, 0.5)), + coef.mvsusie = function(...) fake_coef, + .package = "mvsusieR" + ) + + result <- multivariate_analysis_pipeline( + X = d$X, Y = d$Y, maf = d$maf, + pip_cutoff_to_skip = 0, + data_driven_prior_matrices = prior_mats, + twas_weights = FALSE + ) + expect_true(is.list(result)) + # outcome_names should propagate as context_names + expect_equal(result$context_names, cnames) + # coef and clfsr should be present in susie_result_trimmed + expect_equal(result$susie_result_trimmed$coef, fake_coef) + expect_equal(dim(result$susie_result_trimmed$clfsr), c(L, p, r)) +}) diff --git a/tests/testthat/test_rr_mrmash_mvsusie.R b/tests/testthat/test_rr_mrmash_mvsusie.R index af746d8a..a0c32105 100644 --- a/tests/testthat/test_rr_mrmash_mvsusie.R +++ b/tests/testthat/test_rr_mrmash_mvsusie.R @@ -35,3 +35,47 @@ test_that("mvsusie_weights errors when X and Y are NULL and fit is NULL", { expect_error(mvsusie_weights(mvsusie_fit = NULL, X = NULL, Y = NULL), "Both X and Y must be provided") }) + +test_that("mvsusie_weights fits model and returns coefficients when fit is NULL", { + skip_if_not(requireNamespace("mvsusieR", quietly = TRUE), + "mvsusieR not installed") + set.seed(42) + n <- 30 + p <- 5 + R <- 3 + X <- matrix(rnorm(n * p), n, p) + Y <- matrix(rnorm(n * R), n, R) + fake_coef <- matrix(rnorm((p + 1) * R), nrow = p + 1, ncol = R) + + local_mocked_bindings( + create_mixture_prior = function(...) list(), + mvsusie = function(...) "mock_fit", + coef.mvsusie = function(...) fake_coef, + .package = "mvsusieR" + ) + + result <- expect_message( + mvsusie_weights(X = X, Y = Y), + "mvsusie_fit is not provided" + ) + # Should return coef without intercept row + expect_equal(dim(result), c(p, R)) + expect_equal(result, fake_coef[-1, ]) +}) + +test_that("mvsusie_weights returns coefficients from provided fit", { + skip_if_not(requireNamespace("mvsusieR", quietly = TRUE), + "mvsusieR not installed") + p <- 5 + R <- 3 + fake_coef <- matrix(rnorm((p + 1) * R), nrow = p + 1, ncol = R) + + local_mocked_bindings( + coef.mvsusie = function(...) fake_coef, + .package = "mvsusieR" + ) + + result <- mvsusie_weights(mvsusie_fit = "precomputed_fit") + expect_equal(dim(result), c(p, R)) + expect_equal(result, fake_coef[-1, ]) +}) diff --git a/tests/testthat/test_susie_wrapper.R b/tests/testthat/test_susie_wrapper.R index d162c956..b3b84d94 100644 --- a/tests/testthat/test_susie_wrapper.R +++ b/tests/testthat/test_susie_wrapper.R @@ -602,6 +602,56 @@ test_that("susie_post_processor uses 1:max_L for eff_idx when V is NULL (fSuSiE) expect_null(result$susie_result_trimmed$V) }) +# ============================================================================= +# susie_post_processor: mvsusie mode (outcome_names, coef, clfsr) +# ============================================================================= + +test_that("susie_post_processor stores outcome_names, coef, and clfsr in mvsusie mode", { + skip_if_not_installed("susieR") + skip_if_not_installed("mvsusieR") + p <- 5 + L <- 3 + R <- 2 + vnames <- paste0("chr1:", 1:p, ":A:G") + cnames <- paste0("cond_", 1:R) + fake_coef <- matrix(rnorm((p + 1) * R), nrow = p + 1, ncol = R) + + fake_output <- list( + pip = setNames(rep(0.01, p), vnames), + alpha = matrix(1 / p, nrow = L, ncol = p), + lbf_variable = matrix(0, nrow = L, ncol = p), + sets = list(cs = NULL, requested_coverage = 0.95), + niter = 10, + V = rep(1, L), + outcome_names = cnames, + conditional_lfsr = array(0.5, dim = c(L, p, R)) + ) + + n <- 20 + X <- matrix(rnorm(n * p), n, p) + colnames(X) <- vnames + + local_mocked_bindings( + coef.mvsusie = function(...) fake_coef, + .package = "mvsusieR" + ) + + result <- susie_post_processor( + fake_output, + data_x = X, + data_y = NULL, + X_scalar = 1, y_scalar = 1, + mode = "mvsusie" + ) + + # outcome_names should be stored as context_names + expect_equal(result$context_names, cnames) + # coef should come from mvsusieR::coef.mvsusie + expect_equal(result$susie_result_trimmed$coef, fake_coef) + # conditional_lfsr should be trimmed to eff_idx + expect_equal(dim(result$susie_result_trimmed$clfsr), c(L, p, R)) +}) + # ============================================================================= # susie_wrapper: dynamic-L break when cs is NULL (Tier 2) # =============================================================================