The Logistic Score Test Derivation

Why Test the SNP Without Fitting Its Effect?

For a binary trait, an additive SNP model relates the probability of being a case to genotype dosage and covariates through a logistic regression. A score test asks whether the data provide evidence for a nonzero per-allele log-odds effect while evaluating the model at the null hypothesis that this effect is zero. Its practical feature is that it requires fitting the covariate model under the null, not fitting a separate alternative model for every SNP.

The test is valid only under stated modeling and large-sample conditions. In particular, the covariates must be accounted for in both the null fit and the information calculation. The raw genotype score alone is not enough: its variance must reflect how much genotype variation remains after accounting for covariates.

Model and Null Fit

For \(n\) unrelated individuals, let \(Y_i\) be the binary outcome for individual \(i\), with \(Y_i=1\) for a case and \(Y_i=0\) for a control. Let \(G_i\) be the SNP dosage, coded as 0, 1, or 2 copies of the effect allele, and let \(C_i\) be a vector of covariates. The additive logistic regression model is

Y_i\mid G_i,C_i\sim\operatorname{Bernoulli}(\pi_i),\qquad \operatorname{logit}(\pi_i)=\gamma_0+\gamma^\mathsf{T}C_i+\beta G_i, \qquad \operatorname{logit}(\pi)=\log\left(\frac{\pi}{1-\pi}\right).

Here \(\beta\) is the change in log odds associated with one additional effect-allele copy, holding the covariates fixed. The hypothesis tested is \(H_0:\beta=0\) against \(H_a:\beta\ne 0\). Under \(H_0\), genotype is omitted from the linear predictor, and the nuisance parameters \(\gamma_0,\gamma\) are estimated by fitting the logistic regression of \(Y_i\) on the covariates alone. Denote this null-model estimate by \(\hat\gamma_0\), using \(\hat\gamma_0\) to represent the full vector of null-model coefficients, including the intercept. The fitted null probability for individual \(i\) is \(\hat\pi_i\).

Definition: The genotype score evaluated at the null fit is \(U_\beta=\sum_{i=1}^n G_i(Y_i-\hat\pi_i)\). It measures the direction and size of the genotype-associated residual signal under \(H_0\). The score test divides its square by the efficient information for \(\beta\), which accounts for estimating the covariate effects.

For independent Bernoulli observations, the log-likelihood, apart from no relevant constants, is

\ell(\gamma_0,\gamma,\beta) =\sum_{i=1}^n\left\{Y_i\log(\pi_i)+(1-Y_i)\log(1-\pi_i)\right\}.

Differentiating with respect to \(\beta\) gives the genotype score at any parameter value, \(\partial\ell/\partial\beta=\sum_i G_i(Y_i-\pi_i)\). Evaluating at \(H_0\) and the null-model estimates yields \(U_\beta\). The null fit also satisfies the nuisance score equations: the residuals \(Y_i-\hat\pi_i\) have zero weighted-free sums against the intercept and each included covariate column. This fact is important when deriving the variance of the genotype score.

Efficient Information and the Score Statistic

Let \(D\) be the \(n\)-by-\(q\) design matrix for the null model, including its intercept and covariate columns. Let \(g\) be the \(n\)-vector of genotype dosages, and let \(W\) be the diagonal matrix with entries \(w_i=\hat\pi_i(1-\hat\pi_i)\). At the null fit, the observed information for the full parameter vector has the same block form as the expected information for logistic regression:

I= \begin{pmatrix} D^\mathsf{T}WD & D^\mathsf{T}Wg\\ g^\mathsf{T}WD & g^\mathsf{T}Wg \end{pmatrix}.

The lower-right entry alone, \(g^\mathsf{T}Wg\), treats the nuisance coefficients as known. They are estimated, however, so the information available for \(\beta\) after accounting for those estimates is the Schur complement of the nuisance block:

I_{\mathrm{eff}} =g^\mathsf{T}Wg -g^\mathsf{T}WD(D^\mathsf{T}WD)^{-1}D^\mathsf{T}Wg.

Equivalently, define the weighted residualized genotype \(\tilde g=g-D(D^\mathsf{T}WD)^{-1}D^\mathsf{T}Wg\). Then \(I_{\mathrm{eff}}=\tilde g^\mathsf{T}W\tilde g\). This expression shows that the information comes from genotype variation not explained by the covariates, with each observation weighted according to its fitted Bernoulli variance.

