Linkage Disequilibrium Effects on SKAT

Why linkage disequilibrium changes SKAT calibration

SKAT combines variant-level score contributions through a weighted quadratic statistic. When variants are correlated through linkage disequilibrium (LD), their scores are correlated under the null as well. This changes the eigenvalues that determine the statistic’s null distribution. The observed quadratic statistic is still calculated from the scores and weights, but its reference distribution—and therefore its p-value—must reflect their covariance.

The earlier material on kernel construction and variant weights defines the quadratic form. Here the focus is the effect of LD on its spectrum and calibration. A key distinction is that genotype correlation and score correlation are related, but not always identical: covariate adjustment and the null model affect the covariance of the scores used in the test.

Definition: Let \(U\) be the \(m\)-vector of variant score statistics under the null, with covariance matrix \(C\), and let \(W\) be a diagonal matrix with nonnegative SKAT weights. The SKAT quadratic statistic is \(Q=U^{T}WU\). Under a multivariate normal null approximation, its distribution is a weighted sum of independent chi-square random variables, with weights equal to the eigenvalues of \(W^{1/2}CW^{1/2}\). Correlation among variant scores appears in the off-diagonal entries of \(C\).

Assumptions and the eigenvalue representation

The eigenvalue representation requires that, under the null hypothesis of no association for the tested set, the score vector is approximately multivariate normal with mean zero and covariance \(C\). The null model must be appropriate for the phenotype and design, and the covariance must account for the nuisance parameters and covariates used in constructing the scores. The weights in \(W\) are fixed for this calibration, and \(C\) must be positive semidefinite. The approximation can be less reliable with extremely sparse data or poorly estimated null-model quantities.

If \(C\) is a covariance matrix, it can be factored as \(C^{1/2}C^{1/2}\). Write \(U=C^{1/2}Z\), where \(Z\) has independent standard normal entries under the approximation. Then \(Q=Z^{T}C^{1/2}WC^{1/2}Z\). The symmetric matrix \(C^{1/2}WC^{1/2}\) has the same nonzero eigenvalues as \(W^{1/2}CW^{1/2}\). Diagonalizing the quadratic form therefore gives the null representation below.

Q \ \overset{d}{\approx}\ \sum_{j=1}^{m}\lambda_j \chi^2_{1,j}, \qquad \lambda_j=\text{eigenvalues of }W^{1/2}CW^{1/2},

Here the \(\chi^2_{1,j}\) variables are mutually independent chi-square variables with one degree of freedom. Zero eigenvalues contribute nothing. If the score components are independent, \(C\) is diagonal and the eigenvalues reduce to the weighted score variances. With LD, off-diagonal entries change the eigenvalues, so treating the variants as independent generally gives the wrong null distribution.

How correlation changes the spectrum

For clarity, suppose two score components have unit null variances, correlation \(\rho\), and weights \(w_1,w_2\). Then \(C\) is a correlation matrix and the weighted kernel is as follows. Its trace is the sum of its eigenvalues, and its determinant is their product.

C= \begin{pmatrix} 1 & \rho\\ \rho & 1 \end{pmatrix}, \qquad W^{1/2}CW^{1/2}= \begin{pmatrix} w_1 & \rho\sqrt{w_1w_2}\\ \rho\sqrt{w_1w_2} & w_2 \end{pmatrix}.

The diagonal weights are unchanged by correlation, but the off-diagonal terms affect the eigenvalues. In general, the trace, and hence the null mean, is unchanged when the score variances and weights remain fixed. The spread of the null distribution can nevertheless change substantially. For more than two variants, the full correlation structure—not just one summary of pairwise LD—determines the spectrum.

Worked Example: Two Variants with Unequal Weights

Suppose two standardized null scores have correlation \(0.6\), and their quadratic-form weights are \(1\) and \(4\). Compare the weighted kernel with the one obtained by incorrectly assuming independence.

Correlated kernel. Substitution gives \(K=\begin{pmatrix}1&0.6\sqrt{1\cdot4}\\0.6\sqrt{1\cdot4}&4\end{pmatrix} =\begin{pmatrix}1&1.2\\1.2&4\end{pmatrix}\). Its trace is \(5\), and its determinant is \(4-1.2^2=2.56\). Thus the eigenvalues solve \(\lambda^2-5\lambda+2.56=0\):

\lambda_{1,2} =\frac{5\pm\sqrt{25-4(2.56)}}{2} =\frac{5\pm\sqrt{14.76}}{2} \approx 4.4209,\ 0.5791.

