SuSiE/mvSuSiE interface up to date with susieR and mvsusieR

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

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, ])

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"
   }

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))
+})

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, ])
+})

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