LD, Correlated Predictors, And TWAS Calibration

Why LD Matters for TWAS Calibration

A TWAS combines variant-association evidence with weights that predict a gene’s expression from genotypes. When variants are in linkage disequilibrium (LD), their association statistics are correlated. A weighted sum of those statistics therefore does not generally have the variance it would have if the variants were independent. The same dependence can link scores for different genes: two genes may have correlated predicted expression, or their weights may emphasize variants in the same LD pattern.

The earlier derivation of a TWAS test statistic gives the basic weighted score. Here the focus is its variance and dependence across genes. These calculations are central to calibration: a score divided by the wrong standard deviation can yield misleading p-values. They also matter for interpretation. Several nearby genes with small p-values do not necessarily represent several independent lines of genetic evidence.

LD-Aware Variance of One Gene’s Score

Suppose there are $m$ variants in a locus. Let $z$ be the vector of their marginal GWAS $z$-statistics, and let $R$ be the $m \times m$ LD correlation matrix for the analyzed variants. Let $w_g$ be the fixed vector of expression-prediction weights for gene $g$, in the same variant order as $z$. The unstandardized TWAS score is $S_g=w_g^\top z$.

Definition: Under the null of no genetic association with the outcome, and under the assumptions below, the marginal variant statistics are approximately $z\sim N(0,R)$. The LD-standardized score for gene $g$ is $$ Z_g=\frac{w_g^\top z}{\sqrt{w_g^\top Rw_g}}. $$ The denominator is the null standard deviation of the weighted score. If it is positive and $R$ is the correct LD matrix, then $Z_g$ is approximately standard normal under the null.

This result relies on assumptions, not just on the formula. The variant statistics must be on compatible scales and correspond to the same set of alleles and variants. Their null covariance must be represented by $R$, at least approximately. The weights are treated as fixed and independent of the GWAS association noise, and the normal approximation is adequate. In practice, population structure, relatedness, sample overlap between weight training and GWAS, or LD estimated from a poorly matched reference can complicate these assumptions. A zero denominator means the weights have no variation in the relevant LD geometry, so this standardized score is not defined.

The variance follows directly from the covariance of $z$. For any fixed vector $w_g$,

\operatorname{Var}(S_g) =\operatorname{Var}(w_g^\top z) =w_g^\top\operatorname{Cov}(z)w_g =w_g^\top Rw_g.

Thus LD contributes cross-products between weights. If two variants have positive LD and weights with the same sign, their cross-product typically increases the score variance. If the weights have opposite signs, the cross-product can reduce it. Ignoring LD is not automatically conservative: it can either overstate or understate the standard error.

Worked Example: A Mismatched LD Matrix Changes the Score

Consider two variants with hypothetical GWAS statistics $z=(2.1,1.8)^\top$, gene weights $w_g=(1,1)^\top$, and true LD correlation $0.8$. Suppose an analyst instead uses an identity matrix, which incorrectly treats the variants as uncorrelated.

With the identity matrix, the computed variance is $1^2+1^2=2$, so the score is $$ Z_{\text{identity}}=\frac{2.1+1.8}{\sqrt{2}}=\frac{3.9}{1.414}=2.758. $$

The correct LD matrix is $R=\begin{pmatrix}1&0.8\\0.8&1\end{pmatrix}$. Its variance for these weights is $$ w_g^\top Rw_g=1+0.8+0.8+1=3.6, $$ so the LD-standardized score is $$ Z_{\text{LD}}=\frac{3.9}{\sqrt{3.6}}=\frac{3.9}{1.897}=2.055. $$

For a two-sided standard normal reference, these scores give approximate p-values $0.0058$ and $0.040$, respectively. The example shows how an independence assumption can make this particular score appear more significant than it is under the stated LD model. More generally, the null variance of the incorrectly standardized score is $3.6/2=1.8$, not 1. Its standard deviation is therefore $\sqrt{1.8}=1.342$. This calculation diagnoses variance inflation; it does not assert that every LD mismatch inflates variance.