Under the normal approximation, \(Q\) has null distribution \(4.4209\chi^2_{1,1}+0.5791\chi^2_{1,2}\). Its mean is \(4.4209+0.5791=5\), and its variance is \(2(4.4209^2+0.5791^2)\), approximately \(39.76\).

Independence approximation. Setting the correlation to zero makes the kernel diagonal with eigenvalues \(1\) and \(4\). The null mean remains \(1+4=5\), but the variance is \(2(1^2+4^2)=34\). Thus ignoring this correlation understates the null variance. The test statistic for a particular observed score vector is still calculated using its scores and weights; it is the reference distribution used to calibrate that statistic that changes.

For any eigenvalues \(\lambda_j\), the first two null moments follow from the mean \(1\) and variance \(2\) of a \(\chi^2_1\) variable. Independence of the chi-square terms means that their variances add.

E(Q)=\sum_j\lambda_j, \qquad \operatorname{Var}(Q)=2\sum_j\lambda_j^2.

The variance formula shows why LD matters even when the mean does not change: it changes the eigenvalues and hence the sum of their squares. A common approximation matches these two moments to a scaled chi-square distribution. If \(Q\) is approximated by \(c\chi^2_{\nu}\), equating means gives \(c\nu=\sum_j\lambda_j\); equating variances gives \(2c^2\nu=2\sum_j\lambda_j^2\). Solving yields the following calibration quantities.

c=\frac{\sum_j\lambda_j^2}{\sum_j\lambda_j}, \qquad \nu=\frac{(\sum_j\lambda_j)^2}{\sum_j\lambda_j^2}.

This Satterthwaite approximation is a convenient moment-matching method, not an exact identity for a general mixture. Methods that evaluate the weighted chi-square tail directly, such as numerical algorithms for quadratic forms, use the eigenvalue mixture more fully. Whichever method is used, it must be based on the appropriate score covariance.

Worked Example: LD Among Three Equally Weighted Scores

Suppose three standardized scores have the same pairwise correlation, \(0.5\), and all three weights equal \(1\). This is a valid correlation matrix: its eigenvalues are \(2,0.5,0.5\), all nonnegative.

Because the weights are all \(1\), the weighted kernel is the correlation matrix itself. For an equicorrelation matrix of size three, the eigenvalue in the direction \((1,1,1)\) is \(1+2\rho\), while the two orthogonal directions have eigenvalue \(1-\rho\). Substituting \(\rho=0.5\) gives \(2,0.5,0.5\).

The null mean is \(2+0.5+0.5=3\). The variance is \(2(2^2+0.5^2+0.5^2)=2(4.5)=9\). Under independence, the three eigenvalues would all be \(1\), giving the same mean \(3\) but variance \(2(1+1+1)=6\). Positive LD has concentrated the spectrum into one larger eigenvalue and two smaller ones, increasing the null variance in this example.

For a moment-matched scaled chi-square, \(c=(4+0.25+0.25)/3=1.5\) and \(\nu=3^2/4.5=2\). Thus the approximation is \(1.5\chi^2_2\). Under independence, the same calculation gives \(c=1\) and \(\nu=3\), or \(\chi^2_3\). These are different reference distributions despite having the same mean.

Calculating and checking a calibration

In practice, one obtains the score covariance under the fitted null model, forms the weighted kernel, computes its eigenvalues, and then evaluates the upper-tail probability for the observed \(Q\). A calibration based on the eigenvalues reflects both the weights and the dependence among scores. A diagonal approximation discards the off-diagonal covariance and is justified only when those correlations are negligible for the tested set.

Obtain the null score covariance.
Use the covariance appropriate to the null model and covariate adjustment. Do not automatically substitute raw genotype correlations if the score covariance after adjustment is available.

Form the weighted kernel.
Calculate \(W^{1/2}CW^{1/2}\), keeping the same weights used in the observed quadratic statistic.

Inspect the spectrum.
Check that eigenvalues are nonnegative up to numerical tolerance, that their sum agrees with the trace, and that their squares produce the variance formula.

Evaluate the upper tail.
Use the eigenvalue mixture or a stated approximation to calculate the p-value. The null distribution must match the covariance used for the score vector.

Worked Example: Moment Calibration with the Two-Variant Spectrum

Continue with the two-variant correlated kernel, whose eigenvalues are approximately \(4.4209\) and \(0.5791\). Calculate the scaled chi-square moment approximation and compare it with the independence approximation.

