Constructing the SKAT-O Hybrid Test

Why Combine Burden and SKAT?

Rare-variant association tests can be more or less sensitive to the pattern of variant effects. A burden test aggregates weighted variant scores before squaring the aggregate. It can be effective when effects point in a common direction, but effects in opposite directions can cancel. SKAT instead squares and sums weighted scores variant by variant, so opposing contributions do not cancel in the same way. Its sensitivity can be better when effect directions vary, though it may be less powerful than a burden test when effects are aligned.

SKAT-O constructs a family of score statistics between these two approaches. Its parameter $\rho$ controls how much the statistic emphasizes a burden-like aggregate versus the SKAT-like sum of squared scores. The family provides a way to address uncertainty about effect directions without treating one pattern as known in advance.

Definition and Score-Test Setup

Suppose a gene contains $m$ variants. Let $U=(U_1,\ldots,U_m)'$ be the vector of variant score statistics under the null model, and let $V$ be its covariance matrix under that model. A score is evidence in the direction of association for its variant; the entries of $U$ need not have the same sign. Let $W$ be a fixed diagonal matrix of nonnegative variant weights, with diagonal entries $w_1,\ldots,w_m$. The weights may reflect a prespecified choice of variant weighting.

The usual SKAT and burden-like score statistics in this notation are

Q_{\mathrm{SKAT}}=U'W^2U,\qquad Q_{\mathrm{burden}}=(\mathbf{1}'WU)^2,

where $\mathbf{1}$ is the $m$-vector of ones. The first statistic adds weighted squared scores. The second first adds weighted scores, then squares their sum.

Definition: For $0\leq\rho\leq1$, the SKAT-O score statistic at $\rho$ is $$ Q_\rho=(1-\rho)Q_{\mathrm{SKAT}}+\rho Q_{\mathrm{burden}} =U'W\big[(1-\rho)I+\rho\mathbf{1}\mathbf{1}'\big]WU. $$ At $\rho=0$, it is the SKAT statistic; at $\rho=1$, it is the burden-like statistic.

The restriction to $[0,1]$ gives a convex combination of the two statistics and a positive semidefinite matrix in the quadratic form. This construction is not an average of variant-level estimated effects. It is a family of quadratic forms in the score vector, with different forms corresponding to different assumptions about how variant effects relate to one another.

How the Correlation Parameter Represents Effect Patterns

One motivation for the family is a random-effects model for the variant effects. Let $\beta=(\beta_1,\ldots,\beta_m)'$ denote the effect vector and suppose, as a model for alternatives, that

