---
title: "Modeling and Accounting for LD Reference Mismatch in Summary Statistics Fine-mapping"
author: "Gao Wang"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Modeling and Accounting for LD Reference Mismatch in Summary Statistics Fine-mapping}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE,comment = "#",fig.width = 5,
                      fig.height = 3,fig.align = "center",
                      fig.cap = "&nbsp;",dpi = 120)
```

This vignette gives a practical workflow for using `susie_rss` when the
GWAS summary statistics and the supplied LD from external reference samples 
mismatch, inducing high false positives in fine-mapping. This vignette complements the
[summary statistics vignette](finemapping_summary_statistics.html) and
the older [RSS diagnostic vignette](susierss_diagnostic.html).

Because an effect variant induces
marginal associations at correlated variants due to LD, SuSiE-RSS 
uses the input LD matrix to subtract the expected contribution of 
the fitted effect across the region at each single effect update. 
If that LD matrix mismatches, the subtraction is imperfect.
The leftover residual can then surface as additional sparse association
signals, creating spurious additional association signals. 
The model introduced here targets a specific mechanism by which LD mismatch 
affects SuSiE-RSS: when the LD matrix is estimated from a finite genotype reference panel, 
or from a reference panel whose population differs from the target GWAS population.

The approach here is to detect this LD mismatch pattern, quantify the
uncertainty it introduces, and attenuate signals that depend strongly on
uncertain LD subtraction. In the current `susie_rss` implementation,
`R_finite` accounts for finite LD reference panel uncertainty, and
`R_mismatch = "eb"` estimates additional LD mismatch component --- mainly motivated
as population drift --- from the data. 
A separate diagnostic flags the reliability of fine-mapped CS as well as the region 
by measuring signal attenuation and residual artifacts to caution users to inspect the results.

```{r}
library(susieR)

count_cs <- function(fit) {
  if (is.null(fit$sets) || is.null(fit$sets$cs)) return(0L)
  length(fit$sets$cs)
}
```

## Recap: a "preventable" allele coding artifact

The older
[RSS diagnostic vignette](susierss_diagnostic.html) discusses allele coding artifacts 
 as source of LD mismatch --- specifically the "allele flip". We briefly revisit the same idea here because it provides
a useful contrast with the our finite reference and population-drift LD mismatch model 
developed below.

```{r}
data(SummaryConsistency, package = "susieR")
if (!exists("SummaryConsistency"))
  load(file.path("..", "data", "SummaryConsistency.RData"))

zflip <- SummaryConsistency$z
Rtoy <- SummaryConsistency$ldref
flip_id <- SummaryConsistency$flip_id
signal_id <- SummaryConsistency$signal_id
ntoy <- 10000
btoy <- numeric(length(zflip))
btoy[signal_id] <- zflip[signal_id]
```

With the allele coding problem present, ordinary `susie_rss` reports
the true signal and an additional signal at the mismatched variant.

```{r}
fit_toy_bad <- susie_rss(zflip, Rtoy, n = ntoy)
susie_plot(fit_toy_bad, y = "PIP", b = btoy, add_legend = TRUE,
           main = "with allele coding artifact")
points(SummaryConsistency$flip_id,
       fit_toy_bad$pip[flip_id], col = 7, pch = 16)
```

The kriging diagnostic for allele flip highlights the discordant variant.

```{r}
kr_toy <- kriging_rss(zflip, Rtoy, n = ntoy)
kr_toy$plot
```

After correcting that allele encoding, the extra credible set
disappears.

```{r}
zfix <- zflip
zfix[flip_id] <- -zfix[flip_id]
fit_toy_fixed <- susie_rss(zfix, Rtoy, n = ntoy)
susie_plot(fit_toy_fixed, y = "PIP", b = btoy, add_legend = TRUE,
           main = "after allele correction")
