Skip to contents

Problem Statement

Source: vignettes/articles/problem-statement.Rmd
problem-statement.Rmd

Current RAS algorithm:

\[LPRS_i = \sum_{j=1}^{m} X_{ij} \hat{\beta}_j\] Where in a specific window of \(m\) SNPs, it calculates a Localized Polygenic Risk Score (LPRS) for each individual \(i\). \(X_{ij}\) is the genotype dosage of individual \(i\) at SNP \(j\), and \(\hat{\beta}_j\) is the marginal effect size of SNP \(j\) estimated from the training split.

It then tests the joint significance by taking the vector of LPRS values outputted and runs a regression against the vector of actual phenotypes. The p-value is calculated to determine if the slope of this regression line is statistically significant.

Summary Statistic Equivalent:

Given vector of marginal Z-scores for the \(m\) SNPs in the window: \[Z = (Z_1, Z_2, ..., Z_m)^T\]

To make sure \(Z\) is a valid replacement for the original individual data, summary statistics must be generated from target-cohort GWAS that was adjusted for the same set of covariates (age, sex, etc.) used in the original model (\(y \sim \text{LPRS} + X\)).

Let \(R\) be the \(m \times m\) Linkage Disequilibrium (LD) correlation matrix of these SNPs, which can be estimated from an external reference (e.g., 1000 Genomes). \(R\) must be ancestry-matched.

Let \(w = (w_1, w_2, ..., w_m)^T\) be vector of weights derived from training data. Similar to original algorithm, \(w\) contains the SNP marginal effect sizes (\(\hat{\beta}\)) from standard GWAS procedures ran on the training split.

Under null hypothesis (\(H_0\)) of the SNP window having absolutely no effect on phenotype, a multivariate normal distribution can be created from this Z-score vector (assuming large enough target-cohort \(N\) for this approximation to hold). \[Z \sim N(0, R)\]

The Test Statistic:

Current ras_scan algorithm needs single p-value for whole window, not m different Z-scores, so these m dimensions have to be compressed down to a single number.

Instead of a multi-degree-of-freedom test, the marginal Z-scores are projected onto the weight vector \(w\) to form a 1 degree-of-freedom localized weighted burden statistic. Instead of estimating joint effects (already in \(w\)), the correlation matrix \(R\) is used only to account for linkage disequilibrium (LD) when evaluating the null distribution of the weighted burden statistic.

Equation: \[T_{\text{burden}} = \frac{w^T Z}{\sqrt{w^T R w}}\]

  • \(T_{\text{burden}}\) is the calculated 1 d.f. burden test statistic
  • \(w^T Z\) is the dot product of the weights and the marginal Z-scores (the summary-stat version of the LPRS calculation)
  • \(\sqrt{w^T R w}\) is the std deviation (sqrt of the Quadratic Form) of this weighted sum, adjusting for correlation structure defined by the LD matrix \(R\)

With \(T_{\text{burden}}\) calculated, under the null hypothesis it follows a standard normal distribution: \[T_{\text{burden}} \sim N(0, 1)\]

The two-tailed window p-value can be calculated analytically with the standard normal cumulative distribution function (\(\Phi\)): \[p = 2(1 - \Phi(|T_{\text{burden}}|))\]

The RAS is simply \(-log_{10}(p-value)\).

Quality Control/Regularization:

  1. Quadratic Form Stability: Unlike methods that require matrix inversion, the denominator \(w^T R w\) is a quadratic form. It is numerically stable even if the LD matrix \(R\) is singular or near-singular. Therefore, no ridge regularization is required.
  2. LD Pruning (Optional): While not required for numerical stability, LD pruning at \(r^2 < 0.2\) can optionally be applied to prevent highly redundant variants from over-weighting the burden statistic, replicating the original RAS design.

Weight Source Decision:

In the original individual algorithm, a 50/50 data-splitting protocol provides an independent training set to fit weights and a testing set to evaluate them, eliminating overfitting. Since summary-statistic algorithms don’t have those individual-level splits, must choose between two methods for defining \(w\):

  • Method A: Derive the weight vector \(w\) from an independent, external discovery summary-statistic dataset (e.g., a large public study) and test it against the target cohort’s Z-scores. This is statistically rigorous because it avoids overfitting winner’s curse bias.
  • Method B: If only a single summary-statistic dataset is available, deriving \(w\) from the same marginal effects we are testing (\(w = \hat{\beta}\)) introduces severe winner’s curse bias, inflating the burden statistic. More practical for rarer diseases/phenotypes which don’t have much data available publicly.

While a joint-effect approach with LD reconstruction (e.g., SuSiE-RSS) is a possible alternative, implementing it goes beyond the scope of replicating the original methodology

Validation:

  1. Type I error test: null simulation of the full algorithm (pruning and max-over-\(t\) window selection), reporting empirical type I error, similar to Table 1 of the original paper.
  2. Power test: power simulation of full algorithm, reporting statistical power for comparison
  3. Convergence validation: compare summary-stat RAS output against individual-data RAS on the same simulated data (completed below), varying \(N\) and LD-panel quality (yet to be done)

Validation Results

Setup: 300-SNP toy chromosome, AR(1) LD at \(\rho = 0.8\), true signal at 3 consecutive causal SNPs (150–152, effect size 0.5). Discovery and target cohorts of \(n = 1000\) each (Method A), LD pruned at \(r^2 < 0.2\) (300 -> 100 SNPs survived). The validated workflow is LD pruning + T_burden + the package’s own unmodified ras_detect() / ras_validate() changepoint detection.

