Deriving the SKAT Score Test

Why Derive the Score?

A gene-level score test begins with a model for how several variant effects contribute to a phenotype. In a burden model, variant effects are combined into a single direction; a variance-component model instead allows their effects to vary around zero. The score for that variance component reveals exactly how the genotype data, variant weights, and departures from the null phenotype model combine in the test.

This derivation uses a Gaussian phenotype model with a linear mean and a variance component for variant effects. It establishes the score and its residual and trace terms. It does not require the variance component to be estimated under the alternative: a score test evaluates the derivative of the likelihood at the null value, where the variant-effect variance is zero.

The columns of \(G\) contain the observed genotype counts, conventionally coded as 0, 1, or 2 copies of an allele, subject to the chosen analysis coding. The weight \(w_j\) sets the relative variance assigned to variant \(j\)'s effect. Thus, if \(w_j=0\), that variant contributes no random-effect variance in this model. The matrix \(W\) describes effect variance, not observed genotype variation.

From Variant Effects to a Phenotype Covariance

Marginalizing over the random vector \(b\) gives a model for the phenotype alone. Since \(E(b)=0\), the marginal mean is \(X\gamma\). Independence of \(b\) and \(\varepsilon\), together with their covariance matrices, gives

The \(n\times n\) matrix \(K\) is the covariance pattern contributed by the variants before multiplication by the variance-component parameter \(\tau\). Its entry for individuals \(i\) and \(i'\) is \(K_{ii'}=\sum_{j=1}^m w_jG_{ij}G_{i'j}\). This identity follows by multiplying the matrices in \(GWG^\mathsf{T}\); it also shows why both the genotypes and the weights determine the covariance.

Under \(H_0\), \(\tau=0\), so the variants do not add a covariance component and the phenotype model reduces to the null model with mean \(X\gamma\) and residual variance \(\sigma^2\). Define the null-model residuals as \(r=y-X\hat\gamma_0\), where \(\hat\gamma_0\) is the null-model estimate of the covariate coefficients. The residuals are not generally the raw phenotypes: covariate-associated mean structure must first be accounted for.

Deriving the Likelihood Score

Write \(V(\tau)=\sigma^2I+\tau K\). If the mean \(X\gamma\) were known, the Gaussian log likelihood, apart from a constant, would be \(-\tfrac12\{\log|V|+(y-X\gamma)^\mathsf{T}V^{-1}(y-X\gamma)\}\). Differentiating with respect to \(\tau\) uses \(\frac{d}{d\tau}\log|V|=\operatorname{tr}(V^{-1}K)\) and \(\frac{d}{d\tau}V^{-1}=-V^{-1}KV^{-1}\). Therefore, at any \(\tau\), the derivative is

The first term rewards alignment between the phenotype departures and the covariance pattern \(K\); the trace term accounts for the increase in covariance that the model would impose even when the data do not show such alignment. A score test evaluates this derivative under the null model, using null estimates of nuisance parameters.

When covariate coefficients are estimated, the exact trace term depends on the likelihood formulation. For ordinary likelihood with the mean profiled out, define \(P=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\), assuming \(X\) has full column rank. The profiled likelihood score is

In particular, the quadratic term uses \(P\), but the trace term remains \(\operatorname{tr}(V^{-1}K)\). Restricted maximum likelihood (REML) uses a different likelihood, based on contrasts that remove the fixed-effect mean. Its score is

These are distinct score formulas. The REML trace term is not obtained by merely replacing \(V^{-1}\) by \(P\) in the ordinary profiled-likelihood score; it comes from differentiating the restricted likelihood. Keeping the distinction explicit prevents a seemingly small algebraic change from combining terms from two different likelihoods.

The Null Score and the SKAT Quadratic Component

At \(\tau=0\), \(V=\sigma^2I\), and \(P=\sigma^{-2}M\), where \(M=I-X(X^\mathsf{T}X)^{-1}X^\mathsf{T}\) is the residual-maker matrix. It projects a vector onto the space orthogonal to the columns of \(X\). Under REML, let \(r=My\). Substitution into the REML score gives

The first term, \(r^\mathsf{T}Kr\), is the central phenotype-genotype quadratic component. To see its variant-level form, substitute \(K=GWG^\mathsf{T}\):

where \(G_j\) is column \(j\) of \(G\). Thus, each variant contributes its weight times the square of its null-model residual score. Unlike a signed sum of variant scores, this quadratic component does not cancel simply because some variant scores are positive and others are negative.

The full likelihood score also contains a centering term. Under REML, that term is \(\sigma^2\operatorname{tr}(MK)\) in the displayed null expression; it is the null expectation of \(r^\mathsf{T}Kr\) when the Gaussian residual covariance is \(\sigma^2M\). In many presentations, the SKAT quadratic statistic is written \(Q=r^\mathsf{T}Kr\), while its null calibration accounts for the null distribution of that quadratic form. The score derivation explains where \(Q\) comes from; the centered likelihood score and the reported quadratic statistic should not be silently treated as identical quantities.

Assumptions, Interpretation, and a Common Pitfall

The derivation assumes the specified Gaussian model, independent residual errors with covariance \(\sigma^2I\), a correctly specified null mean \(X\gamma\), and a random variant-effect vector with mean zero and covariance \(\tau W\). The weights are fixed for the test, and \(W\) must be positive semidefinite for it to define a valid covariance. The genotype matrix and phenotype residuals must refer to the same individuals in the same row order. If residual covariance is not \(\sigma^2I\), the corresponding null covariance must be modeled and used in \(V\) and \(P\), rather than applying the identity-covariance expressions without adjustment.

A score is a local likelihood derivative, not an estimated effect size and not, by itself, a p-value. A positive value means the chosen likelihood initially rises in the positive-\(\tau\) direction at the null; it does not establish causation, identify an associated variant, or quantify a phenotype change per allele. The test evaluates a gene-level variance component, and its calibration must reflect the null distribution of the statistic.

Common pitfall: A frequent error is to write the profile-likelihood quadratic term with \(P\) and then replace the trace term by \(\operatorname{tr}(PK)\), as though both changes followed from profiling out the mean. They do not. For ordinary likelihood with the mean profiled out, the score is \(\tfrac12\{y^\mathsf{T}PKPy-\operatorname{tr}(V^{-1}K)\}\). The score with \(\operatorname{tr}(PK)\) is the REML score, derived from the restricted likelihood. State which likelihood is being differentiated before using its formula.

A second related error is to call \(Q=r^\mathsf{T}Kr\) the entire likelihood score. It is the quadratic component that connects residual variant scores to the SKAT statistic, but the variance-component likelihood score also has a trace term and a scale factor. Keeping these objects distinct clarifies both the derivation and the role of null calibration.