```

In practice, it is useful to run this basic check before modeling LD
mismatch. When there is an obvious low level allele coding artifact, 
it is the best practice to fix
it directly through **proper bioinformatics procedures**. 
This can be done using many bioinformatics tools,
for example `pecotmr::allele_qc` from [this package](https://github.com/StatFunGen/pecotmr).
In this toy example, `susie_rss` reports
`r count_cs(fit_toy_bad)` credible sets before the correction and
`r count_cs(fit_toy_fixed)` after the correction, and the kriging rule
identifies the deliberately misaligned variant.


## Example using real GWAS data and a 1000 Genomes European LD reference panel


The example below from real data below shows a different situation where the
kriging outlier check does not find isolated variants to remove, but
model-based LD uncertainty still changes the fine-mapping result.
The summary statistics came from a real GWAS study of a genomic region.
The supplied LD matrix `R` was computed from a 1000 Genomes European
LD reference panel with roughly 500 samples. The full region used in
the original analysis contains more than 5,000 variants.

In the full region, ordinary `susie_rss` reported 29 credible sets,
whereas setting `R_finite = 500` reduced this to 3 credible sets, and
adding `R_mismatch = "eb"` reduced it to 1 credible set. To keep the
vignette small, the example included here uses a subset of 500 variants.
The rest of this vignette walks through this smaller example step by
step, showing how to use `R_finite`, `R_mismatch`, and the accompanying
diagnostics to evaluate LD mismatch in practice.


```{r}
data(rss_mismatch_example, package = "susieR")
if (!exists("rss_mismatch_example"))
  load(file.path("..", "data", "rss_mismatch_example.RData"))
z <- rss_mismatch_example$z
bhat <- rss_mismatch_example$bhat
shat <- rss_mismatch_example$shat
n <- rss_mismatch_example$n
R <- rss_mismatch_example$R
str(rss_mismatch_example)
```

```{r}
length(z)
dim(R)
```

The marginal z-scores show the association pattern before any
fine-mapping model is fitted.

```{r}
susie_plot(z, y = "z")
```

The marginal z-scores show one apparent signal spread over many
variants. In this subset, the smallest marginal p-value is approximately
$6.0 \times 10^{-4}$.

We define a few small helpers for the comparisons below.

```{r}
diag_scalar <- function(fit, name) {
  d <- fit$R_finite_diagnostics
  if (is.null(d) || is.null(d[[name]])) return(NA_real_)
  as.numeric(d[[name]][1])
}

diag_flag <- function(fit, name) {
  d <- fit$R_finite_diagnostics
  if (is.null(d) || is.null(d[[name]])) return(NA)
  as.logical(d[[name]][1])
}

fit_table <- function(fits) {
  data.frame(
    method = names(fits),
    n_cs = vapply(fits, count_cs, integer(1)),
    niter = vapply(fits, function(x) x$niter, integer(1)),
    converged = vapply(fits, function(x) isTRUE(x$converged), logical(1)),
    Q_art = vapply(fits, diag_scalar, numeric(1), name = "Q_art"),
    R_sensitivity_flag = vapply(fits, diag_flag, logical(1),
                                name = "R_sensitivity_flag"),
    R_reliability_flag = vapply(fits, diag_flag, logical(1),
                                name = "R_reliability_flag")
  )
}
```

## Ordinary SuSiE-RSS

First run `susie_rss` without any LD mismatch correction.

```{r}
fit_rss <- susie_rss(bhat = bhat, shat = shat,
                     R = R, n = n)
count_cs(fit_rss)
```

```{r}
susie_plot(fit_rss, y = "PIP", add_legend = TRUE,
           main = "susie_rss")
```

```{r}
fit_rss$sets$purity
cs_corr <- get_cs_correlation(fit_rss, Xcorr = R)
round(cs_corr, 3)
```

In this subset, the ordinary fit reports `r count_cs(fit_rss)`
credible sets and does not converge within `r fit_rss$niter`
iterations. The reported credible sets are tight, with
minimum purity at least
`r round(min(fit_rss$sets$purity$min.abs.corr), 3)`, but several of
them are nearly perfectly correlated with each other. This combination,
many high-purity credible sets that are mutually correlated, is a
typical sign of LD mismatch.

## Kriging diagnostic

The older RSS diagnostic checks whether individual z-scores are
surprising conditional on the other z-scores and the supplied LD matrix.
This is useful for finding isolated allele coding problems.

```{r}
kr <- kriging_rss(z, R, n = n)
kr$plot
```

The default rule labels possible allele switches when `logLR > 2` and
`abs(z) > 2`.

```{r}
outlier <- which(kr$conditional_dist$logLR > 2 &
                 abs(kr$conditional_dist$z) > 2)