Theorem: Under \(H_0\), if the logistic model is correctly specified, observations are independent, the nuisance parameters are identifiable and estimable, and standard large-sample regularity conditions hold, the score statistic \(S=U_\beta^2/I_{\mathrm{eff}}\) has an asymptotic chi-square distribution with one degree of freedom. Its p-value is the upper-tail probability \(P(\chi^2_1\ge S)\). For a two-sided test at significance level \(\alpha\), reject \(H_0\) when \(S\) exceeds the \(1-\alpha\) quantile of \(\chi^2_1\).

A short justification explains why the nuisance adjustment enters through information. A Taylor expansion of the nuisance score equations around the true null nuisance parameter gives the first-order change in its estimate as the inverse nuisance information times the nuisance score. Substituting that change into the genotype score removes the part of the genotype signal explained by the nuisance score. The resulting variance is the Schur complement above. Under the null and regularity conditions, this adjusted score divided by the square root of its variance is asymptotically standard normal; squaring it gives the stated \(\chi^2_1\) reference distribution.

In this calculation, the null nuisance score is zero at the fitted null model. Consequently, the efficient score's numerator is the raw \(U_\beta\) evaluated at that fit, even though its variance must be adjusted for nuisance estimation. The test is one degree of freedom because it tests one additive genotype coefficient. It tests association conditional on the included covariates under the specified model; by itself, it does not establish a causal effect.

Worked Example: Score Test With an Intercept Only

Consider six hypothetical individuals with genotypes \(G=(0,1,1,2,0,2)\) and binary outcomes \(Y=(0,1,1,0,0,1)\). Use an intercept-only null model. There are three cases out of six, so the null-model estimate is \(\hat\pi_i=3/6=0.5\) for every individual. This fit is finite, and the genotype varies in the sample.

The genotype score is \(U_\beta=\sum_iG_iY_i-\hat\pi\sum_iG_i\). Here, \(\sum_iG_iY_i=0+1+1+0+0+2=4\), and \(\sum_iG_i=6\), so

U_\beta=4-(0.5)(6)=1.

For an intercept-only model, efficient information is the weighted sum of squared deviations from the genotype mean. The mean dosage is \(\bar G=6/6=1\), and \(\sum_i(G_i-\bar G)^2=1+0+0+1+1+1=4\). Since each \(w_i=0.5(1-0.5)=0.25\),

I_{\mathrm{eff}}=0.25(4)=1,\qquad S=\frac{1^2}{1}=1.

The asymptotic p-value is \(P(\chi^2_1\ge1)\approx0.3173\). This calculation does not provide statistically significant evidence against \(H_0\) at \(\alpha=0.05\). With only six individuals, however, the large-sample approximation is not dependable; the example illustrates the arithmetic, not a reliable small-sample conclusion.

Worked Example: Covariate-Adjusted Score and Information

Suppose six hypothetical individuals fall into two covariate groups, with \(C_i=0\) for the first three and \(C_i=1\) for the last three. In each group the dosages are \(0,1,2\). The outcomes are \((1,0,0)\) in the first group and \((1,1,0)\) in the second. Fit the null logistic model with an intercept and \(C_i\), but no genotype.

The observed case fractions are \(1/3\) and \(2/3\), respectively. The null model with an intercept and a binary group indicator fits these fractions exactly, so its fitted probabilities are \(1/3\) for the first group and \(2/3\) for the second. The null fit is finite, and both groups contain outcome variation. In each group the dosages vary, so genotype is not entirely determined by group.

In the first group, \(\sum G_iY_i=0\) and \(\sum G_i=3\); its contribution to the score is \(0-(1/3)(3)=-1\). In the second group, \(\sum G_iY_i=1\) and \(\sum G_i=3\); its contribution is \(1-(2/3)(3)=-1\). Therefore \(U_\beta=-2\).

Within each group, the fitted weights are \(w_i=(1/3)(2/3)=2/9\). The weighted projection of genotype onto the intercept and group indicator gives the group mean dosage, 1. Thus the residualized genotypes in each group are \((-1,0,1)\). The efficient information is

I_{\mathrm{eff}} =\frac{2}{9}\left\{(-1)^2+0^2+1^2\right\} +\frac{2}{9}\left\{(-1)^2+0^2+1^2\right\} =\frac{8}{9}.

The score statistic is \(S=(-2)^2/(8/9)=4.5\), giving \(P(\chi^2_1\ge4.5)\approx0.0339\). The asymptotic test rejects \(H_0\) at the 0.05 level. This is still a six-person example, so the approximation may be poor; in a real analysis, the result should not be treated as reliable solely because its approximate p-value is below 0.05.

