SKAT for Quantitative Traits

Why Quantitative-Trait SKAT Uses a Linear Null Model

A gene can contain several rare variants whose effects on a quantitative phenotype differ in direction and magnitude. SKAT tests their joint contribution without requiring all variant effects to point in the same direction. For a quantitative outcome, the test starts with a linear model for the phenotype under the null hypothesis of no association between the gene’s variants and the outcome. Covariates such as age, sex, or ancestry variables are included in this null model so that the gene-level test targets association beyond those modeled covariate effects.

The test is a variance-component test: it asks whether the pattern of variant effects is more variable than expected when all effects are zero. Its statistic is a weighted quadratic form in the null-model residuals. Consequently, its null distribution is generally a mixture of chi-square distributions, not a single chi-square distribution. The covariate adjustment, variant weights, and mixture calibration must be kept together when calculating a p-value.

Quantitative-Trait Null Model and SKAT Statistic

The design matrix must contain the intended covariates, and the linear mean model and error assumptions must be adequate for the chosen calibration. In particular, unexplained relatedness or heteroscedasticity can invalidate the independent, constant-variance error model; such settings require an appropriate covariance model and corresponding test calibration. The genotype matrix and weights are treated as fixed for the null-distribution calculation. Usually the weights are selected before testing, based on the analysis plan and variant information.

Define the residual projection matrix \(P=I-X(X'X)^{-1}X'\) when \(X\) has full column rank. More generally, use the orthogonal projection onto the residual space, \(P=I-XX^+\), where \(X^+\) is the Moore–Penrose inverse. Then \(r=Py\), and the residual degrees of freedom are \(n-\operatorname{rank}(X)\). If estimating residual variance from the null model, the usual estimate is \(s^2=r'r/[n-\operatorname{rank}(X)]\), provided the residual degrees of freedom are positive. This rank-based statement also covers a rank-deficient design matrix.

This tutorial uses \(Q=r'GWG'r\), without a factor of one-half. Other presentations may include a constant factor in the statistic; that is valid only when the null-distribution eigenvalues are scaled by the same factor. Under the Gaussian null model, the relevant eigenvalues are those of \(\sigma^2 W^{1/2}G'PGW^{1/2}\). If \(\sigma^2\) is treated as known, \(Q\) has the mixture distribution given by those eigenvalues. When \(\sigma^2\) is unknown, it is commonly estimated under the null and used in the mixture calibration; the resulting calibration is a plug-in procedure, and software may use finite-sample refinements.

Fitting the Null Model

Calculating the Gene Statistic and Its Null Distribution

The covariate adjustment removes the part of each genotype explained by the null-model covariates. Algebraically, \(G'PG\) is the matrix of cross-products among the residualized genotype columns. With weights, the matrix \(W^{1/2}G'PGW^{1/2}\) determines the mixture’s eigenvalues. If these eigenvalues are \(\lambda_1,\ldots,\lambda_m\), then under the null, \(Q\) has the distribution of \(\sigma^2\sum_{k=1}^{m}\lambda_k\chi^2_{1,k}\), where the \(\chi^2_{1,k}\) are independent chi-square random variables with one degree of freedom. Zero eigenvalues contribute nothing. In practice, the upper-tail probability is calculated using a mixture-distribution method, rather than by comparing \(Q\) with an ordinary chi-square distribution having \(m\) degrees of freedom.

Why the Mixture Appears

A short derivation shows why the eigenvalues, rather than simply the number of variants, determine the null distribution. Under \(H_0\), \(y=X\beta+\varepsilon\), and \(P X=0\). Therefore \(r=Py=P\varepsilon\). Write \(K=GWG'\). The statistic is \(Q=\varepsilon'PKP\varepsilon\). Since \(\varepsilon=\sigma z\) for \(z\sim N(0,I)\), it follows that \(Q=\sigma^2z'PKPz\).

The symmetric matrix \(PKP\) is positive semidefinite. Diagonalizing it gives nonnegative eigenvalues; an orthogonal transformation of \(z\) still has independent standard normal components. Each squared component is therefore an independent \(\chi^2_1\) random variable, weighted by its eigenvalue. The nonzero eigenvalues of \(PKP\) agree with those of \(W^{1/2}G'PGW^{1/2}\), so the two forms yield the same mixture. This justifies both the covariate adjustment and the eigenvalue-based calibration.

Interpretation and Common Pitfalls

A small SKAT p-value is evidence against the null that the tested variants make no contribution to the quantitative phenotype, conditional on the modeled covariates and the assumptions used for calibration. It is a gene-level result. By itself, it does not establish that every variant is associated, identify a causal variant, or show that genetic variation causes a change in phenotype. Nor does the statistic give a single per-allele effect estimate: its purpose is to test a weighted collection of effects through their joint contribution.

• Using unadjusted genotype cross-products. The calibration must use \(G'PG\), not \(G'G\), when covariates are included in the null model.
• Using the wrong residual degrees of freedom. The correct count is \(n-\operatorname{rank}(X)\); it equals \(n-p\) only when \(X\) has full column rank.
• Comparing with an ordinary chi-square distribution. The eigenvalues generally differ, so use the weighted mixture distribution for the p-value.
• Rescaling only one part of the test. If all weights are rescaled, the statistic and all mixture eigenvalues must be rescaled together.
• Overstating the conclusion. A significant result supports gene-level association under the model; it does not prove causality or specify which variant is responsible.