For the correlated spectrum, the eigenvalue sum is \(5\). The sum of squared eigenvalues is \(4.4209^2+0.5791^2\), approximately \(19.88\). Therefore \(c=19.88/5=3.976\), and \(\nu=5^2/19.88\approx1.258\). This gives the moment-matched reference \(3.976\chi^2_{1.258}\). Its mean is approximately \(3.976(1.258)=5.00\), and its variance is approximately \(2(3.976^2)(1.258)=39.76\), agreeing with the mixture moments up to rounding.

Under independence the eigenvalues are \(1\) and \(4\), so their sum is \(5\) and their squared sum is \(17\). Hence \(c=17/5=3.4\) and \(\nu=25/17\approx1.471\), giving \(3.4\chi^2_{1.471}\). Its mean is \(5\) and variance is \(2(3.4^2)(1.471)\approx34\). These approximations have different shapes and variances, so using the independent-variant calibration for correlated scores can change the reported p-value. The direction and size of that change depend on the observed \(Q\); the variance comparison alone does not determine a particular p-value.

Diagnostics, alternatives, and limitations

Several checks can reveal a faulty LD-aware calibration. A materially negative eigenvalue suggests a non-positive-semidefinite covariance estimate or a numerical problem; tiny negative values near computational tolerance may instead be rounding artifacts. The eigenvalue sum should equal the trace of the weighted kernel, and \(2\sum_j\lambda_j^2\) should agree with the variance calculated from that kernel. If changing the assumed correlation leaves a purportedly unchanged spectrum without justification, verify that the off-diagonal terms were included.

LD itself should be assessed in the analyzed sample and variant set. Strong correlation can make the kernel effectively low-rank: one or a few eigenvalues may account for much of the total, even when many variants were included. This does not mean that the number of variants has changed; it means the quadratic statistic has fewer dominant independent directions under its null approximation.

The main contrast is with a burden test, which first combines variant scores into a single weighted score and then tests that aggregate. SKAT instead uses a quadratic form that accumulates weighted squared score contributions. Their treatment of dependence differs: a burden statistic’s variance depends on covariances among its terms, while SKAT’s null mixture is determined by the eigenvalues of its weighted score covariance. Neither method should be calibrated by pretending correlated scores are independent unless that approximation is warranted. Direct mixture-tail evaluation is an alternative to moment matching; it can better represent the mixture but still relies on the score covariance and null approximation being appropriate.

The normal-score mixture can be inaccurate when score distributions are highly discrete or sparse, or when null-model estimation is unstable. Estimating LD imprecisely can also distort small eigenvalues or the upper tail. The eigenvalue calculation accounts for dependence as represented in \(C\); it does not repair a misspecified phenotype model, an unsuitable variant set, or errors in the estimated covariance. Calibration should therefore be checked in the context of the null model and data structure, not treated as an automatic consequence of diagonalizing a matrix.

Common mistake: “The weights and the diagonal of the score covariance are unchanged, so the SKAT null distribution is unchanged.” This is incorrect: off-diagonal covariance changes the kernel eigenvalues, and the mixture variance is \(2\sum_j\lambda_j^2\). Use the covariance-aware spectrum for calibration, and verify its eigenvalue sum and squared sum against the kernel’s trace and null variance.

Check Your Understanding

Try each question, then compare with the solutions below.

For two unit-variance scores with correlation \(0.3\) and equal weights \(1\), find the kernel eigenvalues, null mean, and null variance.
Why can the null mean remain unchanged when LD changes, while the null variance changes?
A weighted kernel has eigenvalues \(3\) and \(1\). Calculate the mean, variance, and moment-matched \(c\) and \(\nu\).
What does a materially negative eigenvalue of an estimated weighted kernel suggest, and what should be checked?

Solutions. (1) The kernel is \(\begin{pmatrix}1&0.3\\0.3&1\end{pmatrix}\), whose eigenvalues are \(1.3\) and \(0.7\). The mean is \(2\), and the variance is \(2(1.3^2+0.7^2)=2(2.18)=4.36\). (2) With fixed weights and unit score variances, the trace is fixed, so the eigenvalue sum and \(E(Q)\) are fixed; correlation can redistribute that sum, changing the sum of squared eigenvalues and hence the variance. (3) The mean is \(3+1=4\), the variance is \(2(9+1)=20\), \(c=10/4=2.5\), and \(\nu=16/10=1.6\). (4) It can indicate a covariance estimate that is not positive semidefinite or a numerical error. Check how the covariance was estimated and formed, and distinguish a material negative value from a tiny rounding artifact.