Correlated Scores Across Genes

For two genes $g$ and $h$, their scores use the same vector of variant statistics. Even if each score is correctly standardized, the scores need not be independent. Under the null model above, their covariance and correlation are

\operatorname{Cov}(Z_g,Z_h) =\operatorname{Corr}(Z_g,Z_h) =\frac{w_g^\top Rw_h} {\sqrt{(w_g^\top Rw_g)(w_h^\top Rw_h)}}.

This is a correlation of test statistics under the stated null model. It depends jointly on the LD pattern and both genes’ weight vectors. It is not determined just by genomic distance, nor by whether the genes are nearby. Strongly correlated scores can arise when weight vectors use overlapping variants or similar combinations of variants. Conversely, highly correlated variants do not force all gene scores to be highly correlated.

Worked Example: Two Genes With Strongly Correlated Scores

Suppose there are two variants with LD matrix $R=\begin{pmatrix}1&0.9\\0.9&1\end{pmatrix}$. Gene A has weights $w_A=(1,0)^\top$, and gene B has weights $w_B=(0,1)^\top$.

Each score has variance 1: $$ w_A^\top Rw_A=1,\qquad w_B^\top Rw_B=1. $$ Their covariance is $$ w_A^\top Rw_B=(1,0) \begin{pmatrix}1&0.9\\0.9&1\end{pmatrix} \begin{pmatrix}0\\1\end{pmatrix}=0.9. $$ Thus $\operatorname{Corr}(Z_A,Z_B)=0.9$. If the hypothetical variant statistics are $z=(2.3,2.3)^\top$, then both gene scores equal 2.3 and each has a two-sided standard normal p-value of about $0.0214$.

Those two p-values are not independent confirmations of the same result: the scores have null correlation 0.9 because the variants are in strong LD. The calculation establishes dependence between the statistics, not that either gene causes the outcome.

Worked Example: Opposite Weights Can Cancel LD Dependence

Keep the same two-variant setting, but let the LD correlation be $0.8$. Take $w_A=(1,1)^\top$ and $w_B=(1,-1)^\top$. The score variances are $$ w_A^\top Rw_A=1+0.8+0.8+1=3.6, $$ and $$ w_B^\top Rw_B=1-0.8-0.8+1=0.4. $$ The covariance is $$ w_A^\top Rw_B=(1,1) \begin{pmatrix}1&0.8\\0.8&1\end{pmatrix} \begin{pmatrix}1\\-1\end{pmatrix}=0. $$ Therefore, the standardized scores are uncorrelated under the null, despite the strong LD between variants. Under the joint normal approximation they are independent. For the hypothetical statistics $z=(2,1)^\top$, the scores are $3/\sqrt{3.6}=1.581$ and $1/\sqrt{0.4}=1.581$. This illustrates why gene-score dependence must be calculated from the weights and LD together, rather than inferred from LD alone.

Using the Wrong LD: A Calibration Diagnostic

Suppose the true covariance of the variant statistics is $R_{\text{true}}$, but the score denominator uses an estimated or reference matrix $R_{\text{ref}}$. For fixed weights, the statistic $\widetilde{Z}_g=w_g^\top z/\sqrt{w_g^\top R_{\text{ref}}w_g}$ has null variance

\operatorname{Var}(\widetilde{Z}_g) =\frac{w_g^\top R_{\text{true}}w_g} {w_g^\top R_{\text{ref}}w_g}.

The ratio equals 1 only when the two matrices give the same variance in the direction of the gene’s weights. They need not be identical in every entry for that particular ratio to equal 1; equally, a small overall difference between matrices can matter if it occurs in a direction heavily weighted by the gene model. A ratio above 1 means the statistic is more variable than a standard normal reference assumes, tending toward anti-conservative p-values. A ratio below 1 means it is less variable, tending toward conservative p-values.

A short derivation makes the distinction clear. Under the null, $\operatorname{Cov}(z)=R_{\text{true}}$, so

