From e6e9a9a2c558fae7687047dbff6d07cd555977b0 Mon Sep 17 00:00:00 2001 From: Daniel Nachun Date: Mon, 12 Oct 2026 12:28:05 -0500 Subject: [PATCH] fix fsusie issues --- NAMESPACE | 5 +- R/MethodOptions.R | 30 ++ R/fineMappingPipeline.R | 22 +- R/fineMappingWrappers.R | 416 +++++++++--------- R/jointEngine.R | 14 +- R/regularizedRegressionWrappers.R | 159 +++++-- R/twasWeightsPipeline.R | 24 +- _pkgdown.yml | 5 +- man/MvsusieMixturePriorOptions.Rd | 34 ++ man/TwasWeightsMethodsParam.Rd | 2 +- man/buildMvsusieMixturePrior.Rd | 64 +++ man/fitFsusie.Rd | 25 +- man/fitGlmnetCv.Rd | 72 +++ man/fitMvsusie.Rd | 9 +- man/fitSusieInfThenSusieRss.Rd | 60 --- man/fsusieWeights.Rd | 6 +- man/fsusieWrapper.Rd | 69 --- man/mvsusieWeights.Rd | 2 +- man/twasEngineOptions.Rd | 14 +- tests/testthat/test_MethodOptions.R | 1 + tests/testthat/test_MethodParam.R | 6 +- tests/testthat/test_fineMappingPipeline.R | 65 ++- tests/testthat/test_fineMappingRow.R | 3 +- tests/testthat/test_fineMappingWrappers.R | 267 +++++------ .../test_regularizedRegressionWrappers.R | 147 ++++++- tests/testthat/test_twasWeightsPipeline.R | 37 +- 26 files changed, 988 insertions(+), 570 deletions(-) create mode 100644 man/MvsusieMixturePriorOptions.Rd create mode 100644 man/buildMvsusieMixturePrior.Rd create mode 100644 man/fitGlmnetCv.Rd delete mode 100644 man/fitSusieInfThenSusieRss.Rd delete mode 100644 man/fsusieWrapper.Rd diff --git a/NAMESPACE b/NAMESPACE index 378c0196..47c4572d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -95,6 +95,7 @@ export(MrmashOptions) export(MrmashPriorParam) export(MrmashRssOptions) export(MultiStudyQtlDataset) +export(MvsusieMixturePriorOptions) export(MvsusieOptions) export(MvsusieRssOptions) export(NcvregOptions) @@ -157,6 +158,7 @@ export(blockMetadata) export(buildLdEigen) export(buildLdScore) export(buildMrmashPriorMatrices) +export(buildMvsusieMixturePrior) export(buildTopLoci) export(calculateFeatureScores) export(causalInferencePipeline) @@ -220,11 +222,11 @@ export(fineMappingResult) export(fineMappingRow) export(finemapCtwasRegions) export(fitFsusie) +export(fitGlmnetCv) export(fitMashContrast) export(fitMvsusie) export(fitMvsusieRss) export(fitRetention) -export(fitSusieInfThenSusieRss) export(flattenTupleRanges) export(formatFinemappingOutput) export(frqData) @@ -232,7 +234,6 @@ export(fsusieAffectedRegions) export(fsusieCredibleBand) export(fsusieGetCs) export(fsusieWeights) -export(fsusieWrapper) export(fullFit) export(genomicRegion) export(genotypeCovariates) diff --git a/R/MethodOptions.R b/R/MethodOptions.R index 44811cd6..98607cf6 100644 --- a/R/MethodOptions.R +++ b/R/MethodOptions.R @@ -390,6 +390,36 @@ setClass("MethodOptions", contains = "SimpleList") )) } +# Refuse names for a reason OTHER than pecotmr supplying them -- a name the +# engine does not have, or one it has that pecotmr will not expose. Those +# cannot go through .configRefuseOwned(), whose message promises the value +# would be overwritten and names where to set it instead; here there is no +# such place, so the reason is stated outright. `unavailable` is a named +# character vector, name = the refused argument, value = why, read as +# "`` cannot be set here: ". +# +# This exists for engines taking `...`, where an unknown name is accepted and +# silently ignored instead of reported (see .engineAcceptedNames). Everywhere +# else the engine's own formals are the check. +# @noRd +.configRefuseUnavailable <- function(extra, unavailable, label) { + clash <- intersect(names(extra), names(unavailable)) + if (length(clash) == 0L) { + return(invisible(NULL)) + } + why <- str_flatten( + sprintf( + "`%s` cannot be set here: %s", + clash, + unname( + unavailable[clash] + ) + ), + collapse = "; " + ) + abort(glue("{label}: {why}.")) +} + # TRUE for a record a constructor produced. Used by the pipelines to refuse a # bare list, which is the whole point of requiring constructors. # @noRd diff --git a/R/fineMappingPipeline.R b/R/fineMappingPipeline.R index d35b1022..0a762fde 100644 --- a/R/fineMappingPipeline.R +++ b/R/fineMappingPipeline.R @@ -1920,7 +1920,7 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { return(.fmFoldWeightsMv(Xtr, Ytr, coverage, userArgs, mvPrior)) } if (token == "fsusie") { - return(.fmFoldWeightsFsusie(Xtr, Ytr, pos, userArgs)) + return(.fmFoldWeightsFsusie(Xtr, Ytr, pos, coverage, userArgs)) } NULL } @@ -1971,9 +1971,11 @@ combineFineMappingResults <- function(..., ldSketch = NULL) { # Per-fold fsusie weights. # @noRd -.fmFoldWeightsFsusie <- function(Xtr, Ytr, pos, userArgs) { +.fmFoldWeightsFsusie <- function(Xtr, Ytr, pos, coverage, userArgs) { + # coverage reaches the fold fit here for the same reason it does in the + # susie and mvsusie fold paths; L is left to susiF's default there too. fsArgs <- .fmMergeUserArgs( - list(X = Xtr, Y = Ytr, pos = pos), + list(X = Xtr, Y = Ytr, pos = pos, covLev = coverage), "fsusie", userArgs ) @@ -3501,14 +3503,18 @@ setMethod( ), csHome["coverage"] ), + # Only what the pipeline ACTUALLY supplies. `prior`, `max_scale` and + # `max_SNP_EM` used to be listed here because this row described the + # retired fsusieWrapper(), which passed them; the pipeline calls + # fitFsusie() and never did, so refusing them closed the only path to + # setting them at all -- the same trap the `residual_variance` note + # above records. "fsusieR::susiF" = c( - set_names(rep(fromData, 4L), c("X", "Y", "pos", "prior")), + set_names(rep(fromData, 3L), c("X", "Y", "pos")), c( L = "`maxNumSingleEffects()` on CredibleSetParam", - min_purity = "the wrapper's `minPurity`", - max_SNP_EM = "the wrapper's `maxSnpEm`", - cov_lev = "the wrapper's `covLev`", - max_scale = "the wrapper's `maxScale`" + cov_lev = "`csCoverage()` on CredibleSetParam", + min_purity = "`minAbsCorr()` on CredibleSetParam" ) ), character(0) diff --git a/R/fineMappingWrappers.R b/R/fineMappingWrappers.R index 5fef0a8e..3912b1b4 100644 --- a/R/fineMappingWrappers.R +++ b/R/fineMappingWrappers.R @@ -267,81 +267,6 @@ fitSusieInfThenSusie <- function( list(susie = susieFit, susieInf = susieInfFit) } -#' Two-stage SuSiE-RSS Fine-mapping -#' -#' RSS analog of \code{fitSusieInfThenSusie}. Fits SuSiE-inf via \code{susieRss} -#' first, then initialises standard SuSiE-RSS from the SuSiE-inf result. The -#' single pair of fits can be used both for fine-mapping post-processing and -#' TWAS weight extraction. -#' -#' @param z Numeric vector of z-scores. -#' @param R LD correlation matrix. -#' @param n Sample size (scalar). -#' @param args Defaults forwarded to both fits, built with -#' \code{\link{SusieOptions}}. A bare list is refused, since it cannot be -#' checked; \code{list()} (the default) means no options. -#' @param SusieInfOptions SuSiE-inf-specific overrides, built with -#' \code{\link{SusieInfOptions}}; they take precedence over \code{args}. -#' @param SusieOptions Standard SuSiE-RSS-specific overrides, built with -#' \code{\link{SusieOptions}}; they take precedence over \code{args}. -#' @param fittedModels Optional list with pre-fitted \code{$susie} and/or -#' \code{$susieInf} objects to skip re-fitting. -#' @return A list with \code{susie} and \code{susieInf} fit objects. -#' @importFrom susieR susie_rss -#' @examples -#' data(eqtlRegionExample) -#' X <- eqtlRegionExample$X[, 1:30] -#' y <- eqtlRegionExample$yRes -#' ss <- lapply(seq_len(ncol(X)), function(j) { -#' coef(summary(lm(y ~ X[, j])))[2, 1:2] -#' }) -#' stat <- list( -#' bhat = vapply(ss, `[`, numeric(1), 1L), -#' shat = vapply(ss, `[`, numeric(1), 2L), -#' z = vapply(ss, function(s) s[1] / s[2], numeric(1)), -#' n = rep(nrow(X), ncol(X))) -#' LD <- cor(X) -#' fitSusieInfThenSusieRss(z = stat$z, R = LD, n = nrow(X)) -#' @importFrom checkmate assertNumeric assertNumber assertList -#' @export -fitSusieInfThenSusieRss <- function( - z, - R, - n, - args = list(), - SusieInfOptions = list(), - SusieOptions = list(), - fittedModels = NULL -) { - assertNumeric(z) - assertNumber(n, lower = 0, finite = TRUE) - .fmAssertChainArgs(args, SusieInfOptions, SusieOptions) - assertList(fittedModels, null.ok = TRUE) - # RSS analog of fitSusieInfThenSusie, built from the shared per-token RSS - # fitter (.fmFitSusieRss). .fmFitSusieRss tags every fit "susieRss", so the - # inf fit is re-tagged "susieInf" to preserve this wrapper's contract. - cached <- fittedModels %||% list() - infRaw <- cached[["susieInf"]] %||% - .fmFitSusieRss( - z, - R, - n, - "susieInf", - userArgs = .fmChainStageArgs(args, SusieInfOptions) - ) - susieInfFit <- .setFinemappingFitClass(infRaw, "susieInf") - susieRaw <- cached[["susie"]] %||% - .fmFitSusieRss( - z, - R, - n, - "susie", - chainFromInf = susieInfFit, - userArgs = .fmChainStageArgs(args, SusieOptions) - ) - susieFit <- .setFinemappingFitClass(susieRaw, "susieRss") - list(susie = susieFit, susieInf = susieInfFit) -} #' Post-process Fine-mapping Fits #' @@ -2357,8 +2282,13 @@ trimFinemappingFit <- function(fit, effectIdx, method, csTables) { } keep <- intersect( c( - "coef", "fitted_func", "cred_band", "outing_grid", "cs", - "csd_X", "lfsr_func" + "coef", + "fitted_func", + "cred_band", + "outing_grid", + "cs", + "csd_X", + "lfsr_func" ), names(populated) ) @@ -2575,99 +2505,6 @@ fsusieGetCs <- function(fsusieObj, X, requestedCoverage = 0.95) { return(sets) } -#' @title Wrapper for fsusie Function with Automatic Post-Processing -#' -#' @description This function serves as a wrapper for the fsusie function, -#' facilitating automatic post-processing such as removing dummy credible sets -#' (cs) that don't meet the minimum purity threshold and calculating -#' correlations for the remaining cs. The function parameters are identical to -#' those of the fSuSiE function. -#' -#' @param X Residual genotype matrix. -#' @param Y Response phenotype matrix. -#' @param pos Genomics position of phenotypes, used for specifying the wavelet -#' model. -#' @param maxNumSingleEffects The maximum number of the credible set. -#' @param prior method to generate the prior. -#' @param maxSnpEm maximum number of SNP used for learning the prior. -#' @param covLev Coverage level for the credible sets. -#' @param maxScale numeric, define the maximum of wavelet coefficients used in -#' the analysis (2^maxScale). Set 10 true by default. -#' @param minPurity Minimum purity threshold for credible sets to be retained. -#' @param methodArgs Options forwarded to \code{fsusieR::susiF}, built with -#' \code{\link{FsusieOptions}}. A bare list is refused: it cannot be checked -#' against the engine, so a misspelled option would be silently ignored. -#' @return A modified fsusie object with the susie sets list, correlations for -#' cs, alpha as df like susie, and without the dummy cs that do not meet the -#' minimum purity requirement. -#' @examples -#' data(eqtlRegionExample) -#' X <- eqtlRegionExample$X[, 1:50] -#' n <- nrow(X) -#' nPos <- 16 -#' base <- sin(seq(0, 2 * pi, length.out = nPos)) -#' Y <- matrix(rep(base, each = n), n, nPos) + -#' X[, 1] %o% (0.5 * cos(seq(0, pi, length.out = nPos))) -#' pos <- seq_len(nPos) -#' fsusieWrapper(X, Y, pos = pos, maxNumSingleEffects = 2, -#' prior = "mixture_normal_per_scale", maxSnpEm = 10, -#' covLev = 0.95, minPurity = 0.5, maxScale = 6) -#' @export -fsusieWrapper <- function( - X, - Y, - pos, - maxNumSingleEffects, - prior, - maxSnpEm, - covLev, - minPurity, - maxScale, - methodArgs = FsusieOptions() -) { - .assertMethodOptions(methodArgs, "FsusieOptions", "methodArgs") - if (!requireNamespace("fsusieR", quietly = TRUE)) { - abort("Package 'fsusieR' is required for this function.") - } - callArgs <- list_modify( - list( - X = X, - Y = Y, - pos = pos, - L = maxNumSingleEffects, - prior = prior, - max_SNP_EM = maxSnpEm, - cov_lev = covLev, - min_purity = minPurity, - max_scale = maxScale - ), - !!!methodArgs - ) - fsusieObj <- exec(fsusieR::susiF, !!!callArgs) - .fsusieWrapperPostprocess(fsusieObj, X, minPurity, covLev) -} - -# Drop dummy credible sets below the purity threshold (else build sets + CS -# correlations), then reshape alpha (per-effect list) into a single data.frame. -.fsusieWrapperPostprocess <- function(fsusieObj, X, minPurity, covLev) { - withSets <- if (all(abs(as.numeric(fsusieObj$purity)) < minPurity)) { - list_assign( - fsusieObj, - cs = list(NULL), - sets = list(cs = list(NULL), requested_coverage = covLev) - ) - } else { - list_assign( - fsusieObj, - sets = fsusieGetCs(fsusieObj, X, requestedCoverage = covLev) - ) - } - list_assign( - withSets, - alpha = bind_rows(map(fsusieObj$alpha, .asDataFrameT)) - ) -} - # ============================================================================= # Uniform fit wrappers for mvSuSiE (individual + RSS) @@ -2686,9 +2523,10 @@ fsusieWrapper <- function( #' @param X Numeric matrix of genotypes (samples x variants). #' @param Y Numeric matrix of multi-trait / multi-context outcomes (samples x #' conditions). -#' @param prior_variance Prior variance matrix; pass the output of -#' \code{mvsusieR::create_mixture_prior(R = ncol(Y))} unless you have a -#' domain-specific prior. +#' @param prior_variance Prior variance; build it with +#' \code{\link{buildMvsusieMixturePrior}}, which is the single interface +#' to \code{mvsusieR::create_mixture_prior()} and defaults to the +#' canonical mixture over the conditions it is given. #' @param coverage Credible set coverage (default 0.95). #' @param methodArgs Options forwarded to \code{mvsusieR::mvsusie}, built #' with \code{\link{MvsusieOptions}}. A bare list is refused: it cannot be @@ -2704,7 +2542,7 @@ fsusieWrapper <- function( #' X <- multiTraitData$X[, 1:60] #' Y <- multiTraitData$Y #' fitMvsusie(X, Y, -#' prior_variance = mvsusieR::create_mixture_prior(R = ncol(Y))) +#' prior_variance = buildMvsusieMixturePrior(conditionNames = colnames(Y))) #' } #' @importFrom checkmate assertNumber #' @export @@ -2785,6 +2623,153 @@ fitMvsusieRss <- function( exec(mvsusieR::mvsusie_rss, !!!callArgs) } +#' @title Arguments For mvsusieR's Mixture-Prior Builder +#' @description Arguments for \code{mvsusieR::create_mixture_prior()}, the +#' function that builds the mvSuSiE prior mixture. The mixture itself, the +#' condition count, the conditions to restrict to and the weight cutoff are +#' all pecotmr's to supply, so they are refused here --- they have formals +#' of their own on \code{\link{buildMvsusieMixturePrior}}. +#' +#' \code{create_mixture_prior()} takes \code{...}, which it forwards to +#' \code{create_cov_canonical()}, so names given here cannot be checked +#' against its formals: a misspelling reaches the engine, which ignores it +#' silently. See \code{\link{MethodOptions}}. +#' +#' No defaults: the engine's own are left in place. +#' @param ... Arguments for \code{mvsusieR::create_mixture_prior()} --- +#' \code{null_weight}, \code{max_mixture_len} --- or the +#' \code{create_cov_canonical()} passthroughs \code{singletons} and +#' \code{hetgrid} that shape the canonical mixture. +#' @return A \code{\link{MethodOptions}} object. +#' @examples +#' MvsusieMixturePriorOptions(singletons = FALSE) +#' @export +MvsusieMixturePriorOptions <- function(...) { + extra <- list(...) + .configRefuseOwned( + extra, + c( + mixture_prior = str_c( + "`mixtureMatrices` / `mixtureWeights` on ", + "buildMvsusieMixturePrior()" + ), + R = "`R` on buildMvsusieMixturePrior()", + include_indices = "`conditionNames` on buildMvsusieMixturePrior()", + weights_tol = "`weightsTol` on buildMvsusieMixturePrior()" + ), + "MvsusieMixturePriorOptions" + ) + .newMethodOptions( + "mvsusieR::create_mixture_prior", + defaults = list(), + extra = extra, + label = "MvsusieMixturePriorOptions", + engine = "mvsusieMixturePrior" + ) +} + +# The `exactly one of mixture_prior / R` arm of create_mixture_prior(). The +# engine tests those two with missing(), so only ONE of the names may appear +# in the call at all -- passing the unused one as NULL still counts as +# supplied and trips its "Require exactly one of" stop. +# @noRd +.mvMixtureArm <- function(R, conditionNames, mixtureMatrices, mixtureWeights) { + if (!is.null(mixtureMatrices)) { + mixture <- compact(list( + matrices = mixtureMatrices, + weights = mixtureWeights + )) + return(list(mixture_prior = mixture)) + } + # The canonical arm labels R x R covariances with `conditionNames` and + # then subsets by those very names, so a count that disagrees with the + # names fails deep inside the engine on a dimnames length. Checking `R` + # here also catches the default standing in for unnamed conditions: + # length(NULL) is 0, which the engine turns into a 0 x 0 subscript error. + assertCount(R, positive = TRUE) + if (!is.null(conditionNames) && length(conditionNames) != R) { + abort(glue( + "buildMvsusieMixturePrior(): `conditionNames` names ", + "{length(conditionNames)} conditions but `R` is {R}. The ", + "canonical mixture labels R x R covariances with those names, ", + "so the two must agree." + )) + } + list(R = R) +} + +#' Build the mvSuSiE mixture prior +#' +#' The single interface to \code{mvsusieR::create_mixture_prior()}, which +#' every mvSuSiE prior in the package goes through. That function takes +#' \strong{exactly one} of \code{mixture_prior} and \code{R} and tests the +#' two with \code{missing()}, so the unused one cannot simply be passed as +#' \code{NULL}; this wrapper owns that branch. Supply +#' \code{mixtureMatrices} for a data-driven mixture --- e.g. the reweighted +#' mr.mash covariances --- or leave it NULL for the canonical uniform +#' mixture over \code{R} conditions. +#' +#' @param conditionNames Character vector naming the conditions the prior is +#' over, forwarded as \code{include_indices}. On the data-driven arm this +#' \emph{restricts} the supplied covariances to these conditions by name, +#' so a mixture spanning more conditions than the fit at hand is subset +#' here. On the canonical arm it only labels the generated covariances. +#' NULL leaves the mixture unrestricted and unlabelled. +#' @param R Number of conditions for the canonical mixture. Defaults to +#' \code{length(conditionNames)}; pass it explicitly when the conditions +#' are unnamed. Unused when \code{mixtureMatrices} is supplied, since the +#' mixture then carries its own dimension. +#' @param mixtureMatrices Named list of prior covariance matrices to use as +#' the mixture. NULL (default) builds the canonical mixture from \code{R} +#' instead. +#' @param mixtureWeights Numeric vector of mixture weights, named to match +#' \code{mixtureMatrices}. NULL leaves the engine's uniform weights. +#' @param weightsTol Mixture components weighted at or below this are +#' dropped. +#' @param methodArgs Options forwarded to +#' \code{mvsusieR::create_mixture_prior}, built with +#' \code{\link{MvsusieMixturePriorOptions}}. A bare list is refused: it +#' cannot be checked against the engine. +#' @return The prior object \code{mvsusieR::create_mixture_prior()} returns, +#' ready to pass as \code{prior_variance} to \code{\link{fitMvsusie}} or +#' \code{\link{fitMvsusieRss}}. +#' @examples +#' \donttest{ +#' # Canonical uniform mixture over two named conditions. +#' buildMvsusieMixturePrior(conditionNames = c("brain", "blood")) +#' } +#' @importFrom checkmate assertCount +#' @export +buildMvsusieMixturePrior <- function( + conditionNames = NULL, + R = length(conditionNames), + mixtureMatrices = NULL, + mixtureWeights = NULL, + weightsTol = 1e-10, + methodArgs = MvsusieMixturePriorOptions() +) { + if (!requireNamespace("mvsusieR", quietly = TRUE)) { + abort("Package 'mvsusieR' is required.") + } + .assertMethodOptions( + methodArgs, + "MvsusieMixturePriorOptions", + "methodArgs" + ) + arm <- .mvMixtureArm(R, conditionNames, mixtureMatrices, mixtureWeights) + callArgs <- list_modify( + c( + arm, + compact(list( + include_indices = conditionNames, + weights_tol = weightsTol + )) + ), + !!!methodArgs + ) + exec(mvsusieR::create_mixture_prior, !!!callArgs) +} + #' Fit fSuSiE on individual-level (X, Y, pos) data #' #' Thin wrapper around \code{fsusieR::susiF}. @@ -2792,6 +2777,18 @@ fitMvsusieRss <- function( #' @param X Numeric matrix of genotypes (samples x variants). #' @param Y Numeric matrix of multi-trait outcomes (samples x traits). #' @param pos Numeric vector of trait positions (length \code{ncol(Y)}). +#' @param maxNumSingleEffects Maximum number of single effects +#' (\code{susiF}'s \code{L}). \code{NULL} leaves susiF's own default. An +#' explicit formal rather than an \code{\link{FsusieOptions}} entry because +#' the pipeline supplies it from \code{maxNumSingleEffects()} on +#' \code{\link{CredibleSetParam}}, the same way \code{\link{fitMvsusieRss}} +#' takes \code{coverage}. +#' @param covLev Credible-set coverage level (\code{susiF}'s \code{cov_lev}). +#' \code{NULL} leaves susiF's own default; the pipeline supplies +#' \code{csCoverage()} on \code{\link{CredibleSetParam}}. +#' @param minPurity Minimum credible-set purity (\code{susiF}'s +#' \code{min_purity}). \code{NULL} leaves susiF's own default; the pipeline +#' supplies \code{minAbsCorr()} on \code{\link{CredibleSetParam}}. #' @param methodArgs Options forwarded to \code{fsusieR::susiF}, built with #' \code{\link{FsusieOptions}}. A bare list is refused: it cannot be checked #' against the engine, so a misspelled option would be silently ignored. @@ -2807,9 +2804,29 @@ fitMvsusieRss <- function( #' pos <- seq_len(nPos) #' fitFsusie(X, Y, pos = pos, methodArgs = FsusieOptions(maxit = 5)) #' @export -fitFsusie <- function(X, Y, pos, methodArgs = FsusieOptions()) { +fitFsusie <- function( + X, + Y, + pos, + maxNumSingleEffects = NULL, + covLev = NULL, + minPurity = NULL, + methodArgs = FsusieOptions() +) { .assertMethodOptions(methodArgs, "FsusieOptions", "methodArgs") - callArgs <- list_modify(list(X = X, Y = Y, pos = pos), !!!methodArgs) + # These three are refused by FsusieOptions() as pipeline-owned, so they + # travel as formals: .splitMethodArgs() routes a formal straight here + # instead of through the constructor that would reject it. compact() + # keeps susiF's own defaults when a caller leaves one NULL. + owned <- compact(list( + L = maxNumSingleEffects, + cov_lev = covLev, + min_purity = minPurity + )) + callArgs <- list_modify( + c(list(X = X, Y = Y, pos = pos), owned), + !!!methodArgs + ) exec(fsusieR::susiF, !!!callArgs) } @@ -3236,7 +3253,7 @@ susieAshRssWeights <- function( #' X <- multiTraitData$X[, 1:60] #' Y <- multiTraitData$Y #' fit <- fitMvsusie(X = X, Y = Y, -#' prior_variance = mvsusieR::create_mixture_prior(R = ncol(Y)), +#' prior_variance = buildMvsusieMixturePrior(R = ncol(Y)), #' methodArgs = MvsusieOptions(max_iter = 50)) #' mvsusieWeights(mvsusieFit = fit) #' } @@ -3338,9 +3355,9 @@ mvsusieWeights <- function(mvsusieFit = NULL) { #' Y <- matrix(rep(base, each = n), n, nPos) + #' X[, 1] %o% (0.5 * cos(seq(0, pi, length.out = nPos))) #' pos <- seq_len(nPos) -#' fit <- fsusieWrapper(X, Y, pos = pos, maxNumSingleEffects = 2, -#' prior = "mixture_normal_per_scale", maxSnpEm = 10, -#' covLev = 0.95, minPurity = 0.5, maxScale = 6) +#' fit <- fitFsusie(X, Y, pos = pos, maxNumSingleEffects = 2, +#' methodArgs = FsusieOptions(prior = "mixture_normal_per_scale", +#' max_SNP_EM = 10, max_scale = 6)) #' fsusieWeights(fsusieFit = fit, X = X, Y = Y, #' variantIds = colnames(X)) #' @export @@ -3437,7 +3454,7 @@ fsusieWeights <- function( } # Normalize alpha to a list of per-effect vectors. fsusieR::susiF returns a -# list; fsusieWrapper reshaping yields an L x nSNP matrix/data.frame. +# list; a caller that has reshaped it hands over an L x nSNP matrix/data.frame. .fsusieAlphaList <- function(alpha) { if (is.list(alpha) && !is.data.frame(alpha)) { return(map(alpha, as.numeric)) @@ -5127,14 +5144,14 @@ extractTopPipInfo <- function(fmRow, sumstats) { } # The canonical (non-reweighted) mvSuSiE mixture prior for residual variance -# `V`: -# create_mixture_prior(R) restricted to the group's conditions. +# `V`: the uniform canonical mixture over the group's conditions, paired with +# `V` in the shape the mvsusie callers consume. # @noRd .fmCanonicalPrior <- function(V, conditionNames, R) { list( - priorVariance = mvsusieR::create_mixture_prior( - R = R, - include_indices = conditionNames + priorVariance = buildMvsusieMixturePrior( + conditionNames = conditionNames, + R = R ), residualVariance = V ) @@ -5143,16 +5160,17 @@ extractTopPipInfo <- function(fmRow, sumstats) { # Rebuild the mvSuSiE data-driven *reweighted* mixture prior + residual variance # from a stored mr.mash fit -- the lean payload # (list(dataDrivenPriorMatrices, w0, V)) that mrmashWeights() attaches at -# fitRetention "slim" and twasWeightsPipeline keeps on the mrmash row. Shared -# by the fine-mapping mvsusie consumer and the twas mvsusie_weights consumer. +# fitRetention "slim" and twasWeightsPipeline keeps on the mrmash row. The +# fine-mapping mvsusie paths are the only consumers: mvsusieWeights() extracts +# from a supplied fit and no longer fits, so it needs no prior of its own. # # Reproduces the deleted multivariate_pipeline.R reweighting bit-identically: # rescaleCovW0(w0) collapses the expanded mr.mash weights onto the original # data-driven covariance matrices ($U), filters to surviving components, and -# create_mixture_prior() wraps them, restricted to the fit's conditions +# buildMvsusieMixturePrior() wraps them, restricted to the fit's conditions # (`conditionNames` = colnames(Y)). `V` becomes mvsusie's residual_variance. # A NULL fit, NULL matrices, or no surviving component falls back to the -# canonical create_mixture_prior(R), matching the legacy `else` branch. +# canonical mixture, matching the legacy `else` branch. # Returns list(priorVariance, residualVariance) (residualVariance NULL only # when no fit was supplied at all). # @noRd @@ -5183,12 +5201,12 @@ extractTopPipInfo <- function(fmRow, sumstats) { if (length(w0Updated) == 0L) { return(.fmCanonicalPrior(fitParts$V, conditionNames, R)) } - mixture <- list(matrices = ddpm$U[names(w0Updated)], weights = w0Updated) list( - priorVariance = mvsusieR::create_mixture_prior( - mixture_prior = mixture, - weights_tol = weightsTol, - include_indices = conditionNames + priorVariance = buildMvsusieMixturePrior( + conditionNames = conditionNames, + mixtureMatrices = ddpm$U[names(w0Updated)], + mixtureWeights = w0Updated, + weightsTol = weightsTol ), residualVariance = fitParts$V ) @@ -5264,7 +5282,8 @@ extractTopPipInfo <- function(fmRow, sumstats) { # Per-fold mvsusie weights. Reuses the data-driven reweighted prior + residual # covariance from the full-data mr.mash fit on every fold -- the prior is over # conditions, identical across folds (only samples are held out). NULL mvPrior -# -> canonical prior (unchanged behavior). +# -> the canonical mixture over this fold's outcome columns, which is the same +# prior .fmCanonicalPrior() builds for the full-data path. # @noRd .fmFoldWeightsMv <- function(Xtr, Ytr, coverage, userArgs, mvPrior) { baseArgs <- list( @@ -5272,7 +5291,10 @@ extractTopPipInfo <- function(fmRow, sumstats) { Y = Ytr, coverage = coverage, prior_variance = if (is.null(mvPrior)) { - mvsusieR::create_mixture_prior(R = ncol(Ytr)) + buildMvsusieMixturePrior( + conditionNames = colnames(Ytr), + R = ncol(Ytr) + ) } else { mvPrior$priorVariance } diff --git a/R/jointEngine.R b/R/jointEngine.R index faa6337e..2a589a22 100644 --- a/R/jointEngine.R +++ b/R/jointEngine.R @@ -288,8 +288,20 @@ setMethod( abort(msg) } verbose <- if (is.null(cfg$verbose)) 1 else cfg$verbose + # Supply what .fmOwnedEngineArgs("fsusieR::susiF") says is supplied: + # these are refused in FsusieOptions, so the pipeline is their only + # source. They are fitFsusie() formals, so .splitMethodArgs() routes them + # past the constructor that refuses them. + csp <- cfg$credibleSetParam fitArgs <- .fmMergeUserArgs( - list(X = Xc, Y = Yc, pos = .jgTraitPos(group)), + list( + X = Xc, + Y = Yc, + pos = .jgTraitPos(group), + maxNumSingleEffects = csp$maxNumSingleEffects, + covLev = csp$coverage, + minPurity = csp$minAbsCorr + ), "fsusie", args$methodArgs[["fsusie"]] ) diff --git a/R/regularizedRegressionWrappers.R b/R/regularizedRegressionWrappers.R index 6018c771..4c880976 100644 --- a/R/regularizedRegressionWrappers.R +++ b/R/regularizedRegressionWrappers.R @@ -778,10 +778,22 @@ initPriorSd <- function(X, y, n = 30) { #' itself is refused --- \code{GlmnetOptions(alpha=)} and #' \code{NcvregOptions(penalty=)} select the penalty family, which is the #' calling wrapper's to set. See \code{\link{MethodOptions}}. +#' +#' \code{GlmnetOptions} rejects two names by hand, which the \code{...} +#' check above would otherwise let through in silence. \code{nfold} is not +#' a \code{cv.glmnet} argument at all --- the fold count is \code{nfolds}, +#' and the shorter spelling only ever worked through R's partial matching. +#' \code{parallel} is a real one, but pecotmr does not expose it: glmnet +#' parallelises over its CV folds through a \code{foreach} backend that +#' this package neither declares nor registers, and weight fitting already +#' runs inside \code{BiocParallel} workers over the same folds, so +#' enabling it would nest two frameworks and oversubscribe the machine. +#' The glasso and univariate-glmnet \code{mr.mash} initialisations are +#' single-threaded for the same reason. #' @param ... Arguments for the engine, under its own names. #' @return A \code{\link{MethodOptions}} object. #' @examples -#' GlmnetOptions(nfold = 10) +#' GlmnetOptions(nfolds = 10) #' QggOptions(updateB = TRUE) #' LassosumOptions(thr = 1e-4) #' @name twasEngineOptions @@ -791,6 +803,28 @@ NULL #' @export GlmnetOptions <- function(...) { extra <- list(...) + # cv.glmnet takes `...`, so neither of these would be reported by the + # name check: the first is not an argument at all, the second is one + # pecotmr will not expose. Both are called out by hand instead. + .configRefuseUnavailable( + extra, + c( + nfold = str_c( + "it is not an argument of glmnet::cv.glmnet(), whose fold ", + "count is `nfolds`. pecotmr passed this misspelling itself ", + "until recently, so it reached the engine only through R's ", + "partial matching" + ), + parallel = str_c( + "glmnet parallelises over its CV folds through a foreach ", + "backend, which pecotmr neither declares nor registers. ", + "Weight fitting already runs inside BiocParallel workers ", + "over the same folds, so a second framework underneath it ", + "would oversubscribe the machine rather than speed it up" + ) + ), + "GlmnetOptions" + ) .configRefuseOwned( extra, c( @@ -1154,6 +1188,82 @@ MrmashRssOptions <- function(...) { ) } +#' Fit a cross-validated glmnet model +#' +#' The single interface to \code{glmnet::cv.glmnet}, which every glmnet fit in +#' the package goes through. It returns the \emph{fit}, not weights: +#' \code{\link{glmnetWeights}} extracts the \code{lambda.min} coefficient +#' vector from it, while the \code{mr.mash} initialisation paths need the +#' coefficient matrix, the intercept row, the lambda sequence and +#' \code{predict()} on held-out rows --- the same fit-versus-extract split as +#' \code{\link{fitMvsusie}} and \code{\link{mvsusieWeights}}. +#' +#' The four knobs below are formals rather than \code{\link{GlmnetOptions}} +#' entries because they are what the package's glmnet call sites genuinely +#' differ on, and leaving them implicit is how three different conventions +#' came to coexist. Each defaults to \code{NULL}, meaning the argument is not +#' passed at all and \code{cv.glmnet}'s own default stands; a value supplied +#' through \code{methodArgs} still wins over the caller's, so none of them is +#' refused by the constructor. +#' +#' @param x Numeric design matrix (samples x features). +#' @param y Response: a vector, or a matrix for \code{family = "mgaussian"}. +#' @param alpha Elastic-net mixing parameter in \code{[0, 1]}: \code{1} = +#' lasso, \code{0} = ridge, \code{0.5} = elastic net. The penalty family is +#' the calling wrapper's to choose, which is why \code{GlmnetOptions()} +#' refuses it. +#' @param family glmnet's response family, e.g. \code{"gaussian"} or +#' \code{"mgaussian"} for a matrix response. \code{NULL} leaves glmnet's +#' own default. +#' @param standardize Logical; standardize the columns of \code{x} before +#' fitting. \code{NULL} leaves glmnet's own default (\code{TRUE}). +#' @param intercept Logical; fit an intercept. \code{NULL} leaves glmnet's +#' own default (\code{TRUE}). +#' @param nfolds Number of cross-validation folds. \code{NULL} leaves +#' \code{cv.glmnet}'s own default (10). +#' @param methodArgs Further arguments forwarded to +#' \code{glmnet::cv.glmnet}, built with \code{\link{GlmnetOptions}}. A bare +#' list is refused, since it cannot be checked against the engine. +#' @return The fit object returned by \code{glmnet::cv.glmnet}. +#' @examples +#' data(eqtlRegionExample) +#' X <- eqtlRegionExample$X[, 1:30] +#' y <- eqtlRegionExample$yRes +#' fit <- fitGlmnetCv(X, y, alpha = 0.5, nfolds = 5) +#' coef(fit, s = "lambda.min")[1:5] +#' @export +fitGlmnetCv <- function( + x, + y, + alpha, + family = NULL, + standardize = NULL, + intercept = NULL, + nfolds = NULL, + methodArgs = GlmnetOptions() +) { + if (!requireNamespace("glmnet", quietly = TRUE)) { + abort("Package 'glmnet' is required for this function.") + } + .assertMethodOptions(methodArgs, "GlmnetOptions", "methodArgs") + # compact() is what makes a NULL formal mean "do not pass it": an + # explicit NULL would reach cv.glmnet and override its default with + # nothing. methodArgs is spliced last, so a user value wins over the + # caller's -- these are literal defaults, not pecotmr-owned arguments. + callArgs <- list_modify( + c( + list(x = x, y = y, alpha = alpha), + compact(list( + family = family, + standardize = standardize, + intercept = intercept, + nfolds = nfolds + )) + ), + !!!methodArgs + ) + exec(glmnet::cv.glmnet, !!!callArgs) +} #' Compute TWAS weights via penalized regression (glmnet) #' #' Fit a cross-validated elastic-net / lasso model with \code{glmnet::cv.glmnet} @@ -1175,28 +1285,21 @@ MrmashRssOptions <- function(...) { #' glmnetWeights(X, y, alpha = 0.5) #' @export glmnetWeights <- function(X, y, alpha, methodArgs = GlmnetOptions()) { - # Check if glmnet is installed - if (!requireNamespace("glmnet", quietly = TRUE)) { - abort("Package 'glmnet' is required for this function.") - } - .assertMethodOptions(methodArgs, "GlmnetOptions", "methodArgs") eff.wgt <- matrix(0, ncol = 1, nrow = ncol(X)) keep <- .dropZeroVariance(X, "glmnetWeights") # `alpha` is the wrapper's own dial and is refused by GlmnetOptions() - # itself, so there is no second check here. Everything else pecotmr - # fixed below is a literal default and stays overridable. - callArgs <- list_modify( - list( - x = X[, keep, drop = FALSE], - y = y, - alpha = alpha, - nfold = 5, - intercept = TRUE, - standardize = FALSE - ), - !!!methodArgs + # itself, so there is no second check here. The three settings below are + # literal defaults and stay overridable through `methodArgs`; the glmnet + # requirement check and the constructor check live in fitGlmnetCv(). + enet <- fitGlmnetCv( + x = X[, keep, drop = FALSE], + y = y, + alpha = alpha, + standardize = FALSE, + intercept = TRUE, + nfolds = 5, + methodArgs = methodArgs ) - enet <- exec(glmnet::cv.glmnet, !!!callArgs) replace(eff.wgt, keep, coef(enet, s = "lambda.min")[2:(sum(keep) + 1)]) } @@ -3687,19 +3790,18 @@ computeCoefficientsGlasso <- function( ) { assertFlag(standardize) assertInt(numThreads) - n <- nrow(X) p <- ncol(X) r <- ncol(Y) conditionNames <- colnames(Y) - # Fit group-lasso - cvfitGlmnet <- glmnet::cv.glmnet( + # Fit group-lasso. `parallel = FALSE` used to be restated here; it is + # cv.glmnet's own default, so transcribing it bought nothing. + cvfitGlmnet <- fitGlmnetCv( x = X, y = Y, - family = "mgaussian", alpha = 1, - standardize = standardize, - parallel = FALSE + family = "mgaussian", + standardize = standardize ) coeffGlmnet <- coef(cvfitGlmnet, s = "lambda.min") @@ -3733,13 +3835,12 @@ computeCoefficientsGlasso <- function( Ynomiss <- Y[samplesKept, i, drop = FALSE] Xnomiss <- X[samplesKept, , drop = FALSE] - cvfit <- glmnet::cv.glmnet( + cvfit <- fitGlmnetCv( x = Xnomiss, y = Ynomiss, - family = "gaussian", alpha = alpha, - standardize = standardize, - parallel = FALSE + family = "gaussian", + standardize = standardize ) coeffic <- as.vector(coef(cvfit, s = "lambda.min")) lambdaSeq <- cvfit$lambda diff --git a/R/twasWeightsPipeline.R b/R/twasWeightsPipeline.R index 80a87bab..811f027e 100644 --- a/R/twasWeightsPipeline.R +++ b/R/twasWeightsPipeline.R @@ -3737,21 +3737,21 @@ setMethod( # @return Normalized coefficient vector of length Kvalid. # @noRd .solveEnsembleGlmnet <- function(Pvalid, yObs, Kvalid, alpha = 1) { - # Engine-call exception: glmnet is a constrained-regression solver here, - # fitting the ensemble's stacking weights -- a different role from its use - # as a weight-learning engine in regularizedRegressionWrappers.R. The - # stacking step is pipeline logic, so the call stays with it. - if (!requireNamespace("glmnet", quietly = TRUE)) { - abort("Package 'glmnet' is required for solver='glmnet'.") - } - + # Here glmnet is a constrained-regression solver fitting the ensemble's + # stacking weights -- a different role from weight learning, but the same + # engine, so the call goes through fitGlmnetCv() like every other. Only + # the stacking logic is pipeline-side. `lower.limits = 0` is what makes + # the regression non-negative; it has no pecotmr formal, so it travels in + # the engine's own options bag. `standardize` is left at glmnet's default + # because the columns of `Pvalid` are already on one scale: they are all + # predictions of the same outcome. fit <- try_fetch( - glmnet::cv.glmnet( + fitGlmnetCv( x = Pvalid, y = yObs, - lower.limits = 0, alpha = alpha, - intercept = FALSE + intercept = FALSE, + methodArgs = GlmnetOptions(lower.limits = 0) ), error = function(cnd) { msg <- glue( @@ -4628,7 +4628,7 @@ setClass("TwasWeightsMethodsParam", contains = "MethodsSelectionParam") #' TwasWeightsMethodsParam(methods = list(lasso = list(), susie = list())) #' TwasWeightsMethodsParam( #' methods = list(susie = list()), -#' qtlDatasetMethods = list(lasso = GlmnetOptions(nfold = 5)), +#' qtlDatasetMethods = list(lasso = GlmnetOptions(nfolds = 5)), #' qtlSumStatsMethods = list(lasso = LassosumOptions()) #' ) #' @export diff --git a/_pkgdown.yml b/_pkgdown.yml index 319d0883..edfb8b7a 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -214,6 +214,7 @@ reference: - MashPosteriorOptions - MrmashCanonicalCovsOptions - MrmashExpandCovsOptions + - MvsusieMixturePriorOptions - PlinkQcOptions - PoolrOptions - QvalueOptions @@ -608,10 +609,9 @@ reference: - fsusieWeights - fitMvsusie - fitMvsusieRss + - buildMvsusieMixturePrior - fitFsusie - - fitSusieInfThenSusieRss - fsusieGetCs - - fsusieWrapper - susieResult - computeCsTables - buildTopLoci @@ -635,6 +635,7 @@ reference: - subtitle: "Penalized regression" contents: + - fitGlmnetCv - glmnetWeights - lassoWeights - enetWeights diff --git a/man/MvsusieMixturePriorOptions.Rd b/man/MvsusieMixturePriorOptions.Rd new file mode 100644 index 00000000..cd9018e9 --- /dev/null +++ b/man/MvsusieMixturePriorOptions.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fineMappingWrappers.R +\name{MvsusieMixturePriorOptions} +\alias{MvsusieMixturePriorOptions} +\title{Arguments For mvsusieR's Mixture-Prior Builder} +\usage{ +MvsusieMixturePriorOptions(...) +} +\arguments{ +\item{...}{Arguments for \code{mvsusieR::create_mixture_prior()} --- +\code{null_weight}, \code{max_mixture_len} --- or the +\code{create_cov_canonical()} passthroughs \code{singletons} and +\code{hetgrid} that shape the canonical mixture.} +} +\value{ +A \code{\link{MethodOptions}} object. +} +\description{ +Arguments for \code{mvsusieR::create_mixture_prior()}, the + function that builds the mvSuSiE prior mixture. The mixture itself, the + condition count, the conditions to restrict to and the weight cutoff are + all pecotmr's to supply, so they are refused here --- they have formals + of their own on \code{\link{buildMvsusieMixturePrior}}. + + \code{create_mixture_prior()} takes \code{...}, which it forwards to + \code{create_cov_canonical()}, so names given here cannot be checked + against its formals: a misspelling reaches the engine, which ignores it + silently. See \code{\link{MethodOptions}}. + + No defaults: the engine's own are left in place. +} +\examples{ +MvsusieMixturePriorOptions(singletons = FALSE) +} diff --git a/man/TwasWeightsMethodsParam.Rd b/man/TwasWeightsMethodsParam.Rd index c389395a..d9128274 100644 --- a/man/TwasWeightsMethodsParam.Rd +++ b/man/TwasWeightsMethodsParam.Rd @@ -46,7 +46,7 @@ Selects the methods \code{\link{twasWeightsPipeline}} runs and TwasWeightsMethodsParam(methods = list(lasso = list(), susie = list())) TwasWeightsMethodsParam( methods = list(susie = list()), - qtlDatasetMethods = list(lasso = GlmnetOptions(nfold = 5)), + qtlDatasetMethods = list(lasso = GlmnetOptions(nfolds = 5)), qtlSumStatsMethods = list(lasso = LassosumOptions()) ) } diff --git a/man/buildMvsusieMixturePrior.Rd b/man/buildMvsusieMixturePrior.Rd new file mode 100644 index 00000000..e40e6e84 --- /dev/null +++ b/man/buildMvsusieMixturePrior.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fineMappingWrappers.R +\name{buildMvsusieMixturePrior} +\alias{buildMvsusieMixturePrior} +\title{Build the mvSuSiE mixture prior} +\usage{ +buildMvsusieMixturePrior( + conditionNames = NULL, + R = length(conditionNames), + mixtureMatrices = NULL, + mixtureWeights = NULL, + weightsTol = 1e-10, + methodArgs = MvsusieMixturePriorOptions() +) +} +\arguments{ +\item{conditionNames}{Character vector naming the conditions the prior is +over, forwarded as \code{include_indices}. On the data-driven arm this +\emph{restricts} the supplied covariances to these conditions by name, +so a mixture spanning more conditions than the fit at hand is subset +here. On the canonical arm it only labels the generated covariances. +NULL leaves the mixture unrestricted and unlabelled.} + +\item{R}{Number of conditions for the canonical mixture. Defaults to +\code{length(conditionNames)}; pass it explicitly when the conditions +are unnamed. Unused when \code{mixtureMatrices} is supplied, since the +mixture then carries its own dimension.} + +\item{mixtureMatrices}{Named list of prior covariance matrices to use as +the mixture. NULL (default) builds the canonical mixture from \code{R} +instead.} + +\item{mixtureWeights}{Numeric vector of mixture weights, named to match +\code{mixtureMatrices}. NULL leaves the engine's uniform weights.} + +\item{weightsTol}{Mixture components weighted at or below this are +dropped.} + +\item{methodArgs}{Options forwarded to +\code{mvsusieR::create_mixture_prior}, built with +\code{\link{MvsusieMixturePriorOptions}}. A bare list is refused: it +cannot be checked against the engine.} +} +\value{ +The prior object \code{mvsusieR::create_mixture_prior()} returns, + ready to pass as \code{prior_variance} to \code{\link{fitMvsusie}} or + \code{\link{fitMvsusieRss}}. +} +\description{ +The single interface to \code{mvsusieR::create_mixture_prior()}, which +every mvSuSiE prior in the package goes through. That function takes +\strong{exactly one} of \code{mixture_prior} and \code{R} and tests the +two with \code{missing()}, so the unused one cannot simply be passed as +\code{NULL}; this wrapper owns that branch. Supply +\code{mixtureMatrices} for a data-driven mixture --- e.g. the reweighted +mr.mash covariances --- or leave it NULL for the canonical uniform +mixture over \code{R} conditions. +} +\examples{ +\donttest{ +# Canonical uniform mixture over two named conditions. +buildMvsusieMixturePrior(conditionNames = c("brain", "blood")) +} +} diff --git a/man/fitFsusie.Rd b/man/fitFsusie.Rd index 7f6f0d75..b8cb2806 100644 --- a/man/fitFsusie.Rd +++ b/man/fitFsusie.Rd @@ -4,7 +4,15 @@ \alias{fitFsusie} \title{Fit fSuSiE on individual-level (X, Y, pos) data} \usage{ -fitFsusie(X, Y, pos, methodArgs = FsusieOptions()) +fitFsusie( + X, + Y, + pos, + maxNumSingleEffects = NULL, + covLev = NULL, + minPurity = NULL, + methodArgs = FsusieOptions() +) } \arguments{ \item{X}{Numeric matrix of genotypes (samples x variants).} @@ -13,6 +21,21 @@ fitFsusie(X, Y, pos, methodArgs = FsusieOptions()) \item{pos}{Numeric vector of trait positions (length \code{ncol(Y)}).} +\item{maxNumSingleEffects}{Maximum number of single effects +(\code{susiF}'s \code{L}). \code{NULL} leaves susiF's own default. An +explicit formal rather than an \code{\link{FsusieOptions}} entry because +the pipeline supplies it from \code{maxNumSingleEffects()} on +\code{\link{CredibleSetParam}}, the same way \code{\link{fitMvsusieRss}} +takes \code{coverage}.} + +\item{covLev}{Credible-set coverage level (\code{susiF}'s \code{cov_lev}). +\code{NULL} leaves susiF's own default; the pipeline supplies +\code{csCoverage()} on \code{\link{CredibleSetParam}}.} + +\item{minPurity}{Minimum credible-set purity (\code{susiF}'s +\code{min_purity}). \code{NULL} leaves susiF's own default; the pipeline +supplies \code{minAbsCorr()} on \code{\link{CredibleSetParam}}.} + \item{methodArgs}{Options forwarded to \code{fsusieR::susiF}, built with \code{\link{FsusieOptions}}. A bare list is refused: it cannot be checked against the engine, so a misspelled option would be silently ignored.} diff --git a/man/fitGlmnetCv.Rd b/man/fitGlmnetCv.Rd new file mode 100644 index 00000000..8929caec --- /dev/null +++ b/man/fitGlmnetCv.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/regularizedRegressionWrappers.R +\name{fitGlmnetCv} +\alias{fitGlmnetCv} +\title{Fit a cross-validated glmnet model} +\usage{ +fitGlmnetCv( + x, + y, + alpha, + family = NULL, + standardize = NULL, + intercept = NULL, + nfolds = NULL, + methodArgs = GlmnetOptions() +) +} +\arguments{ +\item{x}{Numeric design matrix (samples x features).} + +\item{y}{Response: a vector, or a matrix for \code{family = "mgaussian"}.} + +\item{alpha}{Elastic-net mixing parameter in \code{[0, 1]}: \code{1} = +lasso, \code{0} = ridge, \code{0.5} = elastic net. The penalty family is +the calling wrapper's to choose, which is why \code{GlmnetOptions()} +refuses it.} + +\item{family}{glmnet's response family, e.g. \code{"gaussian"} or +\code{"mgaussian"} for a matrix response. \code{NULL} leaves glmnet's +own default.} + +\item{standardize}{Logical; standardize the columns of \code{x} before +fitting. \code{NULL} leaves glmnet's own default (\code{TRUE}).} + +\item{intercept}{Logical; fit an intercept. \code{NULL} leaves glmnet's +own default (\code{TRUE}).} + +\item{nfolds}{Number of cross-validation folds. \code{NULL} leaves +\code{cv.glmnet}'s own default (10).} + +\item{methodArgs}{Further arguments forwarded to +\code{glmnet::cv.glmnet}, built with \code{\link{GlmnetOptions}}. A bare +list is refused, since it cannot be checked against the engine.} +} +\value{ +The fit object returned by \code{glmnet::cv.glmnet}. +} +\description{ +The single interface to \code{glmnet::cv.glmnet}, which every glmnet fit in +the package goes through. It returns the \emph{fit}, not weights: +\code{\link{glmnetWeights}} extracts the \code{lambda.min} coefficient +vector from it, while the \code{mr.mash} initialisation paths need the +coefficient matrix, the intercept row, the lambda sequence and +\code{predict()} on held-out rows --- the same fit-versus-extract split as +\code{\link{fitMvsusie}} and \code{\link{mvsusieWeights}}. +} +\details{ +The four knobs below are formals rather than \code{\link{GlmnetOptions}} +entries because they are what the package's glmnet call sites genuinely +differ on, and leaving them implicit is how three different conventions +came to coexist. Each defaults to \code{NULL}, meaning the argument is not +passed at all and \code{cv.glmnet}'s own default stands; a value supplied +through \code{methodArgs} still wins over the caller's, so none of them is +refused by the constructor. +} +\examples{ +data(eqtlRegionExample) +X <- eqtlRegionExample$X[, 1:30] +y <- eqtlRegionExample$yRes +fit <- fitGlmnetCv(X, y, alpha = 0.5, nfolds = 5) +coef(fit, s = "lambda.min")[1:5] +} diff --git a/man/fitMvsusie.Rd b/man/fitMvsusie.Rd index 3f35d87a..621a8fd8 100644 --- a/man/fitMvsusie.Rd +++ b/man/fitMvsusie.Rd @@ -18,9 +18,10 @@ fitMvsusie( \item{Y}{Numeric matrix of multi-trait / multi-context outcomes (samples x conditions).} -\item{prior_variance}{Prior variance matrix; pass the output of -\code{mvsusieR::create_mixture_prior(R = ncol(Y))} unless you have a -domain-specific prior.} +\item{prior_variance}{Prior variance; build it with +\code{\link{buildMvsusieMixturePrior}}, which is the single interface +to \code{mvsusieR::create_mixture_prior()} and defaults to the +canonical mixture over the conditions it is given.} \item{coverage}{Credible set coverage (default 0.95).} @@ -45,6 +46,6 @@ data(multiTraitData) X <- multiTraitData$X[, 1:60] Y <- multiTraitData$Y fitMvsusie(X, Y, - prior_variance = mvsusieR::create_mixture_prior(R = ncol(Y))) + prior_variance = buildMvsusieMixturePrior(conditionNames = colnames(Y))) } } diff --git a/man/fitSusieInfThenSusieRss.Rd b/man/fitSusieInfThenSusieRss.Rd deleted file mode 100644 index 9162e1b8..00000000 --- a/man/fitSusieInfThenSusieRss.Rd +++ /dev/null @@ -1,60 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fineMappingWrappers.R -\name{fitSusieInfThenSusieRss} -\alias{fitSusieInfThenSusieRss} -\title{Two-stage SuSiE-RSS Fine-mapping} -\usage{ -fitSusieInfThenSusieRss( - z, - R, - n, - args = list(), - SusieInfOptions = list(), - SusieOptions = list(), - fittedModels = NULL -) -} -\arguments{ -\item{z}{Numeric vector of z-scores.} - -\item{R}{LD correlation matrix.} - -\item{n}{Sample size (scalar).} - -\item{args}{Defaults forwarded to both fits, built with -\code{\link{SusieOptions}}. A bare list is refused, since it cannot be -checked; \code{list()} (the default) means no options.} - -\item{SusieInfOptions}{SuSiE-inf-specific overrides, built with -\code{\link{SusieInfOptions}}; they take precedence over \code{args}.} - -\item{SusieOptions}{Standard SuSiE-RSS-specific overrides, built with -\code{\link{SusieOptions}}; they take precedence over \code{args}.} - -\item{fittedModels}{Optional list with pre-fitted \code{$susie} and/or -\code{$susieInf} objects to skip re-fitting.} -} -\value{ -A list with \code{susie} and \code{susieInf} fit objects. -} -\description{ -RSS analog of \code{fitSusieInfThenSusie}. Fits SuSiE-inf via \code{susieRss} -first, then initialises standard SuSiE-RSS from the SuSiE-inf result. The -single pair of fits can be used both for fine-mapping post-processing and -TWAS weight extraction. -} -\examples{ -data(eqtlRegionExample) -X <- eqtlRegionExample$X[, 1:30] -y <- eqtlRegionExample$yRes -ss <- lapply(seq_len(ncol(X)), function(j) { - coef(summary(lm(y ~ X[, j])))[2, 1:2] -}) -stat <- list( - bhat = vapply(ss, `[`, numeric(1), 1L), - shat = vapply(ss, `[`, numeric(1), 2L), - z = vapply(ss, function(s) s[1] / s[2], numeric(1)), - n = rep(nrow(X), ncol(X))) -LD <- cor(X) -fitSusieInfThenSusieRss(z = stat$z, R = LD, n = nrow(X)) -} diff --git a/man/fsusieWeights.Rd b/man/fsusieWeights.Rd index f483508f..786b123f 100644 --- a/man/fsusieWeights.Rd +++ b/man/fsusieWeights.Rd @@ -73,9 +73,9 @@ base <- sin(seq(0, 2 * pi, length.out = nPos)) Y <- matrix(rep(base, each = n), n, nPos) + X[, 1] \%o\% (0.5 * cos(seq(0, pi, length.out = nPos))) pos <- seq_len(nPos) -fit <- fsusieWrapper(X, Y, pos = pos, maxNumSingleEffects = 2, - prior = "mixture_normal_per_scale", maxSnpEm = 10, - covLev = 0.95, minPurity = 0.5, maxScale = 6) +fit <- fitFsusie(X, Y, pos = pos, maxNumSingleEffects = 2, + methodArgs = FsusieOptions(prior = "mixture_normal_per_scale", + max_SNP_EM = 10, max_scale = 6)) fsusieWeights(fsusieFit = fit, X = X, Y = Y, variantIds = colnames(X)) } diff --git a/man/fsusieWrapper.Rd b/man/fsusieWrapper.Rd deleted file mode 100644 index fb7c13bd..00000000 --- a/man/fsusieWrapper.Rd +++ /dev/null @@ -1,69 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/fineMappingWrappers.R -\name{fsusieWrapper} -\alias{fsusieWrapper} -\title{Wrapper for fsusie Function with Automatic Post-Processing} -\usage{ -fsusieWrapper( - X, - Y, - pos, - maxNumSingleEffects, - prior, - maxSnpEm, - covLev, - minPurity, - maxScale, - methodArgs = FsusieOptions() -) -} -\arguments{ -\item{X}{Residual genotype matrix.} - -\item{Y}{Response phenotype matrix.} - -\item{pos}{Genomics position of phenotypes, used for specifying the wavelet -model.} - -\item{maxNumSingleEffects}{The maximum number of the credible set.} - -\item{prior}{method to generate the prior.} - -\item{maxSnpEm}{maximum number of SNP used for learning the prior.} - -\item{covLev}{Coverage level for the credible sets.} - -\item{minPurity}{Minimum purity threshold for credible sets to be retained.} - -\item{maxScale}{numeric, define the maximum of wavelet coefficients used in -the analysis (2^maxScale). Set 10 true by default.} - -\item{methodArgs}{Options forwarded to \code{fsusieR::susiF}, built with -\code{\link{FsusieOptions}}. A bare list is refused: it cannot be checked -against the engine, so a misspelled option would be silently ignored.} -} -\value{ -A modified fsusie object with the susie sets list, correlations for - cs, alpha as df like susie, and without the dummy cs that do not meet the - minimum purity requirement. -} -\description{ -This function serves as a wrapper for the fsusie function, - facilitating automatic post-processing such as removing dummy credible sets - (cs) that don't meet the minimum purity threshold and calculating - correlations for the remaining cs. The function parameters are identical to - those of the fSuSiE function. -} -\examples{ -data(eqtlRegionExample) -X <- eqtlRegionExample$X[, 1:50] -n <- nrow(X) -nPos <- 16 -base <- sin(seq(0, 2 * pi, length.out = nPos)) -Y <- matrix(rep(base, each = n), n, nPos) + - X[, 1] \%o\% (0.5 * cos(seq(0, pi, length.out = nPos))) -pos <- seq_len(nPos) -fsusieWrapper(X, Y, pos = pos, maxNumSingleEffects = 2, - prior = "mixture_normal_per_scale", maxSnpEm = 10, - covLev = 0.95, minPurity = 0.5, maxScale = 6) -} diff --git a/man/mvsusieWeights.Rd b/man/mvsusieWeights.Rd index 0ac56e29..eca38ab8 100644 --- a/man/mvsusieWeights.Rd +++ b/man/mvsusieWeights.Rd @@ -25,7 +25,7 @@ data(multiTraitData) X <- multiTraitData$X[, 1:60] Y <- multiTraitData$Y fit <- fitMvsusie(X = X, Y = Y, - prior_variance = mvsusieR::create_mixture_prior(R = ncol(Y)), + prior_variance = buildMvsusieMixturePrior(R = ncol(Y)), methodArgs = MvsusieOptions(max_iter = 50)) mvsusieWeights(mvsusieFit = fit) } diff --git a/man/twasEngineOptions.Rd b/man/twasEngineOptions.Rd index ac60502c..f0e88157 100644 --- a/man/twasEngineOptions.Rd +++ b/man/twasEngineOptions.Rd @@ -77,9 +77,21 @@ Options forwarded to the engine behind a TWAS weight method, itself is refused --- \code{GlmnetOptions(alpha=)} and \code{NcvregOptions(penalty=)} select the penalty family, which is the calling wrapper's to set. See \code{\link{MethodOptions}}. + + \code{GlmnetOptions} rejects two names by hand, which the \code{...} + check above would otherwise let through in silence. \code{nfold} is not + a \code{cv.glmnet} argument at all --- the fold count is \code{nfolds}, + and the shorter spelling only ever worked through R's partial matching. + \code{parallel} is a real one, but pecotmr does not expose it: glmnet + parallelises over its CV folds through a \code{foreach} backend that + this package neither declares nor registers, and weight fitting already + runs inside \code{BiocParallel} workers over the same folds, so + enabling it would nest two frameworks and oversubscribe the machine. + The glasso and univariate-glmnet \code{mr.mash} initialisations are + single-threaded for the same reason. } \examples{ -GlmnetOptions(nfold = 10) +GlmnetOptions(nfolds = 10) QggOptions(updateB = TRUE) LassosumOptions(thr = 1e-4) } diff --git a/tests/testthat/test_MethodOptions.R b/tests/testthat/test_MethodOptions.R index f5a0a63f..f166208f 100644 --- a/tests/testthat/test_MethodOptions.R +++ b/tests/testthat/test_MethodOptions.R @@ -397,6 +397,7 @@ test_that("a constructor is unchecked only when its engine takes dots", { DprOptions = "RcppDPR::fit_model", GlmnetOptions = "glmnet::cv.glmnet", MashCorEmOptions = "mashr::mash_estimate_corr_em", + MvsusieMixturePriorOptions = "mvsusieR::create_mixture_prior", MvsusieRssOptions = "mvsusieR::mvsusie_rss", NcvregOptions = "ncvreg::cv.ncvreg", PoolrOptions = "poolr::fisher", diff --git a/tests/testthat/test_MethodParam.R b/tests/testthat/test_MethodParam.R index 7abe73f3..0dfbaf56 100644 --- a/tests/testthat/test_MethodParam.R +++ b/tests/testthat/test_MethodParam.R @@ -157,11 +157,11 @@ test_that("malformed slots are rejected with the caller's own error", { test_that("TwasWeightsMethodsParam shares the rules with its own registry", { p <- TwasWeightsMethodsParam( methods = list(susie = list()), - qtlDatasetMethods = list(lasso = GlmnetOptions(nfold = 5)), + qtlDatasetMethods = list(lasso = GlmnetOptions(nfolds = 5)), qtlSumStatsMethods = list(lasso = LassosumOptions()) ) expect_s4_class(p, "TwasWeightsMethodsParam") - expect_equal(p$qtlDatasetMethods$lasso$nfold, 5) + expect_equal(p$qtlDatasetMethods$lasso$nfolds, 5) # The registries differ: prsCs is summary-statistics only, enet # individual only, and both are methods a fine-mapping Param never has. expect_error( @@ -469,7 +469,7 @@ test_that("show renders a slot holding per-method entries", { multi <- paste( capture.output(show(TwasWeightsMethodsParam( methods = list(susie = list()), - qtlDatasetMethods = list(lasso = GlmnetOptions(nfold = 5)), + qtlDatasetMethods = list(lasso = GlmnetOptions(nfolds = 5)), qtlSumStatsMethods = list(lasso = LassosumOptions()) ))), collapse = " " diff --git a/tests/testthat/test_fineMappingPipeline.R b/tests/testthat/test_fineMappingPipeline.R index fb2dd94f..48002ace 100644 --- a/tests/testthat/test_fineMappingPipeline.R +++ b/tests/testthat/test_fineMappingPipeline.R @@ -1462,12 +1462,19 @@ test_that("fineMappingPipeline(QtlDataset): fsusie with single trait rejected", function(R, ...) list(R = R) } .fmp_mockSusiF <- function() { - function(X, Y, pos) { + # L / cov_lev / min_purity are refused by FsusieOptions as pipeline-owned, + # so the pipeline is their only source: record them here so a test can + # assert they actually arrive. `...` keeps the stub tolerant of further + # susiF arguments. + function(X, Y, pos, L = NULL, cov_lev = NULL, min_purity = NULL, ...) { list( token = "fsusie", n_X_cols = ncol(X), n_Y_cols = ncol(Y), - pos = pos + pos = pos, + L = L, + cov_lev = cov_lev, + min_purity = min_purity ) } } @@ -1936,6 +1943,45 @@ test_that("fineMappingPipeline(QtlDataset): jointSpec='trait' with fsusie wires expect_true("jointTraits" %in% pecotmr:::.tupleColumnNames(res)) }) +test_that("fineMappingPipeline(QtlDataset): fsusie receives the CredibleSetParam settings", { + # Regression: .fmOwnedEngineArgs() refuses L / cov_lev / min_purity in + # FsusieOptions on the grounds that the pipeline supplies them, but for a + # while nothing did -- the row described the retired fsusieWrapper() while + # the pipeline called fitFsusie(). `--L` was silently ignored and susiF + # ran at its own default. The mocked susiF records what it was handed. + qd <- .fmp_makeQtlDataset( + contexts = "brain", + traits = c("ENSG_A", "ENSG_B") + ) + local_mocked_bindings( + extractBlockGenotypes = .fmp_mockExtractor(), + .fmPostprocessOne = .fmp_mockPostprocess(), + .package = "pecotmr" + ) + local_mocked_bindings( + susiF = .fmp_mockSusiF(), + .package = "fsusieR" + ) + res <- suppressMessages( + fineMappingPipeline( + qd, + methods = "fsusie", + cisWindow = 1000L, + jointSpecification = "trait", + credibleSetParam = CredibleSetParam( + maxNumSingleEffects = 3L, + coverage = 0.9, + minAbsCorr = 0.4 + ) + ) + ) + # .fmp_mockPostprocess stores the raw fit under `payload`. + recorded <- susieFit(res[1, ])$payload + expect_equal(recorded$L, 3L) + expect_equal(recorded$cov_lev, 0.9) + expect_equal(recorded$min_purity, 0.4) +}) + test_that("fineMappingPipeline(QtlDataset): fsusie multi-trait per context dispatch", { qd <- .fmp_makeQtlDataset( contexts = "brain", @@ -5239,14 +5285,21 @@ test_that("CredibleSetParam validity rejects non-scalar counts and flags", { test_that("SerOptions and FsusieOptions refuse what their wrappers own", { # The susie_ser and susiF arms of .fmOwnedEngineArgs(): ser takes only - # the coverage half of csHome, and susiF owns five wrapper settings on - # top of its data arguments. + # the coverage half of csHome, and susiF owns the three settings the + # pipeline actually supplies to it. expect_error(SerOptions(z = 1), "supplied by pecotmr") expect_error(SerOptions(n = 100L), "supplied by pecotmr") expect_error(SerOptions(coverage = 0.9), "csCoverage") expect_error(FsusieOptions(L = 2), "maxNumSingleEffects") - expect_error(FsusieOptions(min_purity = 0.5), "minPurity") - expect_error(FsusieOptions(max_scale = 6), "maxScale") + expect_error(FsusieOptions(min_purity = 0.5), "minAbsCorr") + expect_error(FsusieOptions(cov_lev = 0.9), "csCoverage") + # ... and ONLY those. max_scale / max_SNP_EM / prior were refused while + # this row still described the retired fsusieWrapper(); nothing supplies + # them on the pipeline path, so refusing them closed the only way to set + # them at all. + expect_s4_class(FsusieOptions(max_scale = 6), "MethodOptions") + expect_s4_class(FsusieOptions(max_SNP_EM = 20), "MethodOptions") + expect_s4_class(FsusieOptions(prior = "mixture_normal"), "MethodOptions") # ser does NOT own L -- susie_ser is single-effect by construction, so # the csHome entry it takes is coverage alone. expect_s4_class(SerOptions(prior_variance = 1), "MethodOptions") diff --git a/tests/testthat/test_fineMappingRow.R b/tests/testthat/test_fineMappingRow.R index dc5b1a03..6e9c3f36 100644 --- a/tests/testthat/test_fineMappingRow.R +++ b/tests/testthat/test_fineMappingRow.R @@ -1007,7 +1007,8 @@ test_that("credibleSetSummary aggregates across a collection with entry identity rep(0.1, length(vids)) } # Canonical CS labels + purity from the fit's CS membership, mirroring what - # fsusieWrapper writes into a real FMR's topLoci (cs_95 / cs_95_purity). + # The fsusie post-processing writes into a real FMR's topLoci + # (cs_95 / cs_95_purity). cs95 <- rep("fsusie_0", length(vids)) for (l in seq_along(fit$cs)) { idx <- fit$cs[[l]] diff --git a/tests/testthat/test_fineMappingWrappers.R b/tests/testthat/test_fineMappingWrappers.R index e85eb091..cbfe9044 100644 --- a/tests/testthat/test_fineMappingWrappers.R +++ b/tests/testthat/test_fineMappingWrappers.R @@ -416,7 +416,11 @@ test_that("postprocessFinemappingFits keeps all effects when V is NULL", { skip_if_not_installed("susieR") p <- 5 L <- 3 - fake_output <- makeFakeSusieOutput(p, maxNumSingleEffects = L, has_V = FALSE) + fake_output <- makeFakeSusieOutput( + p, + maxNumSingleEffects = L, + has_V = FALSE + ) R <- diag(p) colnames(R) <- rownames(R) <- names(fake_output$pip) post <- postprocessFinemappingFits( @@ -1722,8 +1726,6 @@ test_that("posterior_mean equals colSums(alpha*mu); posterior_sd equals sqrt(pma }) -context("fsusieWrapper") - # ---- cal_purity ---- test_that("cal_purity with min method and single element CS", { set.seed(42) @@ -1784,111 +1786,6 @@ test_that("cal_purity with multiple credible sets", { }) # ---- fsusieGetCs ---- -# ---- fsusieWrapper ---- -test_that("fsusieWrapper errors when fsusieR is not installed", { - skip_if( - requireNamespace("fsusieR", quietly = TRUE), - "fsusieR is installed, skipping not-installed test" - ) - set.seed(1) - X <- matrix(rnorm(50), nrow = 10, ncol = 5) - Y <- matrix(rnorm(40), nrow = 10, ncol = 4) - expect_error( - fsusieWrapper( - X = X, - Y = Y, - pos = seq_len(4), - maxNumSingleEffects = 3, - prior = "mixture_normal", - maxSnpEm = 100, - covLev = 0.95, - minPurity = 0.5, - maxScale = 5 - ), - "fsusieR" - ) -}) - -test_that("fsusieWrapper low-purity branch sets cs to list(NULL) and cs_corr to NULL", { - skip_if_not_installed("fsusieR") - fake_fit <- list( - cs = list(c(1, 2), c(3)), - purity = c(0.1, 0.05), # all < min_purity = 0.5 - pip = c(0.1, 0.2, 0.3, 0.05, 0.05), - alpha = list( - matrix(0.1, nrow = 2, ncol = 5), - matrix(0.1, nrow = 2, ncol = 5) - ) - ) - local_mocked_bindings( - susiF = function(...) fake_fit, - .package = "fsusieR" - ) - set.seed(1) - X <- matrix(rnorm(50), nrow = 10, ncol = 5) - Y <- matrix(rnorm(40), nrow = 10, ncol = 4) - out <- fsusieWrapper( - X = X, - Y = Y, - pos = seq_len(4), - maxNumSingleEffects = 3, - prior = "mixture_normal", - maxSnpEm = 100, - covLev = 0.95, - minPurity = 0.5, - maxScale = 5 - ) - expect_equal(out$cs, list(NULL)) - expect_equal(out$sets$cs, list(NULL)) - expect_null(out$cs_corr) -}) - -test_that("fsusieWrapper high-purity branch builds sets", { - skip_if_not_installed("fsusieR") - set.seed(2) - p <- 5 - fake_fit <- list( - cs = list(c(1, 2), c(3, 4)), - purity = c(0.95, 0.9), # all > min_purity = 0.5 - pip = c(0.4, 0.4, 0.6, 0.6, 0.1), - alpha = list( - matrix( - rep(c(0.4, 0.4, 0.05, 0.05, 0.1), each = 2), - nrow = 2, - byrow = FALSE - ), - matrix( - rep(c(0.05, 0.05, 0.45, 0.4, 0.05), each = 2), - nrow = 2, - byrow = FALSE - ) - ) - ) - local_mocked_bindings( - susiF = function(...) fake_fit, - cal_cor_cs = function(obj, X) matrix(c(1, 0.9, 0.9, 1), nrow = 2), - .package = "fsusieR" - ) - X <- matrix(rnorm(10 * p), nrow = 10, ncol = p) - Y <- matrix(rnorm(40), nrow = 10, ncol = 4) - out <- fsusieWrapper( - X = X, - Y = Y, - pos = seq_len(4), - maxNumSingleEffects = 3, - prior = "mixture_normal", - maxSnpEm = 100, - covLev = 0.95, - minPurity = 0.5, - maxScale = 5 - ) - expect_length(out$sets$cs, 2) - expect_equal(names(out$sets$cs), c("L1", "L2")) - # cs_corr is no longer stored on the fit (computed on demand later by - # computeCsCorrelation()). - expect_null(out$cs_corr) - expect_equal(out$sets$requested_coverage, 0.95) -}) test_that("fsusieGetCs creates susie-like sets", { set.seed(42) @@ -4571,49 +4468,6 @@ test_that("mvsusie / fsusie wrappers: numeric guards fire", { ) }) -test_that("fitSusieInfThenSusieRss: argument guards fire", { - expect_error( - fitSusieInfThenSusieRss(z = "nope", R = NULL, n = 100), - "z.*Must be of type 'numeric'" - ) - expect_error( - fitSusieInfThenSusieRss(z = 1, R = NULL, n = -1), - "n.*is not >= 0" - ) - expect_error( - fitSusieInfThenSusieRss(z = 1, R = NULL, n = 1, args = 1L), - "must be built with SusieOptions" - ) - # A bare list is the case the constructors exist to refuse: it cannot be - # checked, so a misspelling would reach susieR and be ignored. - expect_error( - fitSusieInfThenSusieRss( - z = 1, - R = NULL, - n = 1, - SusieInfOptions = list(L = 5) - ), - "must be built with SusieInfOptions" - ) - expect_error( - fitSusieInfThenSusieRss( - z = 1, - R = NULL, - n = 1, - SusieOptions = SusieOptions(nope = 1) - ), - "unknown argument" - ) - expect_error( - fitSusieInfThenSusieRss( - z = 1, - R = NULL, - n = 1, - fittedModels = 1L - ), - "fittedModels.*Must be of type 'list'" - ) -}) test_that(".fmwConcat answers an empty list for no pieces", { expect_identical(pecotmr:::.fmwConcat(list()), list()) @@ -4727,6 +4581,117 @@ test_that(".fmSerScreen: disables on 0, skips no-signal, keeps signal + adaptive expect_true(fn(X, yNull, NA)) # malformed cutoff -> advisory keep }) +# Capture what buildMvsusieMixturePrior() hands the engine, which is the only +# thing the wrapper is responsible for: the engine itself is mvsusieR's. +.fmw_captureMixturePrior <- function(expr) { + captured <- NULL + local_mocked_bindings( + create_mixture_prior = function(...) { + captured <<- list(...) + "PRIOR" + }, + .package = "mvsusieR" + ) + force(expr) + captured +} + +test_that("buildMvsusieMixturePrior: canonical arm names only R", { + args <- .fmw_captureMixturePrior( + buildMvsusieMixturePrior(conditionNames = c("c1", "c2")) + ) + # create_mixture_prior() tests the two with missing(), so the unused name + # must be ABSENT, not NULL -- supplying both trips its own stop(). + expect_false("mixture_prior" %in% names(args)) + expect_equal(args$R, 2L) + expect_equal(args$include_indices, c("c1", "c2")) + expect_equal(args$weights_tol, 1e-10) +}) + +test_that("buildMvsusieMixturePrior: data-driven arm omits R", { + U <- list(compA = diag(2), compB = diag(2) * 2) + args <- .fmw_captureMixturePrior(buildMvsusieMixturePrior( + conditionNames = c("c1", "c2"), + mixtureMatrices = U, + mixtureWeights = c(compA = 0.25, compB = 0.75), + weightsTol = 1e-8 + )) + expect_false("R" %in% names(args)) + expect_equal(names(args$mixture_prior$matrices), c("compA", "compB")) + expect_equal(args$mixture_prior$weights, c(compA = 0.25, compB = 0.75)) + expect_equal(args$include_indices, c("c1", "c2")) + expect_equal(args$weights_tol, 1e-8) +}) + +test_that("buildMvsusieMixturePrior: NULL mixtureWeights stay uniform", { + args <- .fmw_captureMixturePrior( + buildMvsusieMixturePrior(mixtureMatrices = list(compA = diag(2))) + ) + # Absent rather than NULL: create_mixture_prior() fills uniform weights + # when `weights` is missing from the mixture. + expect_false("weights" %in% names(args$mixture_prior)) + expect_equal(names(args$mixture_prior$matrices), "compA") +}) + +test_that("buildMvsusieMixturePrior: unnamed conditions need an explicit R", { + # length(NULL) is 0, which create_cov_canonical() turns into a 0 x 0 + # subscript error deep in the engine. Caught here instead. + expect_error(buildMvsusieMixturePrior(), "Must be >= 1") + args <- .fmw_captureMixturePrior(buildMvsusieMixturePrior(R = 3)) + expect_equal(args$R, 3) + expect_false("include_indices" %in% names(args)) +}) + +test_that("buildMvsusieMixturePrior: conditionNames must agree with R", { + expect_error( + buildMvsusieMixturePrior(conditionNames = c("c1", "c2"), R = 3), + "names 2 conditions but `R` is 3" + ) +}) + +test_that("buildMvsusieMixturePrior: methodArgs forwarded, list refused", { + args <- .fmw_captureMixturePrior(buildMvsusieMixturePrior( + conditionNames = c("c1", "c2"), + methodArgs = MvsusieMixturePriorOptions(singletons = FALSE) + )) + expect_false(args$singletons) + expect_error( + buildMvsusieMixturePrior(R = 2, methodArgs = list(singletons = FALSE)), + "MvsusieMixturePriorOptions" + ) +}) + +test_that("MvsusieMixturePriorOptions refuses the arguments the builder owns", { + for (owned in c("mixture_prior", "R", "include_indices", "weights_tol")) { + expect_error( + exec(MvsusieMixturePriorOptions, !!!set_names(list(1), owned)), + "supplied by pecotmr" + ) + } +}) + +test_that(".fmFoldWeightsMv: NULL mvPrior uses buildMvsusieMixturePrior", { + # Regression: the fold path used to inline its own + # create_mixture_prior(R = ncol(Ytr)) -- the last fragment of the legacy + # hardcoded canonical prior -- so its prior came back unlabelled while + # the full-data path's did not. + Xtr <- matrix(rnorm(20), 10, 2, dimnames = list(NULL, c("v1", "v2"))) + Ytr <- matrix(rnorm(20), 10, 2, dimnames = list(NULL, c("c1", "c2"))) + args <- .fmw_captureMixturePrior({ + local_mocked_bindings( + fitMvsusie = function(...) "FIT", + mvsusieWeights = function(mvsusieFit) { + matrix(0, 2, 2, dimnames = list(c("v1", "v2"), NULL)) + }, + .fmLeanFoldFit = function(fit, token) NULL, + .package = "pecotmr" + ) + pecotmr:::.fmFoldWeightsMv(Xtr, Ytr, 0.95, NULL, NULL) + }) + expect_equal(args$R, 2L) + expect_equal(args$include_indices, c("c1", "c2")) +}) + test_that(".buildMvsusieReweightedPrior: canonical fallback when no usable fit", { bp <- function(...) pecotmr:::.buildMvsusieReweightedPrior(...) # No fit at all -> canonical prior, residualVariance NULL. diff --git a/tests/testthat/test_regularizedRegressionWrappers.R b/tests/testthat/test_regularizedRegressionWrappers.R index c9f60288..027e3257 100644 --- a/tests/testthat/test_regularizedRegressionWrappers.R +++ b/tests/testthat/test_regularizedRegressionWrappers.R @@ -17,7 +17,7 @@ context("SS-TWAS: weights, pipeline, and omnibus combination") # # What remains: tests of the internal SuSiE-RSS weight extractors that # are still present in `R/twasWeights.R` (`.susieRssExtractWeights`, -# `susieRssWeights`, `susieInfRssWeights`, `fitSusieInfThenSusieRss`). +# `susieRssWeights`, `susieInfRssWeights`). # ============================================================================= # SuSiE-RSS weight extraction @@ -64,21 +64,6 @@ test_that("mrmashWeights fitRetention: slim omits the full fit, full keeps it", # Two-stage SuSiE-RSS fitting # ============================================================================= -test_that("fitSusieInfThenSusieRss returns two fits", { - skip_if_not_installed("susieR") - set.seed(42) - p <- 20 - n <- 500 - R <- diag(p) - z <- rnorm(p) - fits <- fitSusieInfThenSusieRss(z, R, n, args = SusieOptions(max_iter = 5)) - expect_true(is.list(fits)) - expect_true("susie" %in% names(fits)) - expect_true("susieInf" %in% names(fits)) - expect_true("susieInf" %in% class(fits$susieInf)) - expect_true("susieRss" %in% class(fits$susie)) -}) - # === Tests migrated from test_mrmashWrapper.R (mr.mash + glasso/glmnet coef helpers) === test_that("computeW0 returns uniform weights when ncomps == 1", { @@ -435,6 +420,136 @@ test_that("computeCoefficientsUnivGlmnet handles NA in Y", { # -------------------------------- individual -------------------------------- +# Capture what reaches cv.glmnet, which is all fitGlmnetCv() is responsible +# for -- the engine itself is glmnet's. The fake fit is not a usable +# cv.glmnet object, so an extractor downstream of the call errors; `try` +# swallows that, since only the captured arguments are under test. +.rrw_captureGlmnet <- function(expr) { + captured <- NULL + local_mocked_bindings( + cv.glmnet = function(...) { + captured <<- list(...) + structure(list(), class = "cv.glmnet") + }, + .package = "glmnet" + ) + try(force(expr), silent = TRUE) + captured +} + +test_that("fitGlmnetCv: NULL formals leave glmnet's own defaults", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + args <- .rrw_captureGlmnet(fitGlmnetCv(f$X, f$y, alpha = 0.5)) + expect_equal(args$alpha, 0.5) + # An explicit NULL would reach cv.glmnet and override its default with + # nothing, so an unset formal must be ABSENT from the call. + for (nm in c("family", "standardize", "intercept", "nfolds")) { + expect_false(is_in(nm, names(args))) + } +}) + +test_that("fitGlmnetCv: set formals use glmnet's own argument names", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + args <- .rrw_captureGlmnet(fitGlmnetCv( + f$X, + f$y, + alpha = 1, + family = "mgaussian", + standardize = TRUE, + intercept = FALSE, + nfolds = 3 + )) + expect_equal(args$family, "mgaussian") + expect_true(args$standardize) + expect_false(args$intercept) + expect_equal(args$nfolds, 3) + expect_false(is_in("nfold", names(args))) +}) + +test_that("fitGlmnetCv: methodArgs win over the caller's defaults", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + args <- .rrw_captureGlmnet(fitGlmnetCv( + f$X, + f$y, + alpha = 0.5, + nfolds = 5, + methodArgs = GlmnetOptions(nfolds = 3) + )) + expect_equal(args$nfolds, 3) + expect_error( + fitGlmnetCv(f$X, f$y, alpha = 0.5, methodArgs = list(nfolds = 3)), + "GlmnetOptions" + ) +}) + +test_that("GlmnetOptions: nfold is refused and names nfolds instead", { + expect_error(GlmnetOptions(nfold = 5), "fold count is `nfolds`") +}) + +test_that("GlmnetOptions: parallel is refused outright", { + # glmnet's `parallel` needs a foreach backend pecotmr does not declare, + # and weight fitting already runs inside BiocParallel workers over the + # same CV folds -- so there is no value a caller could usefully set. + # cv.glmnet takes `...`, so nothing else would report it. + expect_error(GlmnetOptions(parallel = TRUE), "cannot be set here") + expect_error(GlmnetOptions(parallel = FALSE), "foreach backend") + # A refusal must not swallow the names that ARE settable. + expect_equal(GlmnetOptions(nfolds = 3)$nfolds, 3) +}) + +test_that("fitGlmnetCv: the engine is never asked to parallelise", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + args <- .rrw_captureGlmnet(fitGlmnetCv(f$X, f$y, alpha = 0.5)) + expect_false(is_in("parallel", names(args))) +}) + +test_that("glmnetWeights: the fold count reaches the engine as nfolds", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + args <- .rrw_captureGlmnet(glmnetWeights(f$X, f$y, alpha = 0.5)) + expect_equal(args$nfolds, 5) + expect_false(args$standardize) + expect_true(args$intercept) + # Regression: pecotmr passed `nfold = 5`, which reached cv.glmnet only by + # partial matching, so a user who spelled the engine's own `nfolds` + # supplied BOTH and glmnet failed with a message about `family`. + w <- glmnetWeights(f$X, f$y, 0.5, methodArgs = GlmnetOptions(nfolds = 3)) + expect_length(as.numeric(w), f$p) +}) + +test_that("mr.mash glmnet initialisation routes through fitGlmnetCv", { + skip_if_not_installed("glmnet") + f <- .rrwXy() + Y <- cbind(f$y, f$y + rnorm(length(f$y))) + gl <- .rrw_captureGlmnet(computeCoefficientsGlasso( + f$X, + Y, + standardize = TRUE, + numThreads = 1L + )) + expect_equal(gl$family, "mgaussian") + expect_equal(gl$alpha, 1) + expect_true(gl$standardize) + # `parallel = FALSE` was a restatement of cv.glmnet's own default. + expect_false(is_in("parallel", names(gl))) + uv <- .rrw_captureGlmnet(pecotmr:::.linreg( + 1L, + f$X, + Y, + alpha = 0.5, + standardize = FALSE, + Xnew = NULL + )) + expect_equal(uv$family, "gaussian") + expect_equal(uv$alpha, 0.5) + expect_false(uv$standardize) + expect_false(is_in("parallel", names(uv))) +}) + test_that("lassoWeights / enetWeights (glmnet) return length-p weights", { skip_if_not_installed("glmnet") f <- .rrwXy() diff --git a/tests/testthat/test_twasWeightsPipeline.R b/tests/testthat/test_twasWeightsPipeline.R index cb10da74..ba6fe4dd 100644 --- a/tests/testthat/test_twasWeightsPipeline.R +++ b/tests/testthat/test_twasWeightsPipeline.R @@ -1914,7 +1914,12 @@ test_that(".twasCheckMultivariateY: mvsusie/mrmash error with single trait, sing sampleIds = character(), pgenPtr = NULL ) - GwasSumStats(studyName = "g1", entry = list(gr), genome = "hg19", ldSketch = gh) + GwasSumStats( + studyName = "g1", + entry = list(gr), + genome = "hg19", + ldSketch = gh + ) } test_that(".twasAssertQcd: errors when qcInfo is empty", { @@ -3640,6 +3645,29 @@ test_that(".solveEnsembleNnls: solver failure and all-zero fall back to equal we expect_equal(z2, c(0.5, 0.5)) }) +test_that(".solveEnsembleGlmnet: routes through fitGlmnetCv", { + skip_if_not_installed("glmnet") + set.seed(3) + P <- matrix(rnorm(60), 20, 3) + yObs <- rnorm(20) + captured <- NULL + local_mocked_bindings( + cv.glmnet = function(...) { + captured <<- list(...) + structure(list(), class = "cv.glmnet") + }, + .package = "glmnet" + ) + # coef() on the fake fit errors downstream of the call; only the + # arguments that reached the engine are under test here. + try(pecotmr:::.solveEnsembleGlmnet(P, yObs, 3L), silent = TRUE) + # `lower.limits = 0` is what makes the stacking regression non-negative; + # it has no pecotmr formal, so it travels in the engine's options bag. + expect_equal(captured$lower.limits, 0) + expect_false(captured$intercept) + expect_equal(captured$alpha, 1) +}) + test_that(".solveEnsembleGlmnet: solver failure falls back to equal weights", { skip_if_not_installed("glmnet") P <- matrix(rnorm(20), 10, 2) @@ -4405,7 +4433,12 @@ test_that("a multivariate group resumes from the cache only in full", { # there was no per-group question being asked. There is -- it is just # all-or-nothing, because reusing some contexts' weights while refitting # the others would mix two fits inside a group that was fitted together. - w <- .tp_tw(studyName = "s1", context = "c1", trait = "t1", method = "mrmash") + w <- .tp_tw( + studyName = "s1", + context = "c1", + trait = "t1", + method = "mrmash" + ) # One context cached out of two -> no reuse. expect_null(pecotmr:::.twasMvCacheHits(