length(outlier)
```

For this example, the kriging diagnostic identifies `r length(outlier)`
isolated allele switch outliers. This illustrates that 
kriging, designed for checking allele flips, does not work on this
example which is not of that type of artifacts. It is not designed to model the broader
finite LD reference error and mismatch between the LD reference panel
and the target GWAS population that can affect an entire region.

## Finite LD reference panel correction

The LD matrix here was computed from an LD reference panel of about 500
individuals. Setting `R_finite = 500` tells `susie_rss` to account for
this finite LD reference panel uncertainty in each SER update.

```{r}
fit_B500 <- susie_rss(bhat = bhat, shat = shat,
                      R = R, n = n,
                      R_finite = 500)
count_cs(fit_B500)
```

```{r}
susie_plot(fit_B500, y = "PIP", add_legend = TRUE,
           main = "R_finite = 500")
```

```{r}
vapply(fit_B500$sets$cs, length, integer(1))
fit_B500$sets$purity
cs_corr_B500 <- get_cs_correlation(fit_B500, Xcorr = R)
round(cs_corr_B500, 3)
```

In this subset of 500 variants, the finite LD reference correction
reduces the result from `r count_cs(fit_rss)` credible sets to
`r count_cs(fit_B500)`. In the original larger region, using the same
correction reduced the ordinary `susie_rss` result from 29 credible
sets to 3. The credible sets also become less artificially sharp. Their
sizes change from 2, 9, 9, 11 and 9 in the ordinary fit to 20, 3, 46 and
46, and the lowest purity drops from 0.9831 to 0.9496. The CS
correlations are less uniformly extreme, but some remain high: one pair
is still 1.000 and another is 0.978. Thus `R_finite = 500` helps, but it
does not fully resolve the mismatch pattern.

## Empirical Bayes mismatch correction

Finite LD reference uncertainty only accounts for the fact that the LD
matrix was estimated from a finite panel. It does not account for
systematic mismatch due to drift between the LD reference population and the GWAS
population. The option `R_mismatch = "eb"` estimates an additional
variance component shared across the region from the fitted residuals.

```{r}
fit_B500_eb <- susie_rss(bhat = bhat, shat = shat,
                         R = R, n = n,
                         R_finite = 500, R_mismatch = "eb")
count_cs(fit_B500_eb)
```

```{r}
susie_plot(fit_B500_eb, y = "PIP", add_legend = TRUE,
           main = "R_finite = 500, R_mismatch = eb")
```

Here, adding `R_mismatch = "eb"` reduces the result to
`r count_cs(fit_B500_eb)` credible set and the fit converges.

One can also run the mismatch model without specifying a finite
LD reference panel size. This treats the continuous mismatch variance as
the only LD uncertainty component.

```{r}
fit_Binf_eb <- susie_rss(bhat = bhat, shat = shat,
                         R = R, n = n,
                         R_mismatch = "eb")
count_cs(fit_Binf_eb)
```

```{r}
susie_plot(fit_Binf_eb, y = "PIP", add_legend = TRUE,
           main = "R_mismatch = eb")
```

This fit also reports `r count_cs(fit_Binf_eb)` credible set and
converges quickly.

```{r}
vapply(fit_B500_eb$sets$cs, length, integer(1))
fit_B500_eb$sets$purity
max(fit_B500_eb$pip)
```

In this subset, the two EB mismatch fits report the same 95% credible
set. It contains 34 variants, has minimum purity 0.8609, and has
coverage 0.9543; the largest PIP is 0.1142. Compared with the earlier
fits, the EB adjustment keeps one signal but represents it as a broad
credible set rather than several small, highly correlated credible sets.
This is more consistent with the weak marginal association in this
region. 

In practice we recommend setting both `R_finite` and `R_mismatch` unless your
reference LD is based on large samples such as current UK Biobank (300,000 to 500,000)
in which case setting e.g. `R_finite = 300000` will have no practical impact to the results.

## Diagnostic metrics

The warnings above should not be ignored. Even after the EB adjustment,
`susie_rss` reports diagnostic warnings, indicating that the fine
mapping result remain uncertain. The fit object returns several
diagnostic metrics to summarize this uncertainty.

```{r}
fits <- list(susie_rss = fit_rss,
             R_finite_500 = fit_B500,
             R_finite_500_eb = fit_B500_eb,
             R_mismatch_eb = fit_Binf_eb)