\operatorname{Var}(w_g^\top z) =w_g^\top R_{\text{true}}w_g.

Dividing a random variable by the fixed denominator $\sqrt{w_g^\top R_{\text{ref}}w_g}$ divides its variance by $w_g^\top R_{\text{ref}}w_g$, giving the ratio above. It is therefore the LD geometry relevant to each weight vector—not a generic statement that one reference is “close enough”—that determines this part of calibration.

Diagnosing Problems and Choosing an Approach

Several checks help distinguish a real pattern of association statistics from a problem with their variance model:

Check variant alignment. Confirm that alleles, variant order, and the variants included in the weights and LD matrix agree. A swapped allele or mismatched variant can change a weight’s sign or attach it to the wrong LD entries.
Check the LD source and matrix. The LD sample should be relevant to the GWAS population, and the matrix should be symmetric and positive semidefinite up to numerical tolerance. Severe discrepancies, or an unstable matrix from a small reference sample, can make the calculated variance unreliable.
Inspect weight-specific variance. Compare $w_g^\top Rw_g$ across plausible LD estimates. Large changes flag genes whose standardization is sensitive to the reference. This is a sensitivity check, not proof that one estimate is correct.
Look at patterns across loci. Excessive dispersion of null-like scores can indicate calibration problems, but correlated gene tests mean that counting small p-values as if they were independent is not a valid diagnostic by itself. Simulations or resampling that preserve LD and the analysis design can provide a more relevant check.

A marginal, single-gene TWAS score is convenient and answers whether that gene’s weighted predictor is associated with the outcome, with its variance accounting for LD. A joint regression of several predicted-expression predictors asks a different question: the association of each predictor conditional on the others. If those predictors are strongly correlated, the joint coefficients can be imprecise and sensitive to small changes in the data or model. LD-aware simulation or resampling can assess calibration while retaining dependence, but results depend on whether the simulated data and null model represent the study adequately. These approaches are complementary rather than interchangeable.

A practical interpretation limit follows: several nearby genes with significant marginal TWAS scores may be responding to the same underlying variant-association pattern. The scores can be correlated even when the genes’ distinct biological roles are unknown. A cluster of p-values alone does not identify which gene is responsible, establish that altered expression mediates an effect, or provide multiple independent confirmations. LD-aware standardization addresses a variance question; it does not resolve all questions about biological interpretation.

Common pitfall: “The two nearby genes both have small TWAS p-values, so they are independent evidence for two affected genes.” This reasoning ignores the covariance $w_g^\top Rw_h$. Calculate or estimate the score correlation from both weight vectors and LD before treating the results as distinct evidence. Even low score correlation would not, on its own, establish a causal interpretation.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

For $w=(1,1)^\top$ and two-variant LD correlation $0.5$, calculate $w^\top Rw$. Why is it not the variance obtained by treating variants as independent?
If $R_{\text{true}}$ gives a score variance of 5 and $R_{\text{ref}}$ gives a denominator variance of 4, what is the null variance of the reference-standardized score, and which direction of calibration error is suggested?
Two correctly standardized gene scores have null covariance 0.6. What is their correlation, and what does this imply about interpreting two small p-values?
With $R=\begin{pmatrix}1&0.7\\0.7&1\end{pmatrix}$, find the covariance between scores with weights $(1,1)^\top$ and $(1,-1)^\top$.

Solutions. (1) The variance is $1+0.5+0.5+1=3$, whereas independence gives $1+1=2$; the LD cross-products contribute $2(0.5)=1$. (2) The null variance is $5/4=1.25$, so the statistic is more variable than a standard normal reference and p-values may be anti-conservative. (3) Since each score has variance 1, their correlation is 0.6. Their p-values are dependent and should not be treated as independent confirmations. (4) The covariance is $(1,1)R(1,-1)^\top=0$: the first component of $R(1,-1)^\top$ is $0.3$, the second is $-0.3$, and their sum is zero.