Tuning: At package default of 5, the local slope re-check window is located on the flat top of the signal’s plateau, so the slope looks weak even when the global changepoint test runs correctly. Increasing it to 8 increased first-pass acceptance on true signal from ~17% to ~100%. With scw fixed, the outer window_size was swept 6–20 at 100 reps per size. Small windows are with low power because the test can’t see the full causal cluster. Power increases steeply through 8 and 10, then plateaus around 0.93–0.97 from window_size = 12 on, while Type I error stays around the nominal 0.05 line the whole way (never above 0.07). The sweep metric (power \(-\) Type I) thus picks window_size = 12.

Final parameters: window_size = 12, slope_check_window_size = 8, slope.p.values.threshold.left/right = 1e-3, p.value.threshold = 1e-3 (the two thresholds taken from the package’s own @examples for a similarly-sized toy profile).

Results at these parameters (three-run validation, see Validation):

Run Type I error (95% CI) Power (95% CI) Reps
ras-ss (this work) 0.038 (0.023–0.059) 0.935 (0.891–0.965) 500 / 200
Individual-level, same data (paired) 0.036 (0.021–0.056) 0.940 (0.898–0.969) 500 / 200
Individual-level, pure ras() (num_rep = 5) 0.010 (0.003–0.023) 1.000 (0.982–1.000) 500 / 200

CIs are exact binomial 95% intervals. Rows 1 and 2 share the same simulated data per rep (same seed -> same genotypes, phenotypes, discovery weights, pruned windows, and CPD parameters), differing only in the per-window statistic (summary-stat T_burden vs. the individual-level lm(phenotype ~ LPRS) regression). Row 3 runs the package’s shipped ras() with its published default num_rep = 5, exactly as a user would.

Reading the table by rows:

  • Rows 1 vs 2 demonstrates consistency: On identical inputs the summary-stat statistic reproduced the same individual-level results on all but 1 of 200 signal reps (187 vs 188 localized to causal center) and all but 1 of 500 null reps (19 vs 18 false positives). That difference on each metric is the measured sacrifice of the summary-statistic substitution.
  • Rows 2 vs 3 shows the resampling-average benefit: The 6-point power gap (0.940 → 1.000) is between two individual-level methods, so it isn’t a summary-statistic problem. Row 2 has the individual genotypes and phenotypes considered in calculation and still falls short. The true difference is between averaging 5 random splits (num_rep = 5) versus a single pass. (NOTE: row 3 also fits its weights on a 500-person internal split rather than row 2’s 1000-person external discovery set, so this is not a perfectly isolated comparison. That difference gives row 3 noisier weights, yet row 3 still wins by 6 points, which is why the gap is attributed to the averaging, not the weight source.) The portion of this benefit that comes from re-splitting the target individuals is information that summary statistics do not contain, so no summary-statistic method can recover it. A reference-panel resampling analog (Monte Carlo resampling of the reference panel to obtain multiple burden profiles and average) can recover the LD noise portion. However, that portion is negligible in this simulation because the reference LD is estimated from 10,000 individuals, so the 6 points measured here are due to the unrecoverable target-resplit component. On real data with a smaller reference panel the analog may recover more. This is listed under Open items as the planned next experiment.
  • Type I error: ras-ss at 0.038 is within the original paper’s own Table 1 range (0.042–0.050), so it is the well-calibrated row. Pure ras() at 0.010 is actually over-conservative (below nominal 0.05 and below the paper’s range). Its lower false-positive rate also gives it perfect power, on the null side. All three rows are at or below the acceptable 0.05.
  • Peak location (the density plot in Validation): All three runs center on the true causal center (151). The single-pass runs (1 and 2) show small side-bumps at neighboring points where a noisier single-pass profile brought the validated peak one step off, a visual representation of the same averaging benefit, and the same open item. Edge-anchoring for ras-ss over the 200 power reps: 187 centered, 0 landed on the trailing plateau edge without also hitting center, 13 missed/snapped off-center.

Takeaway: ras-ss is not a replacement for the original method when individual-level data is available. In that case, the 5-averaged original remains the best option and is ~6 points more powerful, as per above simulation. ras-ss is the method that runs on the studies the original cannot: the summary-statistic-only majority of published GWAS, where there are no target individuals to resplit. On matched data with matched weights it sacrifices nothing to the individual-level method (rows 1 vs 2). The gap from the original RAS is entirely the resampling average, which is an inevitable loss in this summary-stat only environment.

Open items, not yet addressed:

  • Reference-panel LD is drawn from the same simulated population as the discovery/target cohorts here; ancestry-mismatch sensitivity (real reference panel vs. true cohort LD) has not yet been stress-tested.
  • Calibration performed on a single toy configuration (\(n=300\), \(\rho=0.8\)); has not been re-verified at other chromosome sizes or LD structures. (Done in simulation 3)
  • Summary-stat resampling analog of num_rep: The three-run validation (Validation) measures the cost of being single-pass at ~6 power points vs. the 5-averaged original, all of it due to the resampling average. A reference-panel resampling analog is the planned way to recover the LD-noise portion of that gap on real data (negligible on this toy, where reference is large). Unfortunately, the target-resplit portion is unrecoverable from summary statistics without significant change of the method.