\beta\sim N\!\left(0,\tau^2W\big[(1-\rho)I+\rho\mathbf{1}\mathbf{1}'\big]W\right), \qquad \tau^2\geq0,\quad 0\leq\rho\leq1.

Here the weights scale the effects, and $\rho$ controls their dependence. If the weights are positive, the standardized effects $\beta_j/w_j$ have variance $\tau^2$ and pairwise correlation $\rho$. A zero weight gives that variant zero modeled effect variance, so its standardized effect and correlation are not defined. At $\rho=0$, effects for variants with nonzero weights are independent under this model. At $\rho=1$, they share a common random direction after scaling by their weights. The effect covariance then has rank at most one; it has rank one when $\tau^2>0$ and at least one weight is nonzero, and otherwise has rank zero.

The score statistic associated with this covariance structure uses the same matrix between $W$ and $W$ in its quadratic form. This correspondence explains why changing $\rho$ changes the balance of the test. It does not mean that the alternative random-effects model is known to be true; rather, it gives an interpretable route to constructing a family of tests.

Worked Example: A Middle Value of the Parameter

Consider three hypothetical variant scores $U=(2,-1,1)'$ and weights $W=\operatorname{diag}(1,2,1)$. Suppose these scores came from a correctly specified null model, so the score vector and its covariance are appropriate for the score-test calculation. The weighted scores are $WU=(2,-2,1)'$.

The SKAT statistic is

Q_{\mathrm{SKAT}}=2^2+(-2)^2+1^2=9.

The weighted aggregate is $2-2+1=1$, so

Q_{\mathrm{burden}}=(1)^2=1.

At $\rho=0.4$, the hybrid statistic is

Q_{0.4}=(1-0.4)(9)+0.4(1)=5.4+0.4=5.8.

This is a statistic, not a p-value. Its value depends on the chosen $\rho$; assessing evidence requires comparing it with its null distribution at that same $\rho$.

Why the Formula Interpolates

The interpolation follows directly by expanding the quadratic form. Let $A_\rho=W[(1-\rho)I+\rho\mathbf{1}\mathbf{1}']W$. Then

\begin{aligned} U'A_\rho U &=(1-\rho)U'W^2U+\rho U'W\mathbf{1}\mathbf{1}'WU\\ &=(1-\rho)Q_{\mathrm{SKAT}}+\rho(\mathbf{1}'WU)^2. \end{aligned}

The second equality uses $U'W\mathbf{1}\mathbf{1}'WU=(\mathbf{1}'WU)(\mathbf{1}'WU)$, a scalar squared. Thus the two endpoint statistics arise from the same quadratic-form construction: the identity matrix selects separate squared scores, while $\mathbf{1}\mathbf{1}'$ makes every weighted score contribute to one aggregate.

The matrix is positive semidefinite for every allowed $\rho$. For any vector $x$,

x'\big[(1-\rho)I+\rho\mathbf{1}\mathbf{1}'\big]x =(1-\rho)x'x+\rho(\mathbf{1}'x)^2\geq0.

Since $W$ is diagonal, $U'A_\rho U$ is also nonnegative. This property is important for interpreting the statistic as a variance-component score statistic and for deriving its null distribution.

Worked Example: Opposing Scores and Cancellation

Take two hypothetical variants with $U=(3,-2)'$ and equal weights $W=I$. Assume the scores were obtained under a valid null model with covariance accounted for in the subsequent calibration. The SKAT statistic is

Q_{\mathrm{SKAT}}=3^2+(-2)^2=13.

The burden aggregate is $3+(-2)=1$, giving $Q_{\mathrm{burden}}=1^2=1$. At $\rho=0.8$,

Q_{0.8}=0.2(13)+0.8(1)=2.6+0.8=3.4.

The difference illustrates cancellation in the aggregate: the negative score offsets much of the positive score for the burden-like component, while SKAT retains both squared contributions. It does not establish which test is more significant; their null distributions must also be considered.

Null Distribution and Calibration at a Fixed Parameter

For inference, the null model must be correctly specified, including the treatment of covariates and the covariance of the score vector. Under the null and standard large-sample score-test conditions, $U$ is approximately multivariate normal with mean zero and covariance $V$. At a fixed $\rho$, $Q_\rho=U'A_\rho U$, where $A_\rho=W[(1-\rho)I+\rho\mathbf{1}\mathbf{1}']W$. Its null distribution is generally a weighted sum of independent chi-square variables, not necessarily a single chi-square distribution.

In particular, let $\lambda_1,\ldots,\lambda_m$ be the eigenvalues of $V^{1/2}A_\rho V^{1/2}$. Under the stated normal approximation,

Q_\rho\ \overset{d}{\approx}\ \sum_{j=1}^{m}\lambda_j Z_j^2 =\sum_{j=1}^{m}\lambda_j\chi^2_{1,j},

where the $Z_j$ are independent standard normal variables, and zero eigenvalues contribute zero. This follows by writing $U=V^{1/2}Z$, with $Z$ having independent standard normal entries, and diagonalizing the symmetric matrix $V^{1/2}A_\rho V^{1/2}$. This is the same quadratic-form principle used for SKAT; the eigenvalues change with $\rho$.

At $\rho=1$, the statistic has the single squared-aggregate form $Q_{\mathrm{burden}}$, but its null variance still depends on the covariance among variant scores. For example, if the aggregate score is $B=\mathbf{1}'WU$, then under the null its variance is $\mathbf{1}'WVW\mathbf{1}$. Ignoring correlations among scores can therefore give an incorrect reference distribution even when the observed aggregate is computed correctly.

Worked Example: Two Variants Under an Illustrative Null Covariance

For a transparent hypothetical calculation, suppose $U=(1,1)'$, $W=I$, and the null covariance is $V=I$. These assumptions make the scores independent standard normal variables under the null. The component statistics are $Q_{\mathrm{SKAT}}=1^2+1^2=2$ and $Q_{\mathrm{burden}}=(1+1)^2=4$. At $\rho=0.5$,

Q_{0.5}=0.5(2)+0.5(4)=3.

Here $A_{0.5}=\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}$, whose eigenvalues are $1.5$ and $0.5$. Consequently, under this illustrative null,

