A Worked SKAT Analysis from Genotypes to Inference

From a Sequenced Gene to a Test

A gene-level rare-variant analysis brings together several choices that are easy to blur: which variants enter the test, how their genotypes are represented, which covariates are adjusted for, and how the resulting statistic is calibrated. A SKAT p-value is meaningful only relative to those choices. This tutorial follows one hypothetical quantitative-trait analysis from its input data to its inference.

The example uses a small sample so that every step can be checked by hand. It is an arithmetic illustration, not a study design or a reliable basis for scientific conclusions. In practice, analysis also requires careful sequencing quality control and an appropriate calibration for the sample size and trait model.

For the worked analysis, assume individuals are unrelated and their errors are independent with a common variance. The quantitative trait is adequately modeled by a linear mean in the stated covariates. Genotypes are correctly called, the included variants and weights are fixed without looking at the trait values, and relevant confounding is addressed by the covariates. These assumptions matter: for example, relatedness or heteroskedasticity may require a different null model and corresponding calibration.

Define the Variants, Weights, and Kernel

A gene-level test begins by specifying its variant set. Quality and inclusion decisions should be made using genotype data and prespecified criteria, not by selecting variants because they appear associated with the trait. The dosage \(G_{ij}\) is the number of copies of the coded allele for individual \(i\) at variant \(j\). A direct kernel weight \(a_j\) is nonnegative; in the convention used here, the kernel and test statistic are

Here \(r\) is the vector of residuals from the covariate-only null model, and \(U_j\) is the score for variant \(j\). Some software describes weights through square-root weights applied to the genotype columns; that is a different parameterization of the same kernel if it produces these \(a_j\). Naming the direct kernel weights avoids ambiguity.

Weights determine how much each variant contributes to the test. Equal weights are a transparent choice for a worked calculation. In an applied analysis, weights may depend on allele frequency or annotation, but their rationale and exact parameterization should be reported. A weight of zero excludes a variant from the kernel; a larger weight gives its score more influence.

Fit the Covariate-Adjusted Null Model

For a quantitative trait with an intercept and one covariate \(x\), the null model is \(y_i=\beta_0+\beta_1x_i+\epsilon_i\), with no variant term. Write \(X\) for the columns of the intercept and covariate, and define the residual-maker matrix \(M=I-X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\). The fitted null-model residual vector is \(r=My\). The score vector \(U=G^\mathsf{T}r\) measures how each genotype column aligns with the part of the trait not explained by the covariates.

Fitting covariates is not an optional final adjustment to a score computed from the raw trait. The residuals must come from the specified null model, and the null distribution must account for the same covariate adjustment. The following hypothetical dataset has two covariate groups and includes one observation at each of the four genotype combinations in each group.

Calibrate the SKAT Statistic

Under the quantitative-trait null model, assume \(\epsilon\sim N(0,\sigma^2I)\). Then \(r=M\epsilon\) under the null and \(U=G^\mathsf{T}M\epsilon\). Consequently, \(U\) has mean zero and covariance matrix \(\sigma^2G^\mathsf{T}MG\). The quadratic statistic \(Q=U^\mathsf{T}WU\), where \(W=\operatorname{diag}(a_1,\ldots,a_m)\), has a mixture chi-square null distribution. Its weights are the nonnegative eigenvalues of \(\sigma^2W^{1/2}G^\mathsf{T}MGW^{1/2}\).

To see why, write \(C=\sigma^2G^\mathsf{T}MG\) and represent the score vector as \(U=C^{1/2}Z\), where \(Z\) has independent standard normal components on the nonzero-eigenvalue subspace. Then \(Q=Z^\mathsf{T}C^{1/2}WC^{1/2}Z\). Diagonalizing this symmetric matrix by an orthogonal transformation preserves the standard normal distribution of \(Z\), so the quadratic form becomes a sum of eigenvalues times independent squared standard normals:

Here \(\lambda_\ell\) are the eigenvalues just specified, and the \(\chi^2_{1,\ell}\) variables are independent chi-square random variables with one degree of freedom. Zero eigenvalues contribute zero. In an ordinary analysis, \(\sigma^2\) is generally unknown and must be estimated or handled by the chosen calibration procedure. The calculation below uses the residual-variance estimate as a plug-in, making the resulting p-value an approximation.

What the Result Does—and Does Not—Say

A SKAT p-value is a statement about how extreme the observed weighted score statistic is under the specified null model and its calibration. It is not an estimate of the size or direction of a single gene effect. Because the test combines squared scores, opposing variant effects do not cancel in the way they can in a simple burden sum. The result also depends on the variant set, weights, covariates, and null-model assumptions; these choices belong in the analysis report.

The calculation above illustrates the logic of a SKAT test rather than a general rule for selecting a test. If a study instead uses SKAT-O, the reported p-value must come from the calibrated hybrid procedure, including its parameter search, rather than from substituting a selected SKAT statistic into the SKAT mixture calculation. Likewise, testing many genes requires an appropriate interpretation of the gene-level p-values in light of the testing plan.

Common Mistake: Using the Wrong Null Distribution

A tempting but incorrect shortcut is to compute \(Q\) and compare it with an ordinary chi-square distribution with \(m\) degrees of freedom simply because there are \(m\) variants. That reasoning ignores both the covariate-adjusted genotype correlations and the variant weights. The SKAT reference distribution is generally a weighted sum of chi-square variables; its eigenvalues depend on the adjusted genotype matrix and the kernel weights. Only in special cases, such as the equal-eigenvalue case in the worked example, can the mixture be reduced to a scaled ordinary chi-square distribution.