knitr::kable(fit_table(fits), digits = 3)
```

There are two relevant reliability diagnostics. `Q_art` measures
whether the final residual points into directions where the supplied LD
matrix has little support. `R_sensitivity_flag` measures whether reported
credible sets depend strongly on the LD uncertainty penalty. If
either diagnostic is triggered, `R_reliability_flag` is `TRUE`, and
the fit should be interpreted cautiously.

In this subset, `R_finite = 500` is flagged by the credible set
sensitivity diagnostic. The two EB mismatch fits are flagged by the
residual space diagnostic, even though they report only one credible
set. This is why the warning points users to
`fit$R_finite_diagnostics` rather than treating the corrected fit as
automatically reliable.

## Diagnostic with one effect

When the reliability flag is raised, the fit object provides a
diagnostic model using single-effect regression (SER):

```{r}
ser <- fit_Binf_eb$R_finite_diagnostics$ser_model
names(ser)
```

This is a credible set model with one effect as in Maller et al.
(2012). For this model, the single row of `alpha` is the PIP.

```{r}
identical(as.numeric(ser$alpha[1, ]), as.numeric(ser$pip))
count_cs(ser)
vapply(ser$sets$cs, length, integer(1))
ser$sets$purity
max(ser$pip)
identical(ser$sets$cs, fit_Binf_eb$sets$cs)
isTRUE(all.equal(as.numeric(ser$pip), as.numeric(fit_Binf_eb$pip)))
```

```{r}
susie_plot(ser, y = "PIP", add_legend = TRUE,
           main = "diagnostic with one effect")
```

In this example, the one-effect diagnostic gives the same result as the
EB mismatch fit above: one 34-variant credible set, minimum purity
0.8609, and maximum PIP 0.1142. This exact match is merely a coincidence, not
a general property of the method. When an EB mismatch fit produces a
reliability warning, compare the EB result with this one-effect
diagnostic rather than treating the EB fit as automatically reliable.

## Summary

For routine analyses with GWAS summary statistics and an external LD
reference panel:

- First run basic allele and variant alignment quality control and checks. 
  The kriging
  diagnostic can be useful for such obvious, localized allele coding artifacts.
- Fit ordinary `susie_rss` and look for signs of instability, such as
  nonconvergence, multiple highly correlated credible sets, or many
  small credible sets in a weak association region.
- Use `R_finite` when the LD matrix was estimated from a finite LD
  reference panel; add `R_mismatch = "eb"` when broader LD mismatch remains after
  accounting for finite reference uncertainty.
- Inspect `fit$R_finite_diagnostics`, including `Q_art`,
  `R_sensitivity_flag`, and `R_reliability_flag`.
- If a reliability warning is raised, compare the EB result with the
  SER diagnostic model, and interpret the fine-mapping
  result conservatively with the SER model as a safer fall-back.

## Caution: allele coding artifacts are outside this model

The old RSS vignette example is based on the out-of-sample LD reference example
from Zou et al. (2022). Because its LD matrix was computed from a finite
reference panel of 500 individuals, it is appropriate to set
`R_finite = 500` when using the new diagnostics. This setting can raise
a reliability warning, but that warning is not an allele flip correction
and it does not replace allele alignment QC.

```{r, message=FALSE, warning=FALSE}
fit_toy_bad_B500_eb <- susie_rss(zflip, Rtoy, n = ntoy,
                                 R_finite = 500,
                                 R_mismatch = "eb",
                                 track_fit = TRUE)
d_toy <- fit_toy_bad_B500_eb$R_finite_diagnostics

knitr::kable(data.frame(
  n_cs = count_cs(fit_toy_bad_B500_eb),
  flipped_variant_pip =
    unname(fit_toy_bad_B500_eb$pip[SummaryConsistency$flip_id]),
  Q_art = d_toy$Q_art,
  residual_artifact_flag = d_toy$artifact_flag,
  CS_sensitivity_flag = d_toy$R_sensitivity_flag,
  overall_reliability_flag = d_toy$R_reliability_flag
), digits = 3, row.names = FALSE)
```

```{r}
susie_plot(fit_toy_bad_B500_eb, y = "PIP", b = btoy, add_legend = TRUE,
           main = "allele artifact with R_finite = 500, R_mismatch = eb")
points(SummaryConsistency$flip_id,
       fit_toy_bad_B500_eb$pip[SummaryConsistency$flip_id],
       col = 7, pch = 16)