Q_{0.5}\ \overset{d}{=}\ 1.5\chi^2_{1,1}+0.5\chi^2_{1,2}.

The observed value $3$ alone does not determine a p-value. The weighted chi-square reference distribution is part of the calculation, and the simple $V=I$ assumption is specific to this example.

From a Fixed Statistic to the SKAT-O Family

A fixed value of $\rho$ corresponds to one specified balance between the two components. In practice, the effect-direction pattern is often uncertain, so SKAT-O considers a family of values: small values emphasize heterogeneous effects, while large values emphasize a common weighted direction. The family can therefore adapt to patterns for which either endpoint, or an intermediate statistic, is relatively informative.

Considering several values creates an additional inference issue. The most favorable statistic in a family cannot simply be reported with its fixed-$\rho$ p-value: selecting among correlated tests can make that p-value too small. The search must be calibrated to account for the dependence among the statistics. The construction here defines the family and its component statistics; calibration of the parameter search is a separate step.

Worked Example: Aligned Scores Favor the Aggregate Component

Suppose two hypothetical variants have scores $U=(3,2)'$, equal weights, and a correctly specified null model for calculating their covariance. Then $Q_{\mathrm{SKAT}}=3^2+2^2=13$, whereas $Q_{\mathrm{burden}}=(3+2)^2=25$. At $\rho=0.6$,

Q_{0.6}=0.4(13)+0.6(25)=5.2+15=20.2.

The weighted scores point in the same direction, so aggregation reinforces rather than cancels them. The burden-like statistic is larger in this example, and the hybrid statistic moves toward it as $\rho$ increases. This comparison describes the observed statistics only; significance still depends on each statistic's null distribution and, if $\rho$ is selected, on calibration for that selection.

Common Mistakes and Practical Checks

Assuming the hybrid is a new endpoint unrelated to the components. It is exactly $(1-\rho)Q_{\mathrm{SKAT}}+\rho Q_{\mathrm{burden}}$. Checking the two endpoints is a useful way to catch a reversed or misapplied $\rho$.
Calling $\rho$ the fraction of variants with effects. In the random-effects formulation, it controls dependence among the standardized effects; it is not a proportion of causal variants.
Treating a larger statistic as automatically stronger evidence. The null distribution changes with $A_\rho$, so comparisons of raw $Q_\rho$ values across parameter values do not by themselves compare p-values.
Reporting the smallest fixed-$\rho$ p-value as the family p-value. Searching across correlated tests requires calibration for that search.
Assuming the effect covariance at $\rho=1$ must have rank one. It has rank at most one. It has rank one only if $\tau^2>0$ and at least one weight is nonzero.

Key takeaway: SKAT-O connects a sum of squared weighted scores to the square of their weighted sum through $Q_\rho=(1-\rho)Q_{\mathrm{SKAT}}+\rho Q_{\mathrm{burden}}$. The parameter represents an assumption about effect-direction dependence; valid inference also requires the appropriate null covariance and, when searching across $\rho$, calibration for that search.

Check Your Understanding

Answer each question using the definitions and assumptions in the tutorial.

For what values of $\rho$ does $Q_\rho$ equal the SKAT statistic and the burden-like statistic?
With weighted scores $WU=(2,-2,1)'$, calculate $Q_{\mathrm{SKAT}}$ and $Q_{\mathrm{burden}}$. Which statistic reflects cancellation?
In the random-effects model, what does $\rho$ describe, and what does it not describe?
Why is the null distribution of a fixed-$\rho$ statistic generally a weighted sum of chi-square variables?
Why can the smallest p-value among several $\rho$ values not automatically be reported as the calibrated SKAT-O p-value?