Worked Example: Reversing the Effect Allele

Use the adjusted example above, but encode genotype as \(G_i^*=2-G_i\), the dosage of the other allele. In the model with an intercept and the same covariate, the constant 2 is already represented by the intercept. The residualized genotype therefore changes sign: \(\tilde g^*=-\tilde g\). The adjusted information remains \(I_{\mathrm{eff}}^*=8/9\).

Because the null residuals sum to zero (the intercept score equation), the new score is \(U_\beta^*=\sum_i(2-G_i)(Y_i-\hat\pi_i)=2\sum_i(Y_i-\hat\pi_i)-U_\beta=-(-2)=2\). Thus \(S^*=(2)^2/(8/9)=4.5\), the same two-sided test statistic and p-value as before. The sign of the score, like the sign of the per-allele coefficient, depends on which allele is counted; the two-sided evidence does not change under this consistent recoding.

Assumptions, Diagnostics, and Practical Limits

The reference distribution is an approximation, not a guarantee for every data set. Its use depends on independent observations, a suitable logistic mean model under \(H_0\), estimable nuisance parameters, and enough information for the score to have an approximately normal distribution. The null-model design must have full column rank, and the genotype must retain variation after its weighted projection onto the covariates. If \(I_{\mathrm{eff}}=0\), there is no adjusted genotype information and the statistic is undefined.

Before interpreting the result, verify that the null model converged and has finite coefficient estimates. Check that the genotype was coded as intended, that the covariate design is identifiable, and that the efficient information is positive and not nearly zero. Very small fitted Bernoulli variances, extreme fitted probabilities, or a small effective amount of residualized genotype variation can make the approximation unstable. A large raw genotype variation does not guarantee a large \(I_{\mathrm{eff}}\): covariates may explain nearly all of that variation.

The score test is useful when a null model can be fit once and tested against many genotype vectors, as in a single-variant scan. A direct alternative is to fit the full logistic regression including genotype and estimate its coefficient, which requires fitting the alternative model. The score calculation avoids that alternative fit, but still relies on the null model and the score-test approximation. An unadjusted genotype test is not a substitute when the analysis calls for covariate adjustment: it uses a different model and can reflect covariate-associated differences rather than the conditional genotype association of interest.

Common Mistakes

Using \(g^\mathsf{T}Wg\) as the adjusted information. This ignores estimation of covariate effects. Use the Schur complement or the weighted residualized-genotype expression.
Fitting the null model without the intercept or required covariates. The null fit must represent the covariate-adjusted model being tested, with a properly specified design.
Calling the score statistic a signed effect estimate. The score has a sign, but the test statistic \(S\) is nonnegative and tests evidence against \(\beta=0\); it is not an estimate of the per-allele log odds ratio.
Treating an asymptotic p-value as exact in a tiny or unstable sample. The chi-square calibration relies on large-sample regularity and adequate information.
Interpreting rejection as proof of causation. The result is evidence of conditional association under the model, not proof that the SNP causes the outcome.

Key takeaway: Fit the logistic model under \(H_0:\beta=0\), calculate \(U_\beta=\sum_iG_i(Y_i-\hat\pi_i)\), and scale its square by the efficient information \(I_{\mathrm{eff}}\), which removes the genotype variation explained by covariates. Under regular large-sample conditions, \(U_\beta^2/I_{\mathrm{eff}}\) is compared with \(\chi^2_1\).

Check Your Understanding

Answer each question, then check the solutions below.

In a null logistic model with an intercept and covariates, what is the genotype score evaluated at the null fit?
Why is \(g^\mathsf{T}Wg\) generally not the efficient information for the SNP coefficient?
Suppose \(U_\beta=3\) and \(I_{\mathrm{eff}}=6\). Calculate the score statistic and give its approximate chi-square reference distribution.
What happens to the two-sided score statistic if the effect allele is reversed, with dosage recoded from \(G_i\) to \(2-G_i\)?

Solutions. (1) \(U_\beta=\sum_iG_i(Y_i-\hat\pi_i)\), with \(\hat\pi_i\) obtained from the null model containing the intercept and covariates. (2) That expression treats nuisance coefficients as known and does not remove genotype variation explained by the covariates; the efficient information is the Schur complement, equivalently the weighted sum of squared residualized dosages. (3) \(S=3^2/6=1.5\), compared with a \(\chi^2\) distribution with one degree of freedom; its upper-tail p-value is about 0.2207. (4) The score changes sign and the efficient information is unchanged, so its square divided by the information, and hence the two-sided test statistic, is unchanged.