Deriving the SKAT Quadratic Score Statistic

Why a Quadratic Score?

A burden test combines variant scores into one signed quantity. That is useful when variant effects are expected to point in a common direction, but signed contributions can cancel when effects differ in direction. SKAT instead retains the squared, weighted contributions of individual variants. The resulting statistic is a quadratic form in the variant-level score vector.

The derivation has two linked parts: first adjust variant scores for nuisance covariates, then ask whether their joint variation is larger than expected under no association. The random-effects formulation makes the second step precise. It leads both to the commonly reported nonnegative SKAT statistic and to a centered score for a variance component. They are related, but are not the same quantity.

The random-effects distribution is a model for how effects may vary across variants; it does not assert that variant effects are literally sampled anew in each study. The test is designed to detect departures from the null consistent with that variation. The weights in \(A\) are fixed for the test being derived. The construction and choice of variant weights are separate modeling decisions.

Efficient Variant Scores With Covariates

Begin with independent observations from a correctly specified canonical generalized linear model with unit dispersion. Conditional on covariates and genotypes, suppose \(Y_i\) has mean \(\mu_i\) and variance \(v_i\). The null model includes a nuisance-covariate matrix \(X\), with full column rank, and excludes the variants' effects. Write \(\hat{\mu}\) for the null fitted means and \(W=\operatorname{diag}(v_1,\ldots,v_n)\), evaluated at the null fit. The null nuisance parameters are assumed to be estimated regularly, with an interior fit and nonsingular nuisance information.

For a canonical link and unit dispersion, the genotype score before adjustment is \(G^\mathsf{T}(y-\hat{\mu})\). Since nuisance covariates are estimated, the relevant information for the variant scores must account for their overlap with \(X\). Define the nuisance-adjusted information projection \(P=W-WX(X^\mathsf{T}WX)^{-1}X^\mathsf{T}W\). At the null maximum-likelihood fit, \(X^\mathsf{T}(y-\hat{\mu})=0\), so the numerical efficient score can be written \(U=G^\mathsf{T}(y-\hat{\mu})\). Its efficient information matrix is \(I=G^\mathsf{T}PG\).

The projection has a direct role: genotype directions explained by the nuisance covariates do not contribute as if they were independent information about the phenotype. Under the null, the efficient score has mean zero and covariance \(I\), to the usual likelihood-based score-test approximation. For small \(\beta\), its mean changes to first order by \(I\beta\). These are the local score properties used in the derivation; they do not require fitting a separate effect for every variant.

From Random Effects to the Quadratic Form

Assume \(I\) is nonsingular on the variant-score subspace being tested. Under the local score approximation, conditional on \(\beta\), \(U\) is approximately multivariate normal with mean \(I\beta\) and covariance \(I\). Now assign the random-effects model \(\beta\sim N(0,\tau A)\), independently of the sampling variation in the score. Integrating over \(\beta\) gives a score vector with mean zero and covariance \(I+\tau IAI\). The extra term represents the score variation induced by heterogeneous variant effects.

For a mean-zero normal vector with covariance \(C(\tau)\), the log density, up to a constant, is \(-\tfrac12\log|C(\tau)|-\tfrac12U^\mathsf{T}C(\tau)^{-1}U\). Differentiate with respect to \(\tau\) at the null, where \(C(0)=I\) and \(C'(0)=IAI\). The derivative of the log determinant is \(\operatorname{tr}(C^{-1}C')\), and the derivative of the inverse is \(-C^{-1}C'C^{-1}\). Substitution therefore gives the variance-component score \(\tfrac12[U^\mathsf{T}I^{-1}(IAI)I^{-1}U-\operatorname{tr}(I^{-1}IAI)]\). Using \(I^{-1}I=I I^{-1}\) equal to the identity matrix simplifies this expression.

The conventional SKAT quadratic statistic is the uncentered quantity \(Q=U^\mathsf{T}AU\). It is nonnegative when \(A\) is positive semidefinite. The centered expression \(S_\tau(0)\) is the score for increasing the variance component from zero; it is not identical to \(Q\). Since \(\operatorname{tr}(AI)\) is fixed under the null fit, large values of \(Q\) correspond to large values of this centered score. This derivation explains why SKAT uses a weighted quadratic form rather than adding signed scores into one burden.

Worked Examples

Why Score Scaling Matters

The unit-dispersion condition is consequential. In a Gaussian identity-link model with residual variance \(\phi\), an unscaled residual score and a likelihood score differ by a factor of \(\phi\). Mixing their mean and covariance formulas would give an incorrect random-effects covariance and hence an incorrectly scaled statistic. One consistent approach is to use the likelihood score divided by \(\phi\), with its information scaled by the same likelihood convention.

For example, let \(M_X=I_n-X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\), assume \(\phi\) is known, and set \(J=G^\mathsf{T}M_XG\). The raw residual score \(U_{\rm raw}=G^\mathsf{T}M_Xy\) has conditional mean \(J\beta\) and null covariance \(\phi J\). The likelihood score \(U=U_{\rm raw}/\phi\) instead has mean \(J\beta/\phi\) and covariance \(I_{\rm eff}=J/\phi\). Under \(\beta\sim N(0,\tau A)\), its marginal covariance is \(I_{\rm eff}+\tau I_{\rm eff}AI_{\rm eff}\), so the corresponding quadratic statistic is formed from the scaled score \(U\), not \(U_{\rm raw}\).

Diagnostics, Alternatives, and Limitations

Before interpreting a SKAT score, verify that the null model is appropriate and that its fitted means and variance weights are finite and sensible. Check that the covariate design has full rank and that \(I=G^\mathsf{T}PG\) contains information in the directions being tested. A variant with no residual variation after adjustment contributes no information; a rank-deficient score system can make inverse-based derivations invalid unless the test is restricted to its informative subspace. Extremely sparse binary outcomes or near separation can also make ordinary likelihood-based score approximations unreliable.

The assumptions used above include independent observations conditional on the model, a correctly specified null mean and variance, regular nuisance-parameter estimation, and the specified dispersion convention. If observations are correlated, if the null model has an overdispersion parameter, or if a noncanonical link is used, the score covariance and projection must be derived for that model rather than copied unchanged. A mismatch between the score's scale and its information matrix is a particularly consequential implementation error.

The main conceptual alternative is a burden score, which combines variant-level scores linearly. It can be efficient when effects align in direction, but signed contributions may offset. The quadratic SKAT statistic squares and weights contributions, avoiding that cancellation in the statistic. This does not guarantee greater power in every setting: the tests encode different alternatives. A likelihood-ratio or Wald approach is another general testing route, but fitting many variant effects can be unstable when variants are rare or numerous. The score construction tests the null using the null-model fit and a set-level quadratic form.

A quadratic statistic alone does not establish the strength of evidence: its null calibration depends on the information structure and the model assumptions. Nor does rejection identify which variant is responsible, determine the sign of an individual effect, or establish a causal relationship. Those questions require additional analyses and assumptions.

• Incorrect: “The raw residual score has mean \(I\beta\) and covariance \(I\) in any generalized linear model.”
• Correction: Those statements here use the unit-dispersion score convention. With a different dispersion or score scaling, derive the mean derivative and covariance consistently before forming the quadratic statistic.
• Incorrect: “Because \(Q\) is positive, the variance-component score must be positive.”
• Correction: The score is \(\tfrac12[Q-\operatorname{tr}(AI)]\); the centering term can exceed \(Q\), so the score can be negative.