```

This fit still reports 2 credible sets, and the flipped variant still
has PIP close to 1. The residual artifact diagnostic does not flag this
fit: `Q_art` is small and `residual_artifact_flag` is `FALSE`. However,
the CS-level sensitivity diagnostic is `TRUE`, so the overall
reliability flag is also `TRUE`. The warning therefore comes from CS
sensitivity, not from `Q_art`.

This behavior is expected. The model in this vignette is designed for
finite reference error and population drift in the LD matrix. It does
not include an allele flip parameter, and it does not try to decide
whether an individual z-score should have its sign reversed. That is
deliberate: allele alignment is a data-processing problem that should
be handled before fitting the LD mismatch model.

The reason is visible from the z-scores. Given the true signal and the
LD matrix, the flipped variant should have a positive marginal
association. Instead its observed z-score has the opposite sign.

```{r}
flip_expected <- Rtoy[SummaryConsistency$flip_id,
                      SummaryConsistency$signal_id] *
                 zflip[SummaryConsistency$signal_id]
knitr::kable(data.frame(
  flipped_variant_z = zflip[SummaryConsistency$flip_id],
  LD_implied_z_from_signal = flip_expected,
  residual_after_signal = zflip[SummaryConsistency$flip_id] - flip_expected
), digits = 3)
```

This is a deterministic sign error at one variant, not the stochastic LD
uncertainty modeled by `R_finite` or `R_mismatch = "eb"`.

The compact track history can then be used to follow this specific
artifact. We do not need to inspect the full object. For this toy
example, the relevant rows are the single effects whose top variable is
either the true signal or the known flipped variant. After the first
IBSS iteration, one effect is already assigned almost entirely to the
flipped variant.

```{r}
track_toy <- fit_toy_bad_B500_eb$trace
track_effect <- subset(track_toy$effect,
                       iteration %in% 1:3 &
                       top_variable %in% c(SummaryConsistency$signal_id,
                                           SummaryConsistency$flip_id))
track_story <- data.frame(
  iteration = track_effect$iteration,
  effect = track_effect$effect,
  variant = paste0(track_effect$top_variable_name, " (", ifelse(
    track_effect$top_variable == SummaryConsistency$signal_id,
    "true signal",
    "flipped variant"
  ), ")"),
  alpha_at_variant = track_effect$top_alpha,
  max_variable_lbf = track_effect$top_lbf,
  single_effect_lbf = track_effect$effect_lbf
)

knitr::kable(track_story,
             digits = 3)
```

Because this is a toy example, we know which variant is the true signal
and which variant was deliberately flipped. The track table itself only
shows the fitted single effect summaries: one effect is already highly
concentrated on the flipped variant after the first iteration.

The CS-level diagnostic also shows why this happens. The fitted log BF
and the log BF attenuation are different quantities: the former
describes the fitted signal strength, whereas the latter measures how
much that signal depends on uncertain LD.

```{r}
cs_story <- track_toy$cs_sensitivity
final_effect <- track_toy$effect[
  track_toy$effect$iteration == max(track_toy$effect$iteration),
  c("effect", "top_lbf")
]
cs_story <- merge(cs_story, final_effect, by = "effect",
                  all.x = TRUE, sort = FALSE)
cs_story <- cs_story[cs_story$is_sensitive,
                     c("cs", "effect", "cs_size", "top_lbf",
                       "cs_max_bf_attenuation", "threshold",
                       "cs_sensitivity_label",
                       "top_sensitive_variable_name")]
names(cs_story) <- c("cs", "effect", "cs_size", "max_log_BF",
                     "log_BF_attenuation", "threshold", "sensitivity",
                     "sensitive_variant")
knitr::kable(cs_story,
             digits = 3)
```

For both credible sets, the attenuation is above the threshold
`log(20) = 3.0`, so both are labeled `highly_sensitive`. This is why
`CS_sensitivity_flag` is `TRUE`, and why the overall reliability flag
is `TRUE`, even though `Q_art = 0.012` and the residual artifact flag is
`FALSE`.

Thus, an allele flip can create a convincing sparse effect before the
LD mismatch model has anything useful to adjust. In this example,
`R_finite = 500` still raises a reliability warning through CS-level
sensitivity, but that warning should be interpreted as a cue to inspect
alignment, not as an automatic correction of the allele coding problem.
In practice this warning is not guaranteed: for example, fitting this
same toy example with `R_mismatch = "eb"` but without `R_finite = 500`
does not raise the reliability warning.

```{r cleanup, include=FALSE}
invisible(NULL)
```

## Session information

```{r session-info}
sessionInfo()
```

## Reference

Maller 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).
