Adjusting SKAT for Covariates

Why Covariate Adjustment Is Part of the Test

A gene-level test asks whether a set of genetic variants is associated with a trait, while accounting for other measured variables that may also relate to the trait. Age, sex, ancestry variables, and study indicators are common examples of covariates. In SKAT, adjustment is not simply a preliminary step that changes the trait values: it also changes the covariance of the genetic score and therefore the null distribution used to calculate a p-value.

The central quantities are the residuals from a model with no genetic effects, the genetic score calculated from those residuals, and the covariance of that score after estimating the covariate effects. For a quantitative trait, the residual-maker matrix removes the linear contribution of the covariates. For a binary trait, the corresponding adjustment uses the fitted Bernoulli variances and the information matrix for the covariates.

This tutorial treats \(W\) as fixed and focuses on how covariate adjustment determines \(U\) and its null covariance. The weights are not selected by looking at the observed association statistic.

Quantitative Traits: Residualize Both the Trait and the Genotypes

For a quantitative trait, use the linear model \(y=X\alpha+G\beta+\varepsilon\), where \(y\) is the \(n\)-vector of trait values, \(\alpha\) is the \(p\)-vector of covariate coefficients, and \(\beta\) is the \(m\)-vector of genetic effects. Under the null hypothesis \(H_0:\beta=0\), assume \(\varepsilon\sim N(0,\sigma^2 I_n)\), with independent errors, constant variance, and a full-column-rank \(X\). The covariate effects are estimated under this null model.

Define the ordinary least-squares projection onto the columns of \(X\) and the residual-maker matrix as follows. The latter is symmetric and removes any component lying in the covariate space.

The null-model residual vector is \(e=M_Xy\). Since \(M_XX=0\), these residuals contain no fitted linear component from the columns of \(X\). The SKAT residual-score vector is \(U=G'e\). Equivalently, \(U=(M_XG)'y\): the same score results from first residualizing each genotype column against the covariates and then relating those residualized genotypes to the trait.

Here \(U\) is the unscaled residual-score vector. In a Gaussian likelihood, the score for \(\beta\) is \(U/\sigma^2\); the unscaled convention makes the variance factor explicit in the covariance and in the null distribution below.

The factor \(\sigma^2\) is essential. The weights are eigenvalues of the covariance of the weighted score, not just of the genotype information matrix. Omitting this factor would give the wrong scale unless \(\sigma^2=1\).

Why the Quantitative-Trait Null Is a Mixture

The mixture representation follows directly from normality and a symmetric-matrix decomposition. Under \(H_0\), \(U=G'M_X\varepsilon\) is multivariate normal with covariance \(\Sigma_U=\sigma^2G'M_XG\). Let \(A=W^{1/2}\Sigma_UW^{1/2}\), and write \(Q=(W^{1/2}U)'(W^{1/2}U)\). Diagonalize the symmetric positive-semidefinite matrix \(A\) as \(A=V\Lambda V'\), where \(V\) is orthogonal and \(\Lambda\) has eigenvalues \(\lambda_1,\ldots,\lambda_m\), including any zeros.

An orthogonal transformation preserves the squared length of a vector. The transformed, standardized normal coordinates are independent standard normal variables \(Z_1,\ldots,Z_m\). Therefore the quadratic form has the stated mixture representation.

This argument also explains what happens if some variants become redundant after adjustment: zero eigenvalues contribute zero to the sum. If \(\sigma^2\) is estimated rather than known, a common calibration substitutes the null-model estimate into the eigenvalues. That plug-in procedure is not the same exact known-variance distribution; its accuracy depends on the sample size and model conditions.

Binary Traits: Use the Logistic Null Model's Information Projection

For a binary trait, write \(Y_i\in\{0,1\}\) and fit a logistic null model with linear predictor \(X_i'\alpha\) and fitted probability \(\hat{\mu}_i\). The model assumes independent Bernoulli outcomes conditional on the covariates, with \(\operatorname{logit}(\mu_i)=X_i'\alpha\). Assume the null fit is finite and \(X'VX\) is invertible, where \(V=\operatorname{diag}\{\hat{\mu}_i(1-\hat{\mu}_i)\}\). The invertibility condition ensures that the covariate coefficients have estimable information.

Let \(e=Y-\hat{\mu}\) be the fitted null residual vector. At the null-model maximum-likelihood estimate, \(X'e=0\), so the observed genetic score is \(U=G'e\). However, its covariance for testing genetic effects must account for estimating \(\alpha\). The relevant information projection is

This matrix adjusts the Bernoulli variance information for the covariate parameters. It is generally not an ordinary idempotent residual-maker like \(M_X\). The covariance of the efficient genetic score is \(\operatorname{Cov}(U)\approx G'P_XG\); the approximation refers to the usual score-test, large-sample normal approximation.

Several Variants: The Full Weighted Covariance Matters

With multiple variants, the score components need not be independent. Correlation can arise because variants are correlated in the sample, because their genotype columns overlap after covariate adjustment, or both. The correct calibration uses the eigenvalues of the entire weighted score covariance matrix. Using the eigenvalues of the unweighted covariance alone discards the effect of \(W\); using the genotype information alone for a quantitative trait also discards \(\sigma^2\).

Practical Checks and a Common Mistake

The adjustment matrices encode which parts of the data remain informative about genetic effects after the covariates have been fitted. If a genotype column is fully explained by the covariates, its adjusted quantitative information is zero. For a binary trait, the analogous loss of information is reflected in \(G'P_XG\). Such a variant cannot contribute independent evidence beyond the covariate model.

A common error is to calculate \(Q\) from an adjusted score but use unadjusted eigenvalues, or to use the eigenvalues of \(\operatorname{Cov}(U)\) without incorporating \(W\). For quantitative traits, another frequent error is to omit the residual variance scale. The correction is to form the covariance appropriate to the fitted null model first, then form the weighted covariance and obtain its eigenvalues.

The essential distinction is between the observed score and its null reference distribution. Covariates affect both. A p-value computed from an unadjusted covariance does not become valid simply because the observed residuals were adjusted.