Why Derive the Score Statistic?
A mixed-model association test asks whether a variant's dosage is associated with a quantitative trait while accounting for covariance among individuals. Relatedness changes the covariance of the observations: a residual in one individual cannot always be treated as independent evidence from a residual in another. The score statistic makes this adjustment explicit. Its numerator is not simply a raw dosage-residual correlation; it is a correlation with covariance-adjusted residuals.
The derivation below uses an intercept-only mean model to isolate the variant score. The variant dosage \(g_i\) is treated as fixed, and \(\beta\) is its effect on the trait per additional allele copy. The model is Gaussian, and its covariance matrix is assumed known and positive definite for the exact distributional result. In practice, variance components are usually estimated under the null and substituted, giving a plug-in test whose calibration is generally approximate.
In a linear mixed model, \(V\) can represent both genetic relatedness and residual variation. The earlier material on covariance decomposition explains how such a covariance can arise from random genetic effects. Here, \(V\) is the object needed for the score derivation; write \(W=V^{-1}\) for its precision matrix.
Removing the Null Intercept
Under \(H_0\), the generalized least-squares estimate of the intercept is $$ \hat{\mu}_0=\frac{\mathbf{1}^{\mathsf T}Wy}{\mathbf{1}^{\mathsf T}W\mathbf{1}}. $$ The corresponding covariance-adjusted residual vector is \(Py\), where $$ P=W-W\mathbf{1}(\mathbf{1}^{\mathsf T}W\mathbf{1})^{-1}\mathbf{1}^{\mathsf T}W. $$ This is a precision-weighted residual vector: it first removes the fitted intercept and then weights the remaining residuals by \(W\). In general, \(Py\) is not the ordinary residual vector \(y-\mathbf{1}\hat{\mu}_0\).
The word efficient matters: \(I\) is not merely \(g^{\mathsf T}Wg\). The latter treats the intercept as if it were known. The subtraction in \(P\) removes the component of the variant dosage that can be explained by the intercept under the covariance-weighted geometry.
Deriving the Score and Its Null Distribution
For fixed \(V\), the Gaussian log-likelihood, apart from a constant, is $$ \ell(\mu,\beta)=-\frac{1}{2}\log|V|-\frac{1}{2}(y-\mathbf{1}\mu-g\beta)^{\mathsf T}W(y-\mathbf{1}\mu-g\beta). $$ Differentiating with respect to \(\beta\) gives \(g^{\mathsf T}W(y-\mathbf{1}\mu-g\beta)\). At \(\beta=0\), evaluate this derivative at the null intercept estimate: $$ g^{\mathsf T}W(y-\mathbf{1}\hat{\mu}_0) =g^{\mathsf T}Wy-\frac{(g^{\mathsf T}W\mathbf{1})(\mathbf{1}^{\mathsf T}Wy)}{\mathbf{1}^{\mathsf T}W\mathbf{1}} =g^{\mathsf T}Py=U. $$ Thus the score is the slope of the log-likelihood in the variant-effect direction at the null fit. A positive value indicates that moving toward a positive \(\beta\) initially increases the likelihood; a negative value indicates the opposite direction.
To obtain the information, the likelihood's information matrix for \((\mu,\beta)\) has entries \(\mathbf{1}^{\mathsf T}W\mathbf{1}\), \(\mathbf{1}^{\mathsf T}Wg\), and \(g^{\mathsf T}Wg\). Removing the intercept's contribution gives the Schur complement $$ g^{\mathsf T}Wg-(g^{\mathsf T}W\mathbf{1})(\mathbf{1}^{\mathsf T}W\mathbf{1})^{-1}(\mathbf{1}^{\mathsf T}Wg)=g^{\mathsf T}Pg=I. $$ This explains why the score is standardized by \(\sqrt I\), rather than by \(\sqrt{g^{\mathsf T}Wg}\).
The null distribution follows directly when \(V\) is known. The matrix \(P\) is symmetric, \(P\mathbf{1}=0\), and \(PVP=P\). The last identity follows by substituting \(V=W^{-1}\) into the definition of \(P\), or by noting that \(P\) is the precision-weighted projection that removes the intercept. Under \(H_0\), \(y=\mathbf{1}\mu+\varepsilon\), so \(U=g^{\mathsf T}P\varepsilon\) has mean zero and variance $$ \operatorname{Var}(U)=g^{\mathsf T}PVPg=g^{\mathsf T}Pg=I. $$ Because \(\varepsilon\) is Gaussian, \(U\) is Gaussian, and therefore \(Z\sim N(0,1)\) and \(Q\sim\chi^2_1\) under the null. With estimated covariance parameters substituted for \(V\), these are reference distributions for a plug-in test, not generally exact finite-sample distributions.
Worked Example: Computing the Score With Unequal Precision
Suppose a hypothetical three-person data set has \(g=(0,1,2)^{\mathsf T}\), \(y=(1,2,4)^{\mathsf T}\), and \(W=\operatorname{diag}(1,0.5,1)\). This positive diagonal precision matrix corresponds to \(V=\operatorname{diag}(1,2,1)\). Treat \(V\) as known and use the intercept-only null model.
First calculate the weighted sums: $$ \mathbf{1}^{\mathsf T}W\mathbf{1}=1+0.5+1=2.5,\quad \mathbf{1}^{\mathsf T}Wy=1+1+4=6,\quad g^{\mathsf T}W\mathbf{1}=0+0.5+2=2.5,\quad g^{\mathsf T}Wy=0+1+8=9. $$ The null intercept estimate is \(\hat{\mu}_0=6/2.5=2.4\). The score is $$ U=g^{\mathsf T}Wy-\frac{(g^{\mathsf T}W\mathbf{1})(\mathbf{1}^{\mathsf T}Wy)}{\mathbf{1}^{\mathsf T}W\mathbf{1}} =9-\frac{(2.5)(6)}{2.5}=3. $$ For the information, \(g^{\mathsf T}Wg=0+0.5+4=4.5\), so $$ I=4.5-\frac{(2.5)^2}{2.5}=2. $$ Consequently, \(Q=3^2/2=4.5\) and \(Z=3/\sqrt{2}\approx2.121\). The two-sided standard normal p-value is \(2P(Z_{\!0}\ge 2.121)\approx0.03389\), where \(Z_{\!0}\sim N(0,1)\). Under the stated known-covariance model, this is evidence against no variant effect at the 0.05 level, but not at the 0.01 level. This small hypothetical example is for calculation practice, not a genome-wide significance claim.
What the Adjusted Residuals Contribute
The same score can be read as \(U=g^{\mathsf T}(Py)\). Each individual's trait residual contributes according to the precision matrix, and the intercept component is removed. With correlated observations, \(W\) need not be diagonal: an individual's adjusted residual can then depend on residuals from other individuals. This is why a raw residual correlation can give a different answer from the mixed-model score.
Worked Example: Reading the Covariance Adjustment
Consider a hypothetical four-person data set with \(g=(0,0,1,1)^{\mathsf T}\), \(y=(2,3,2,5)^{\mathsf T}\), and \(W=\operatorname{diag}(1,0.5,1,0.5)\). Again suppose \(V=W^{-1}\) is known. Here \(\mathbf{1}^{\mathsf T}W\mathbf{1}=3\) and \(\mathbf{1}^{\mathsf T}Wy=8\), so \(\hat{\mu}_0=8/3\).
The unweighted null residuals would be \(y-\mathbf{1}(8/3)\), but the score uses their precision-weighted form: $$ Py=W(y-\mathbf{1}\hat{\mu}_0) =\left(-\frac{2}{3},\frac{1}{6},-\frac{2}{3},\frac{7}{6}\right)^{\mathsf T}. $$ Therefore, $$ U=g^{\mathsf T}Py=-\frac{2}{3}+\frac{7}{6}=\frac{1}{2}. $$ For its information, \(g^{\mathsf T}Wg=1.5\) and \(g^{\mathsf T}W\mathbf{1}=1.5\), giving $$ I=1.5-\frac{(1.5)^2}{3}=0.75. $$ Thus \(Q=0.25/0.75=0.3333\), \(Z=0.5/\sqrt{0.75}\approx0.5774\), and the two-sided normal p-value is approximately \(0.5637\). The positive score points toward a positive dosage effect, but its magnitude is small relative to its null standard deviation.
For comparison, ignoring the weights gives an ordinary intercept-adjusted dosage-residual cross-product of \(g^{\mathsf T}y-(\mathbf{1}^{\mathsf T}g)(\mathbf{1}^{\mathsf T}y)/4=7-(2)(12)/4=1\), not \(0.5\). That raw quantity is not the mixed-model score. The covariance model changes both the score and its information.
Practical Use, Diagnostics, and Alternatives
A score scan fits the null model and evaluates the variant score without fitting a separate alternative model for every variant. This is useful when the null covariance is expensive to estimate but can be reused across many variants. It does not eliminate the need for a sound null fit: a misspecified covariance, poor convergence, or inappropriate treatment of estimated variance components can undermine the resulting p-values.
For a single dosage vector, \(I\) must be positive. If \(I=0\), the dosage is indistinguishable from the null intercept under the covariance-weighted projection, so this score test cannot identify a separate effect. Very small \(I\) can also make the standardized statistic numerically unstable. In software or hand calculations, useful checks include confirming \(P\mathbf{1}=0\), checking that \(I\ge0\) up to numerical tolerance, and verifying that \(U=g^{\mathsf T}Py\) agrees with the expanded weighted-sum formula.
The score test, Wald test, and likelihood-ratio test answer the same null-versus-variant question through different calculations. The score evaluates the likelihood slope at the null; the Wald test uses the fitted effect and its standard error under the alternative; the likelihood-ratio test compares maximized likelihoods. With a known Gaussian covariance and a single linear effect, their results can coincide. In practical mixed-model scans, covariance estimation, approximations, or model details can make their computed statistics differ. The score test's central computational advantage is that it does not require fitting the variant effect under the alternative for each marker.
A common mistake is to use \(g^{\mathsf T}Wg\) as the information even though the intercept was estimated. That omits the intercept adjustment and usually overstates information. The correction is to use the Schur complement \(g^{\mathsf T}Pg\), with the matching score \(g^{\mathsf T}Py\). Another mistake is to regard a plug-in p-value as exact merely because the Gaussian known-\(V\) derivation is exact: once covariance parameters have been estimated, the stated exact result no longer applies automatically.
Check Your Understanding
Work each problem, then compare with the solutions.
- For \(W=\operatorname{diag}(1,1)\), \(g=(0,1)^{\mathsf T}\), and \(y=(1,3)^{\mathsf T}\), calculate \(U\), \(I\), and \(Q\) with an intercept in the null model.
- Why is \(g^{\mathsf T}Wg\) not generally the efficient information when the null intercept is estimated?
- State the null distribution of \(Z\) when \(V\) is known and the Gaussian model assumptions hold. What changes when variance components are estimated and substituted?
- If \(I=0\), what does that mean for the variant test?
Solutions. (1) Here \(\mathbf{1}^{\mathsf T}W\mathbf{1}=2\), \(\mathbf{1}^{\mathsf T}Wy=4\), \(g^{\mathsf T}W\mathbf{1}=1\), and \(g^{\mathsf T}Wy=3\). Thus \(U=3-(1)(4)/2=1\). Also \(g^{\mathsf T}Wg=1\), so \(I=1-1/2=0.5\), and \(Q=1^2/0.5=2\). (2) That expression treats the intercept as known; estimating it requires subtracting the intercept component, giving \(I=g^{\mathsf T}Wg-(g^{\mathsf T}W\mathbf{1})^2/(\mathbf{1}^{\mathsf T}W\mathbf{1})\). (3) \(Z\sim N(0,1)\) exactly under the stated known-\(V\) Gaussian null. With estimated variance components, the reference distribution is generally approximate. (4) The dosage has no remaining variation after covariance-weighted adjustment for the intercept, so the test has no information to distinguish its effect.