Multivariate adaptive shrinkage across contexts with mash
pecotmr authors
2026-10-22
Source:vignettes/mash-pipeline.Rmd
mash-pipeline.RmdOverview
A QTL measured in several contexts — tissues, cell types, conditions — is usually not a set of independent experiments. Effects are shared, and an effect that looks null in one context is often a real effect that context was underpowered to see.
mash learns the
covariance structure of effect sharing from the data and shrinks each
variant’s per-context estimate toward it. mashPipeline()
drives the whole fit; the chainable steps behind it let you swap out any
stage.
The input is a named list of summary-statistic collections split into partitions, each with a distinct job:
| Partition | Role |
|---|---|
strong |
discovery variants; the data-driven covariances are learned here |
random |
unbiased background; the mixture weights are fit here |
null |
null variants; the residual correlation is estimated here |
strong and random are required.
null is optional and used only by the residual-correlation
estimators that need it.
Fitting chunks below are run with results = "hide":
mashr and flashier narrate their progress on
stdout, which knitr cannot filter, so each fit is split into a compute
chunk and a display chunk. The code is otherwise exactly what you would
type.
## QtlSumStats: 3 entries, genome build hg19
## 1 studies, 3 contexts, 1 traits
## LD reference: plink1 @ pecotmr://extdata/toy_canonicalThe one-shot fit
sumStatsList <- list(strong = ss, random = ss, null = ss)
fit <- mashPipeline(
sumStatsList,
alpha = 0,
prior = MashPriorParam(nPcs = 2L)
)
length(fit$U)## [1] 17## null equal_effects muscle ED_PCA_1 ED_PCA_2
## 0.6708 0.3145 0.0147 0.0000 0.0000fit$U is the prior covariance list and
fit$w the estimated mixture weights over it. A large weight
on null is normal and healthy: most variants are not
QTLs.
alpha must match the scale you are on:
alpha = 0 for effect sizes (BETA/SE),
alpha = 1 for Z-scores. inputScale controls
which pair is extracted from the sumstats and defaults to
"auto" — BETA + SE when every entry carries both, otherwise
Z with unit standard errors.
Residual correlation
Conditions measured on overlapping samples have correlated residuals,
and mash needs an estimate of that correlation (Vhat) to
avoid reading shared noise as shared signal.
residualCorrelationMethod chooses the estimator:
| Value | Estimator | Needs |
|---|---|---|
"identity" (default) |
no residual correlation | — |
"simple" |
mashr::estimate_null_correlation_simple() |
null |
"simpleSpecific" |
nearest PD correlation of cov(nullZ)
|
null |
"corshrink" |
adaptive shrinkage of the null correlation |
null, CorShrink |
"mle" |
EM refinement against the prior |
random, priorCovariances
|
The default is "identity" because it always works,
not because it is the best choice. Multi-context QTL
data is the overlapping-sample case par excellence — the same donors
appear in every context — so the identity assumption is usually wrong.
Prefer a data-driven estimate whenever you have a null set:
fitSimple <- mashPipeline(
sumStatsList,
alpha = 0,
residualCorrelationMethod = "simple",
prior = MashPriorParam(priorCovariances = fit$U)
)## null equal_effects blood muscle ED_PCA_1
## 0.7091 0.2350 0.0315 0.0244 0.0000Passing fit$U back in holds the prior fixed, so
Vhat is the only thing that differs from the fit above and
the shift in the weights is attributable to it alone. Re-estimating the
prior here would confound the comparison – and cost a second data-driven
estimate for a prior you already have.
Supplying a null partition while leaving the default in
place reports that the partition is going unused, since that combination
is more often an oversight than an intent.
You can also skip estimation entirely and pass a Vhat
computed elsewhere — useful when it was estimated on a larger reference
and shipped as a static artefact:
mashPipeline(sumStatsList, alpha = 0, residualCorrelation = myVhat)The chain, step by step
mashPipeline() is a wrapper over four steps. Take them
apart when you want to reuse an intermediate, try several priors against
one Vhat, or compute posteriors for a set other than the
one you fit on.
vhat <- mashResidualCorrelation(sumStatsList, alpha = 0, method = "simple")
prior <- mashPriorCovariances(sumStatsList, alpha = 0, vhat = vhat, nPcs = 2L)
model <- mashModelFit(
sumStatsList,
alpha = 0,
priorCovariances = prior,
vhat = vhat
)
class(model)## [1] "mash"mashPriorCovariances() returns
list(U, w, loglik), while
mashCovarianceComponents() returns a bare covariance list.
Either can be handed to any consumer — the wrapper is unwrapped for you,
so there is no need to reach for $U.
mashModelFit() returns a mashr model
object. That is not the same thing as the MashPrior class,
which packages a prior for use as a TWAS weight method; see TWAS weights.
Posteriors
mashPosterior() applies the fitted model to a target
set. Fitting on random and computing posteriors on
strong is the standard split — the weights should not be
learned from the same variants they are applied to.
post <- mashPosterior(model, ss, alpha = 0, vhat = vhat)
names(post)## [1] "PosteriorMean" "PosteriorSD" "lfdr" "NegativeProb"
## [5] "lfsr" "PosteriorCov"
dim(post$PosteriorMean)## [1] 200 3| Matrix | Meaning |
|---|---|
PosteriorMean / PosteriorSD
|
shrunk effect and its uncertainty |
lfsr |
local false sign rate — the quantity to threshold on |
lfdr |
local false discovery rate |
NegativeProb |
posterior probability the effect is negative |
PosteriorCov |
per-variant across-context covariance |
lfsr is the usual significance measure here: it asks
whether the sign of an effect is confidently known, which is a
more honest question than whether it is exactly zero.
## [1] 24Choosing the covariance components
The data-driven part of the prior is built from
components and then refined by an engine:
mashPriorCovariances(
sumStatsList,
alpha = 0,
vhat = vhat,
components = c("canonical", "pca", "flash", "flashNonneg"),
engine = "covEd", # or "covUdr", or "none"
nPcs = 2L
)"covEd" is the default and is mashr’s own extreme
deconvolution. "covUdr" is opt-in and has known numerical
issues; its TED update additionally assumes i.i.d. (Z-scale) data, so it
is requested as
CovUdrOptions(fit = UdFitOptions(unconstrained.update = "ted"))
rather than as a separate engine name. udr is two functions –
ud_init() then ud_fit() – so each one’s
arguments travel in their own bundle and are checked against their own
source. Its dominant cost is
CovUdrOptions(init = UdInitOptions(n_unconstrained = ))
(default 50), which is worth cutting sharply when you have few
conditions.
The engine refines only the data-driven components.
Canonical components are fixed structural hypotheses – identity,
singletons, equal effects – and reach the prior unmodified, matching
mashr’s own workflow. engine = "none" uses the data-driven
components as generated, without refinement.
To bypass the covariance chain entirely and fit weights against a prior you already have:
mashPipeline(sumStatsList, alpha = 0, priorCovariances = myU)Posterior contrasts and feature scores
Posterior means answer “how big is the effect in each context”. A different question — “which contexts does this feature deviate in, relative to its own average” — is what the contrast layer computes.
data(mashPosteriorExample)
p <- mashPosteriorExample
contrast <- mashPosteriorContrast(
p$posterior$PosteriorMean,
p$posterior$PosteriorCov,
p$orig$bhat
)
scores <- calculateFeatureScores(contrast)
scores## # A tibble: 8 × 2
## condition zScore
## <chr> <dbl>
## 1 ALL 3.41
## 2 Ast 2.58
## 3 End 2.79
## 4 Exc 3.34
## 5 Inh 3.11
## 6 Mic 3.09
## 7 OPC 1.85
## 8 Oli 2.81Each row is one condition’s meta-analysed deviation score across
features — a per-context summary of how distinctive that context’s
effects are. The bundled example is a brain cell-type panel, so the
conditions are cell types (Ast, Exc,
Mic, …) plus the ALL pooled column.
Three ways to score a feature
calculateFeatureScores() is one of three scoring routes
over the same contrast result, and they answer different questions:
# 1. meta-analysed deviation z, across features (shown above)
head(scores, 3)## # A tibble: 3 × 2
## condition zScore
## <chr> <dbl>
## 1 ALL 3.41
## 2 Ast 2.58
## 3 End 2.79
# 2. the fraction of features whose deviation is significant
nSignificantScore(contrast, pCutoff = 1e-5)## # A tibble: 8 × 2
## condition ratio
## <chr> <dbl>
## 1 ALL 0.176
## 2 Ast 0.118
## 3 End 0.118
## 4 Exc 0.176
## 5 Inh 0.118
## 6 Mic 0.176
## 7 OPC 0.118
## 8 Oli 0.118nSignificantScore() asks how often a condition
deviates rather than how strongly — a condition can carry a
large mean deviation driven by a handful of features, or a modest one
that is consistent across many, and only the second is captured
here.
The third route brings fine-mapping in. scoreFromCs()
restricts to the lead variant of each credible set, so the score
reflects deviation at variants that are plausibly causal rather than at
whatever tagged them:
scoreFromCs(p$fineMapping, contrast, condition = "Exc")## [1] 8.347042e-05It takes the least-significant credible-set lead and reports the largest pairwise |mean|/se there — a deliberately conservative summary, since one strong credible set should not carry a condition on its own.
The contrast design itself is built from
makePairwiseContrastCol(), which sets +1 and −1 for a pair
of conditions in a zero vector, and fitMashContrast(),
which fits one feature’s contrast. mashPosteriorContrast()
calls both across all features and pairs, so you only need them directly
when building a non-pairwise design.
Utilities
The mash track ships a set of helpers for getting real data into shape:
| Function | Use |
|---|---|
mashInput() |
build the (strong, random, null) matrices from
fine-mapping output |
sanitizeMashData() |
replace NaN/Inf in
bhat/sbhat with safe values |
filterInvalidSummaryStat() |
drop rows mash cannot use |
mergeMashData() / sliceMashData()
|
combine or subset the matrix bundles |
mashRandNullSample() |
draw the random/null partitions |
filterMixtureComponents() |
drop components with negligible weight |
metaAnalysisPerCondition() |
per-condition fixed-effect meta-analysis |
updateMashModelCov() |
swap a fitted model’s covariance list |
Next steps
- Building QtlSumStats objects — assembling the multi-context input.
-
TWAS weights — where a mash prior is
consumed, via
MashPriorand mr.mash.
Session info
## R version 4.5.3 (2026-03-11)
## Platform: x86_64-conda-linux-gnu
## Running under: Ubuntu 24.04.5 LTS
##
## Matrix products: default
## BLAS/LAPACK: /home/runner/work/pecotmr/pecotmr/.pixi/envs/default/lib/libopenblasp-r0.3.34.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
## [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
## [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] pecotmr_0.9.3
##
## loaded via a namespace (and not attached):
## [1] bitops_1.1-0 rlang_1.3.0
## [3] magrittr_2.0.5 horseshoe_0.2.0
## [5] otel_0.2.0 matrixStats_1.5.0
## [7] susieR_0.16.6 compiler_4.5.3
## [9] systemfonts_1.3.2 vctrs_0.7.3
## [11] quadprog_1.5-8 stringr_1.6.0
## [13] pkgconfig_2.0.3 crayon_1.5.3
## [15] fastmap_1.2.0 backports_1.5.1
## [17] XVector_0.50.0 utf8_1.2.6
## [19] Rsamtools_2.26.0 rmarkdown_2.32
## [21] tzdb_0.5.0 UCSC.utils_1.6.1
## [23] ragg_1.5.2 bit_4.6.0
## [25] purrr_1.2.2 xfun_0.61
## [27] Rfast_2.1.5.2 MultiAssayExperiment_1.36.2
## [29] cachem_1.1.0 GenomeInfoDb_1.46.2
## [31] trust_0.1-9 jsonlite_2.0.0
## [33] reshape_0.8.10 DelayedArray_0.36.1
## [35] tictoc_1.2.1 BiocParallel_1.44.0
## [37] irlba_2.4.1 parallel_4.5.3
## [39] R6_2.6.1 bslib_0.12.0
## [41] stringi_1.8.9 RColorBrewer_1.1-3
## [43] SQUAREM_2026.1 GenomicRanges_1.62.1
## [45] jquerylib_0.1.4 numDeriv_2016.8-1.1
## [47] Rcpp_1.1.2 Seqinfo_1.0.0
## [49] assertthat_0.2.1 SummarizedExperiment_1.40.0
## [51] knitr_1.52 snow_0.4-4
## [53] readr_2.2.0 IRanges_2.44.0
## [55] BiocBaseUtils_1.12.0 Matrix_1.7-6
## [57] splines_4.5.3 tidyselect_1.2.1
## [59] abind_1.4-8 yaml_2.3.12
## [61] codetools_0.2-20 metafor_5.0-1
## [63] lattice_0.23-1 tibble_3.3.1
## [65] plyr_1.8.9 withr_3.0.3
## [67] Biobase_2.70.0 S7_0.2.2
## [69] evaluate_1.0.5 flashier_1.0.21
## [71] archive_1.1.13 desc_1.4.3
## [73] RcppParallel_5.1.11-2 Biostrings_2.78.0
## [75] pillar_1.11.1 MatrixGenerics_1.22.0
## [77] metadat_1.6-0 checkmate_2.3.3
## [79] stats4_4.5.3 softImpute_1.4-3
## [81] generics_0.1.4 vroom_1.7.1
## [83] invgamma_1.2 mathjaxr_2.0-0
## [85] truncnorm_1.0-9 S4Vectors_0.48.1
## [87] hms_1.1.4 ggplot2_4.0.3
## [89] scales_1.4.0 ashr_2.2-63
## [91] glue_1.8.1 tools_4.5.3
## [93] fs_2.1.0 mvtnorm_1.4-2
## [95] grid_4.5.3 tidyr_1.3.2
## [97] nlme_3.1-170 mashr_0.2.79
## [99] deconvolveR_1.2-1 cli_3.6.6
## [101] zigg_0.0.2 ebnm_1.0-55
## [103] textshaping_1.0.5 mixsqp_0.3-54
## [105] S4Arrays_1.10.1 dplyr_1.2.1
## [107] gtable_0.3.6 sass_0.4.10
## [109] digest_0.6.39 BiocGenerics_0.56.0
## [111] SparseArray_1.10.10 htmlwidgets_1.6.4
## [113] farver_2.1.2 htmltools_0.5.9
## [115] pkgdown_2.2.1 lifecycle_1.0.5
## [117] httr_1.4.9 rmeta_3.0
## [119] bit64_4.8.6 MASS_7.3-66