Kernel Tests and SKAT Methods · Tutorial 446 of 989
SKAT for Binary Traits
Build and interpret the logistic score test used by SKAT for case-control phenotypes, including the role of covariates, genotype weights, and fitted probabilities.
Advanced
11 min read
What You'll Learn
Specify the logistic null model for a binary phenotype and case-control sample
Construct the weighted SKAT score statistic from genotypes and fitted residuals
Account for nuisance covariates in the score covariance and null distribution
Calculate a one-variant score statistic and its asymptotic p-value
Explain how case-control balance changes fitted probabilities and score information
Avoid interpreting sample fitted probabilities as population disease risks
Why Binary-Trait SKAT Needs a Logistic Null Model
For a quantitative phenotype, the null model describes a conditional mean and its residual variation. A binary case-control phenotype instead has a probability of being a case, so SKAT uses a logistic regression null model. The fitted probabilities determine the residuals in the score statistic and, through their Bernoulli variances, how much information the sample contributes.
The null model includes covariates but excludes the genetic variants being tested. SKAT then asks whether the variant scores, taken together and weighted, are unusually large under that fitted null. This is a score-based test: the genetic effects are evaluated at zero without fitting a separate effect estimate for every variant under the alternative.
Definition: Let \(Y_i\) be the binary phenotype for participant \(i\), with \(Y_i=1\) for a case and \(Y_i=0\) for a control. Let \(C_i\) be the vector of covariates, including an intercept, and let \(\mu_i=P(Y_i=1\mid C_i)\) under the null. The logistic null model is
$$
\operatorname{logit}(\mu_i)=\log\left(\frac{\mu_i}{1-\mu_i}\right)=C_i^\top\gamma,
$$
where \(\gamma\) is estimated using the phenotype and covariates, but not the tested genotypes. For \(m\) variants, \(G\) is the \(n\times m\) genotype matrix, with dosage \(G_{ij}\) for participant \(i\) and variant \(j\). The null hypothesis is that all tested genetic effects are zero, conditional on the covariates.
The model assumes that, conditional on the covariates under the null, each phenotype has a Bernoulli distribution with probability \(\mu_i\). Standard score-test calculations also rely on independent participants and on a sufficiently well-behaved fitted null model. Relatedness, clustering, or other dependence requires an appropriate model for that structure; the ordinary independent-participant formula should not simply be applied unchanged.
From Logistic Residuals to the SKAT Statistic
Let \(\hat{\mu}_i\) be the probability fitted from the logistic null model, and let \(Y\) be the vector of observed phenotypes. The residual vector is \(Y-\hat{\mu}\). For variant \(j\), the score for a genetic effect at zero is
$$
U_j=\sum_{i=1}^n G_{ij}(Y_i-\hat{\mu}_i).
$$
Thus \(U_j\) compares the observed case pattern among carriers with the pattern expected under the fitted null. Put the scores into the vector \(U=G^\top(Y-\hat{\mu})\). Let \(W\) be a diagonal matrix of prespecified nonnegative variant weights. The SKAT statistic is
$$
Q=U^\top W^2U
=(Y-\hat{\mu})^\top GW^2G^\top(Y-\hat{\mu}).
$$
A larger \(Q\) indicates greater weighted squared score evidence against the null, but its scale depends on the genotype data, weights, and fitted probabilities. It is not compared with a single universal cutoff. Its null distribution is generally a weighted sum of independent chi-square variables, with weights determined by the eigenvalues of the score covariance matrix after accounting for fitted nuisance parameters.
For clarity, write \(V=\operatorname{diag}\{\hat{\mu}_i(1-\hat{\mu}_i)\}\), and let \(C\) be the \(n\times q\) covariate design matrix. The covariance of the efficient genetic score under the null is \(G^\top P G\), where
$$
P=V-VC(C^\top VC)^{-1}C^\top V.
$$
Accordingly, the mixture weights for \(Q\) are the eigenvalues of \(W G^\top P G W\). This adjustment matters because the fitted covariate effects are nuisance parameters: uncertainty in those estimates affects the score covariance. The logistic variance contribution for participant \(i\) is \(\hat{\mu}_i(1-\hat{\mu}_i)\), which is largest at \(\hat{\mu}_i=0.5\) and smaller when the fitted probability is close to zero or one.
Why the Covariate Adjustment Appears
The covariance expression can be understood as the part of the genetic score variation that remains after accounting for the nuisance covariate score. The argument is a standard efficient-score calculation. Under the null, the working covariance of the phenotype residuals is \(V\). The score for the covariates has covariance \(C^\top VC\), and its covariance with the genetic score is \(G^\top VC\).
Removing the component predictable from the covariate score subtracts its covariance contribution from the unadjusted genetic covariance \(G^\top VG\):
This is the Schur-complement form of the nuisance-parameter adjustment. In practice, the null model is fitted first, the genetic score is evaluated at zero, and the resulting covariance is used to obtain the mixture distribution. One does not treat the fitted probabilities as known constants when calculating that covariance.
Worked Example: A One-Variant Logistic Score Test
Suppose a hypothetical sample has eight unrelated participants, four cases and four controls. The null model contains only an intercept, so its fitted probability is \(\hat{\mu}_i=4/8=0.5\) for everyone. A hypothetical variant has dosage 1 for two cases and one control, and dosage 0 for everyone else. Use weight \(w=1\).
Score. The total dosage among cases is 2 and the total dosage in the sample is 3. Therefore
$$
U=\sum_i G_i(Y_i-\hat{\mu}_i)
=2-0.5(3)=0.5.
$$
Score variance. With an intercept-only model, the covariate adjustment subtracts the fitted mean component. Here \(V=0.25I\), and the genotype sum is 3, so
$$
G^\top PG
=0.25(3)-\frac{0.25^2}{8(0.25)}(3)^2
=0.75-0.28125
=0.46875.
$$
Statistic and reference distribution. The statistic is \(Q=w^2U^2=0.5^2=0.25\). With one variant, its asymptotic null distribution is \(0.46875\chi^2_1\). The standardized value is \(0.25/0.46875=0.5333\), giving an approximate upper-tail p-value of 0.465. This small hypothetical sample is used to show the calculation, not to claim that the asymptotic approximation is reliable at this sample size.
Fitted Probabilities and Case-Control Balance
In a case-control sample, the proportion of cases is set by recruitment and need not equal disease prevalence in the source population. The null model must nevertheless describe the sampled participants used for the test. In the simplest intercept-only model, the fitted probability is the sample case fraction. With covariates, fitted probabilities vary across participants according to the logistic regression.
The fitted probabilities enter the score through \(Y_i-\hat{\mu}_i\), and enter its variance through \(\hat{\mu}_i(1-\hat{\mu}_i)\). Therefore two samples with the same genotypes and case count can yield different scores and score covariances if their covariate patterns, fitted probabilities, or case-control balance differ. In particular, a fitted probability near 0.5 contributes more Bernoulli variance than one near 0.05 or 0.95. That does not by itself establish that a particular sample design is always more powerful: information also depends on genotype distribution, covariate adjustment, and how cases and controls are sampled.
Under conventional case-control sampling in which selection depends on phenotype but not genotype within phenotype and covariate strata, logistic regression slopes can be used to assess conditional odds-ratio associations, while the fitted intercept reflects the sampled case fraction rather than population prevalence. The SKAT test concerns genetic association conditional on the model; its sample fitted probabilities should not be presented as population disease risks without a suitable prevalence correction and sampling model.
Worked Example: A Weighted Multi-Variant Statistic
Consider a hypothetical intercept-only sample of eight participants with four cases and four controls. The genotypes are shown below. The first variant has weight 1 and the second has weight 0.5.
Participant
Phenotype \(Y\)
Variant 1
Variant 2
1
Case
1
0
2
Case
1
1
3
Case
0
0
4
Case
0
0
5
Control
1
1
6
Control
0
1
7
Control
0
1
8
Control
0
0
The fitted probability is 0.5 for every participant. Variant 1 has total dosage 3, including 2 among cases; variant 2 has total dosage 4, including 1 among cases. Thus
$$
U_1=2-0.5(3)=0.5,\qquad
U_2=1-0.5(4)=-1.
$$
The negative second score indicates fewer observed case dosages than expected under the fitted null; squaring means either direction contributes to SKAT. With \(W=\operatorname{diag}(1,0.5)\),
$$
Q=U^\top W^2U
=1^2(0.5)^2+0.5^2(-1)^2
=0.25+0.25
=0.5.
$$
The statistic is not enough by itself to decide significance. Its mixture reference distribution depends on the covariance between the two scores as well as their individual variances. In this example, the adjusted covariance matrix is
$$
G^\top PG=
\begin{pmatrix}
0.46875 & 0.125\\
0.125 & 0.5
\end{pmatrix},
$$
so the mixture weights are the eigenvalues of \(W G^\top PG W\), approximately 0.4798 and 0.1140. A p-value must be calculated from the resulting weighted chi-square distribution; comparing \(Q=0.5\) directly with a one-degree-of-freedom chi-square cutoff would use the wrong reference distribution.
Worked Example: An Imbalanced Case-Control Sample
Suppose a hypothetical sample has 100 unrelated participants: 20 cases and 80 controls. In an intercept-only null model, \(\hat{\mu}_i=0.2\). A variant is present in 10 participants, of whom 3 are cases and 7 are controls, with dosage 1 for carriers and 0 otherwise.
The score and its adjusted variance are
$$
U=3-0.2(10)=1,
\qquad
G^\top PG
=0.2(0.8)\left\{10-\frac{10^2}{100}\right\}
=0.16(9)=1.44.
$$
With weight 1, \(Q=1^2=1\), and the one-variant asymptotic reference is \(1.44\chi^2_1\). The standardized statistic is \(1/1.44=0.6944\), yielding an approximate p-value of 0.405. The fitted probability of 0.2 describes the sample under this null model; it should not be interpreted as the population disease prevalence merely because it was estimated from the sample.
Common Mistakes and Practical Interpretation
Using a linear-probability residual model. A binary phenotype has Bernoulli variance \(\mu_i(1-\mu_i)\), not a constant Gaussian residual variance. Fit the logistic null and use its fitted probabilities in the score and covariance.
Ignoring nuisance covariates in the variance. Fitting covariates and then using an unadjusted genotype covariance treats their estimated effects as if they had no uncertainty. Use the adjusted covariance \(G^\top PG\).
Using the same chi-square reference for every SKAT statistic. The multi-variant null distribution is generally a mixture determined by the genotype weights and adjusted score covariance. Calculate its tail probability using the mixture eigenvalues.
Equating a balanced sample with population prevalence. A 50:50 case-control sample may be deliberately assembled. Its fitted intercept does not automatically estimate population risk or prevalence.
Calling a significant score a causal genetic effect. A small p-value is evidence against the specified no-association null under the model and sampling assumptions. It does not establish causation, identify a causal variant, or show that every tested variant has an effect.
Key takeaway: Binary-trait SKAT fits a covariate-adjusted logistic null, scores each variant using \(Y_i-\hat{\mu}_i\), and evaluates a weighted sum of squared scores against a mixture distribution. The fitted probabilities—and especially their Bernoulli variances—are essential to valid inference in the case-control sample.
Check Your Understanding
Answer these questions using the definitions and calculations above.
What quantities enter the logistic null model, and which quantities are excluded when testing a set of variants?
For one variant with case dosage total 4, sample dosage total 10, and fitted probability 0.3 for every participant, calculate its score.
Why does the covariance used for SKAT include an adjustment for the covariate design matrix?
What determines the mixture chi-square weights for a multi-variant binary-trait SKAT test?
Why should the fitted intercept from an artificially balanced case-control sample not automatically be interpreted as population disease prevalence?