susie_rss fine-maps a
region from summary statistics z-scores and an LD matrix. When the LD
matrix does not match the summary statistics, the fit can report extra,
spurious credible sets (CSs). This vignette walks through the workflow
we recommend for handling this, step by step, on simulated data:
susie_rss function by
setting R_finite = B and
R_mismatch = "eb_mix", where B is the sample
size of your LD reference panel. This accounts for an LD matrix
estimated from an external reference panel, drawn from a population that
can differ from the summary statistics population.SER_fallback. If it is
TRUE, the corrected fit is still not reliable; report the
one-effect (SER) model instead.susie_rss a
list of reference panels. It estimates how much weight
each panel gets, then fine-maps with the weighted (mixed) LD.All examples are simulated summary statistics for regions with one causal variant. The first four cases simulate a phenotype from GTEx genotypes (n = 838). The multi-panel case simulates phenotypes for an African-ancestry cohort and a European-ancestry cohort and meta-analyzes them.
| Case | Step | Summary statistics | LD reference | Issue |
|---|---|---|---|---|
allele_flip |
1 | GTEx, n = 838 | in-sample GTEx LD | 1% of the variants have the sign of their z-score reversed (2 in this region) |
ld_mismatch |
2 | GTEx, n = 838 | 1000 Genomes EUR LD, B = 489 | external LD panel |
ld_mismatch_sketch |
2 | GTEx, n = 838 | ADSP sketch, B = 500 | external LD panel, given as a genotype sketch instead of an LD matrix |
ser_fallback |
3 | GTEx, n = 838 | 1000 Genomes EUR LD, B = 489 | external LD panel |
multi_panel |
4 | 1000 Genomes AFR (n = 400) + UK Biobank unrelated EUR (n = 400), total n = 800 | AFR and EUR panel LD | no single panel matches the mixed summary statistics |
The inputs for all five cases are packed in one list. Each element
holds what susie_rss needs:
z: the summary statistics z-scores;R: the LD matrix (a named list of LD matrices for
multi_panel), or X: a factor matrix (genotype
sketch) of the reference panel, used instead of R;n: the summary statistics sample size;B: the reference panel sample size, used as
R_finite.Because the data are simulated, each element also stores
causal, the index of the true causal variant. The
external-reference cases also store z_insample and
R_insample: the z-scores and the in-sample LD on the same
variants. allele_flip stores z_clean (the
z-scores before flipping) and flip_id (the indices of the
flipped variants); multi_panel stores prop,
the true proportion of summary statistics samples from each
ancestry.
A few small helpers for the comparisons below.
fit_table() puts several fits side by side, including how
many CSs contain the causal variant; cs_summary() lists the
CSs of one fit; plot_region() draws PIPs only in the part
of the region around the CSs and the causal variant.
cs_summary <- function(fit, causal = NULL) {
cs <- fit$sets$cs
if (length(cs) == 0) return(data.frame(note = "no credible sets"))
data.frame(
cs = names(cs),
size = lengths(cs),
purity = fit$sets$purity$min.abs.corr,
top_variant = vapply(cs, function(v) v[which.max(fit$pip[v])], numeric(1)),
max_pip = vapply(cs, function(v) max(fit$pip[v]), numeric(1)),
has_causal = if (is.null(causal)) NA else
vapply(cs, function(v) any(v %in% causal), logical(1)),
row.names = NULL
)
}
diag_value <- function(fit, name) {
d <- fit$R_finite_diagnostics
if (is.null(d) || is.null(d[[name]])) NA else d[[name]][1]
}
fit_row <- function(fit, causal = NULL) {
if (is.null(fit))
return(data.frame(n_cs = NA, n_cs_with_causal = NA, converged = NA,
Q_art = NA, R_sensitivity_flag = NA, SER_fallback = NA))
cs <- fit$sets$cs
data.frame(
n_cs = length(cs),
n_cs_with_causal = if (is.null(causal)) NA else
sum(vapply(cs, function(v) any(v %in% causal), logical(1))),
converged = if (is.null(fit$converged)) NA else as.logical(fit$converged),
Q_art = as.numeric(diag_value(fit, "Q_art")),
R_sensitivity_flag = as.logical(diag_value(fit, "R_sensitivity_flag")),
SER_fallback = as.logical(diag_value(fit, "SER_fallback"))
)
}
fit_table <- function(fits, causal = NULL) {
out <- do.call(rbind, lapply(fits, fit_row, causal = causal))
out <- data.frame(method = names(fits), out, row.names = NULL)
out[, colSums(!is.na(out)) > 0, drop = FALSE]
}
kriging_outliers <- function(kr) {
which(kr$conditional_dist$logLR > 2 & abs(kr$conditional_dist$z) > 2)
}
plot_region <- function(fit, x, main, pad = 30) {
idx <- c(unlist(fit$sets$cs), x$causal)
if (length(idx) == 0) idx <- which.max(fit$pip)
pos <- max(1, min(idx) - pad):min(length(fit$pip), max(idx) + pad)
susie_plot(fit, y = "PIP", pos = pos, add_legend = TRUE, main = main)
# Mark the true causal variant in red. Passing b = to susie_plot()
# misplaces this marker when pos covers only part of the region.
points(x$causal, fit$pip[x$causal], col = 2, pch = 16)
}In this region the LD matrix is in-sample: it was computed from the same 838 GTEx individuals as the summary statistics. The sign of the z-score was then reversed at 1% of the variants, as happens when the effect allele is coded differently in the summary statistics and the LD reference. Since the LD itself is correct, any problem here comes from allele flips.
Before running kriging, we strongly encourage checking allele consistency between the summary statistics and the LD reference, for example with pecotmr. If the inconsistency is widespread, for example because of different variant naming systems, statistical methods can fail arbitrarily badly. Kriging is a check for any flipped variants left over after that QC.
The kriging diagnostic predicts each z-score from all the others using the LD matrix. A variant whose observed z-score has the opposite sign of its prediction is flagged as a likely allele flip and drawn in red. See the RSS diagnostic vignette for details.
out_flip <- kriging_outliers(kr_flip)
c(flagged = length(out_flip),
truly_flipped = sum(out_flip %in% af$flip_id),
total_flipped = length(af$flip_id))
# flagged truly_flipped total_flipped
# 1 1 2Kriging flags 1 variant, and it is a true flip. The other flip is at a variant with a small z-score, where a sign error changes little, so it matters much less for fine-mapping.
For reference, the same check on the unflipped z-scores, with the same in-sample LD, flags no variants:
We reverse the sign of the flagged variant and fine-map with the original, flipped and corrected z-scores.
z_fix <- af$z
z_fix[out_flip] <- -z_fix[out_flip]
fit_flip <- susie_rss(z = af$z, R = af$R, n = af$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# WARNING: IBSS algorithm did not converge in 50 iterations!
# Warning: IBSS algorithm did not converge in 50 iterations!
fit_fix <- susie_rss(z = z_fix, R = af$R, n = af$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
fit_clean <- susie_rss(z = af$z_clean, R = af$R, n = af$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
knitr::kable(fit_table(list(original = fit_clean,
flipped = fit_flip,
kriging_corrected = fit_fix),
causal = af$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged |
|---|---|---|---|
| original | 1 | 1 | TRUE |
| flipped | 5 | 0 | FALSE |
| kriging_corrected | 1 | 1 | TRUE |
With the flipped z-scores, susie_rss reports 5 CSs, none
containing the causal variant, and does not converge. After the
correction it reports 1 CS containing the causal variant, the same as
with the original z-scores.
In this case the alleles are aligned, but the LD matrix comes from an external panel (489 Europeans from 1000 Genomes), rather than from the 838 GTEx individuals in the summary statistics. GTEx samples are of mixed ancestry, so this LD differs from the LD of the summary statistics samples. Running kriging first finds nothing to correct:
However, fine-mapping with the external reference is still misled:
fit_in <- susie_rss(z = lm$z_insample, R = lm$R_insample, n = lm$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
fit_lm <- susie_rss(z = lm$z, R = lm$R, n = lm$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
knitr::kable(fit_table(list(in_sample_LD = fit_in, ordinary = fit_lm),
causal = lm$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged |
|---|---|---|---|
| in_sample_LD | 1 | 1 | TRUE |
| ordinary | 2 | 1 | TRUE |
With the 1000 Genomes LD, ordinary susie_rss reports 2
CSs, only one of which contains the causal variant. With the in-sample
LD there is a single CS.
R_finite and R_mismatchSuSiE-rRSS accounts for this with two options of
susie_rss:
R_finite = B tells susie_rss that the LD
was estimated from B reference individuals, here 489, and
inflates the uncertainty of each update accordingly.R_mismatch = "eb_mix" estimates from the data on how
much additional LD mismatch the region has. It then applies that extra
uncertainty to each variant in proportion to how likely the variant’s
residual is to be mismatch rather than true signal.fit_lm_rrss <- susie_rss(z = lm$z, R = lm$R, n = lm$n,
R_finite = lm$B, R_mismatch = "eb_mix")
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# HINT: Switching to PIP-based convergence: R uncertainty inflation modifies per-variant SER likelihoods, which prevents a consistent model-level ELBO.
knitr::kable(fit_table(list(in_sample_LD = fit_in, ordinary = fit_lm,
`SuSiE-rRSS` = fit_lm_rrss),
causal = lm$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged | Q_art | R_sensitivity_flag | SER_fallback |
|---|---|---|---|---|---|---|
| in_sample_LD | 1 | 1 | TRUE | NA | NA | NA |
| ordinary | 2 | 1 | TRUE | NA | NA | NA |
| SuSiE-rRSS | 1 | 1 | TRUE | 0.011 | FALSE | FALSE |
With SuSiE-rRSS (R_finite = 489,
R_mismatch = "eb_mix"), the fit reports 1 CS and drops the
spurious second CS of the ordinary fit. SER_fallback
(explained in Step 3) is FALSE, so the SuSiE-rRSS fit is
the result to report.
Instead of a LD matrix, the reference panel can be given as a factor
matrix (sketch) X with p columns (p = number of variants).
The workflow is unchanged: pass X instead of
R.
Here the reference is a sketch of the ADSP panel with 500 rows:
Kriging again flags nothing. For a sketch, R_finite is
the sketch dimension, the number of rows of X, here
500.
fit_sk_in <- susie_rss(z = sk$z_insample, R = sk$R_insample, n = sk$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
fit_sk <- susie_rss(z = sk$z, X = sk$X, n = sk$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# WARNING: IBSS algorithm did not converge in 50 iterations!
# Warning: IBSS algorithm did not converge in 50 iterations!
fit_sk_rrss <- susie_rss(z = sk$z, X = sk$X, n = sk$n,
R_finite = sk$B, R_mismatch = "eb_mix")
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# HINT: Switching to PIP-based convergence: R uncertainty inflation modifies per-variant SER likelihoods, which prevents a consistent model-level ELBO.
knitr::kable(fit_table(list(in_sample_LD = fit_sk_in, ordinary = fit_sk,
`SuSiE-rRSS` = fit_sk_rrss),
causal = sk$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged | Q_art | R_sensitivity_flag | SER_fallback |
|---|---|---|---|---|---|---|
| in_sample_LD | 1 | 1 | TRUE | NA | NA | NA |
| ordinary | 4 | 1 | FALSE | NA | NA | NA |
| SuSiE-rRSS | 1 | 1 | TRUE | 0.027 | FALSE | FALSE |
With the sketch, ordinary susie_rss reports 4 CSs and
does not converge. With SuSiE-rRSS (R_finite = 500,
R_mismatch = "eb_mix"), it reports 1 CS with a single
variant, and SER_fallback is FALSE.
SER_fallback
points to the SER modelThis region uses the same 1000 Genomes reference (B = 489) on another simulated dataset. As before, start with kriging:
Kriging flags nothing, so we fit susie_rss with the
in-sample LD, the 1000 Genomes LD, and the 1000 Genomes LD with
SuSiE-rRSS:
fit_sf_in <- susie_rss(z = sf$z_insample, R = sf$R_insample, n = sf$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
fit_sf <- susie_rss(z = sf$z, R = sf$R, n = sf$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# WARNING: IBSS algorithm did not converge in 50 iterations!
# Warning: IBSS algorithm did not converge in 50 iterations!
fit_sf_rrss <- susie_rss(z = sf$z, R = sf$R, n = sf$n,
R_finite = sf$B, R_mismatch = "eb_mix")
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# HINT: Switching to PIP-based convergence: R uncertainty inflation modifies per-variant SER likelihoods, which prevents a consistent model-level ELBO.
# WARNING: IBSS algorithm did not converge in 50 iterations!
# Warning: IBSS algorithm did not converge in 50 iterations!
# WARNING: Possible summary statistics and R reference mismatch detected. Fine-mapping results may be unreliable; inspect fit$R_finite_diagnostics for details. A one-effect diagnostic is available at fit$R_finite_diagnostics$ser_model; it is a one-effect credible-set model as in Maller et al. 2012, and its alpha row is the PIP.
# Warning: Possible summary statistics and R reference mismatch detected.
# Fine-mapping results may be unreliable; inspect fit$R_finite_diagnostics for
# details. A one-effect diagnostic is available at
# fit$R_finite_diagnostics$ser_model; it is a one-effect credible-set model as in
# Maller et al. 2012, and its alpha row is the PIP.
knitr::kable(fit_table(list(in_sample_LD = fit_sf_in, ordinary = fit_sf,
`SuSiE-rRSS` = fit_sf_rrss),
causal = sf$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged | Q_art | R_sensitivity_flag | SER_fallback |
|---|---|---|---|---|---|---|
| in_sample_LD | 1 | 1 | TRUE | NA | NA | NA |
| ordinary | 4 | 0 | FALSE | NA | NA | NA |
| SuSiE-rRSS | 4 | 0 | FALSE | 0.082 | TRUE | TRUE |
The ordinary fit reports 4 CSs, none containing the causal variant,
and does not converge. The SuSiE-rRSS fit does not converge either and
still reports 4 CSs, none of which contain the causal variant.
susie_rss prints a warning that the result may be
unreliable: SER_fallback is TRUE.
SER_fallback checksSER_fallback is TRUE when
either of two checks fails:
artifact_flag):
Q_art is the share of the fitted residual that lies in
directions the LD matrix barely supports. A value above 0.1 means part
of the residual cannot be explained by the LD.R_sensitivity_flag): for each CS, how much its log Bayes
factor drops once LD uncertainty is accounted for. A drop above 30 means
the CS depends on uncertain LD subtraction, which is how imperfect
conditioning creates spurious signals.round(fit_sf_rrss$R_finite_diagnostics$bf_attenuation$cs_max, 2)
# L2 L5 L3 L4
# 1198.61 467.01 875.20 742.06Q_art is 0.082, below 0.1, so the residual check passes.
But every CS loses far more than 30 in log Bayes factor, so the
sensitivity check fails.
When SER_fallback is TRUE, the fit provides
a one-effect model in fit$R_finite_diagnostics$ser_model.
It is a single-effect regression (SER), as in Maller et al. (2012): it
assumes one causal variant in the region, so it cannot report extra CSs
created by LD mismatch. Its PIPs are the single row of
alpha. Using the SER model helps avoid false
discoveries.
Conclusion: report the full fit from step 2 unless
SER_fallback is TRUE, in which case report the
SER model.
final_fit <- function(fit) {
if (isTRUE(fit$R_finite_diagnostics$SER_fallback))
fit$R_finite_diagnostics$ser_model
else fit
}
knitr::kable(fit_table(list(in_sample_LD = fit_sf_in, ordinary = fit_sf,
`SuSiE-rRSS` = fit_sf_rrss,
SER_model = final_fit(fit_sf_rrss)),
causal = sf$causal), digits = 3)| method | n_cs | n_cs_with_causal | converged | Q_art | R_sensitivity_flag | SER_fallback |
|---|---|---|---|---|---|---|
| in_sample_LD | 1 | 1 | TRUE | NA | NA | NA |
| ordinary | 4 | 0 | FALSE | NA | NA | NA |
| SuSiE-rRSS | 4 | 0 | FALSE | 0.082 | TRUE | TRUE |
| SER_model | 1 | 1 | NA | NA | NA | NA |
The SER model reports 1 CS of 3 variants, which contains the causal variant, with a top PIP of 0.33.
These summary statistics are a meta-analysis of an African-ancestry cohort (400 samples) and a European-ancestry cohort of unrelated UK Biobank samples (400 samples), so each ancestry makes up 50% of the 800 samples. Neither the AFR nor the EUR reference panel alone matches the LD of this mixed sample.
Pass both panels as a named list. susie_rss then works
in two stages:
omega of each panel, so that
the mixed LD omega_AFR * R_AFR + omega_EUR * R_EUR best
fits the z-scores.mp <- cases$multi_panel
fit_mp <- susie_rss(z = mp$z, R = mp$R, n = mp$n)
# HINT: Setting max_iter = 50 for the SuSiE RSS model because slow convergence is often a sign of unstable summary-statistics fitting. To disable this message, explicitly set max_iter = 50 or another value in the susie_rss() call.
# HINT: Switching to PIP-based convergence: multi-panel mixture weight updates change R(omega) each iteration, which prevents ELBO monotonicity.
knitr::kable(data.frame(panel = names(mp$R),
estimated_weight = fit_mp$omega_weights,
true_proportion = mp$prop,
row.names = NULL), digits = 3)| panel | estimated_weight | true_proportion |
|---|---|---|
| AFR | 0.409 | 0.5 |
| EUR | 0.591 | 0.5 |
| cs | size | purity | top_variant | max_pip | has_causal |
|---|---|---|---|---|---|
| L1 | 1 | 1 | 61 | 0.962 | TRUE |
The estimated weights, 0.409 for AFR and 0.591 for EUR, are close to the true proportions of 0.5 and 0.5. The fit reports 1 CS: a single variant, which is the causal variant.
For comparison, the fit also keeps one fit per panel, which is used to choose where the weight estimation starts:
spf <- fit_mp$single_panel_fits
names(spf) <- paste0(names(mp$R), "_only")
knitr::kable(fit_table(c(spf, list(mixed = fit_mp)), causal = mp$causal),
digits = 3)| method | n_cs | n_cs_with_causal | converged |
|---|---|---|---|
| AFR_only | 5 | 1 | TRUE |
| EUR_only | 3 | 1 | FALSE |
| mixed | 1 | 1 | TRUE |
With a single panel, the fit reports extra CSs: 5 with the AFR panel and 3 with the EUR panel (which does not converge), and in each case only one of them contains the causal variant. With the mixed LD there is a single CS, the causal variant.
If the panels are given as factor matrices instead of LD matrices,
pass X = list(AFR = X_AFR, EUR = X_EUR).
For fine-mapping with summary statistics and external LD:
pecotmr::summaryStatsQc().susie_rss with
R_finite and R_mismatch = "eb_mix") when the
LD comes from a reference panel, given either as an LD matrix
R or as a genotype sketch X.fit$R_finite_diagnostics$SER_fallback. If it is
TRUE, report
fit$R_finite_diagnostics$ser_model instead of the full
fit.susie_rss estimates the panel weights
(fit$omega_weights) and fine-maps with the mixed LD. This
can be combined with R_finite and R_mismatch
as above.sessionInfo()
# R version 4.6.1 (2026-06-24)
# Platform: x86_64-pc-linux-gnu
# Running under: Ubuntu 26.04.1 LTS
#
# Matrix products: default
# BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
# LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#
# locale:
# [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
# [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
# [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
# [9] LC_ADDRESS=C LC_TELEPHONE=C
# [11] LC_MEASUREMENT=en_US.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] L0Learn_2.1.0 susieR_0.16.6 rmarkdown_2.32
#
# loaded via a namespace (and not attached):
# [1] Matrix_1.7-6 gtable_0.3.6 jsonlite_2.0.0 compiler_4.6.1
# [5] crayon_1.5.3 Rcpp_1.1.2 stringr_1.6.0 jquerylib_0.1.4
# [9] scales_1.4.0 yaml_2.3.12 fastmap_1.2.0 lattice_0.23-1
# [13] plyr_1.8.9 ggplot2_4.0.3 R6_2.6.1 labeling_0.4.3
# [17] mixsqp_0.3-54 knitr_1.52 MASS_7.3-66 maketools_1.3.2
# [21] zigg_0.0.2 bslib_0.12.0 RColorBrewer_1.1-3 rlang_1.3.0
# [25] stringi_1.8.9 cachem_1.1.0 reshape_0.8.10 xfun_0.61
# [29] sass_0.4.10 sys_3.4.3 S7_0.2.2 RcppParallel_6.2.1
# [33] otel_0.2.0 cli_3.6.6 withr_3.0.3 magrittr_2.0.5
# [37] digest_0.6.39 grid_4.6.1 irlba_2.4.1 lifecycle_1.0.5
# [41] vctrs_0.7.3 Rfast_2.1.5.2 evaluate_1.0.5 glue_1.8.1
# [45] farver_2.1.2 buildtools_1.0.0 reshape2_1.4.5 matrixStats_1.5.0
# [49] tools_4.6.1 htmltools_0.5.9Maller JB, McVean G, Byrnes J, Vukcevic D, Palin K, Su Z, et al. Bayesian refinement of association signals for 14 loci in 3 common diseases. Nature Genetics 44, 1294–1301 (2012).