GLS Derivation for Correlated Phenotypes

Why Correlated Phenotypes Change the Regression Calculation

In a quantitative-trait genome-wide association study, individuals can be related, so their phenotypes may share genetic contributions. A regression that treats every phenotype as independent can then assign the wrong relative weight to observations and misstate the uncertainty in estimated effects. Generalized least squares (GLS) addresses this by using the phenotype covariance matrix in the regression calculation.

The starting point is the covariance decomposition in a GWAS linear mixed model: a random genetic contribution has covariance determined by the genomic relationship matrix, and an independent residual contribution adds residual variance. Here the focus is what follows from that covariance: how to derive the GLS coefficient estimator and its variance, conditional on the covariance being specified.

Model and Assumptions

Here \(K\) is the genomic relationship matrix, \(I\) is the \(n \times n\) identity matrix, and \(\sigma_g^2\) and \(\sigma_e^2\) are the genetic and residual variance components. For the GLS formula, treat \(V\) as known and positive definite. A positive-semidefinite \(K\) and \(\sigma_e^2>0\) are sufficient for positive definiteness in this model. Also require that \(X\) have full column rank, so its \(q\) columns are linearly independent.

The fixed-effects columns may include an intercept and a genotype dosage column. For example, coding the number of effect alleles as 0, 1, or 2 gives a coefficient in phenotype units per additional effect allele, conditional on the other columns in \(X\). The GLS derivation itself does not require the errors to be normally distributed: the mean and covariance assumptions suffice for unbiasedness and the covariance formula. If \(y\) is multivariate normal, then the GLS estimator is normally distributed as well.

Deriving the GLS Estimator

When errors are independent with equal variance, ordinary least squares minimizes the sum of squared residuals. With correlated errors, residuals in directions of high variance should count less than residuals in directions of low variance. GLS therefore minimizes the residual quadratic form weighted by \(V^{-1}\).

Because \(V\) is positive definite, \(V^{-1}\) exists and is positive definite. Differentiate \(Q(b)\) with respect to \(b\), set the derivative equal to zero, and solve the resulting normal equations:

Full column rank of \(X\) makes \(X^\mathsf{T}V^{-1}X\) positive definite and invertible. Thus the unique minimizer is:

This is a linear transformation of the phenotype vector. Its expectation and covariance follow directly from the model assumptions.

For the covariance result, write \(A=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\), so \(\hat{\beta}_{GLS}=Ay\). Then \(AVA^\mathsf{T}=(X^\mathsf{T}V^{-1}X)^{-1}\), using symmetry of \(V\). A coefficient contrast \(c^\mathsf{T}\hat{\beta}_{GLS}\) consequently has variance \(c^\mathsf{T}(X^\mathsf{T}V^{-1}X)^{-1}c\). The square root of a coefficient’s variance is its standard error, in the same units as that coefficient.

GLS in a Genotype Regression

For a genotype association model, \(X\) can contain an intercept and genotype dosage. The following calculation shows how a non-diagonal covariance affects both coefficients and their covariance. As in the preceding example, all data and covariance values here are hypothetical.

Why GLS Is Efficient, and What Its Assumptions Buy

Among linear unbiased estimators, GLS has the smallest covariance in the positive-semidefinite ordering when \(V\) is the true covariance. This is the generalized least-squares form of the Gauss–Markov result. To see the justification, let \(B=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}\) be the GLS coefficient matrix, and let \(A y\) be any other linear unbiased estimator. Unbiasedness for every \(\beta\) requires \(AX=I_q\), and \(BX=I_q\). Set \(D=A-B\); then \(DX=0\). Expanding the covariance difference gives:

The cross terms vanish because \(VB^\mathsf{T}=X(X^\mathsf{T}V^{-1}X)^{-1}\), so \(DVB^\mathsf{T}=DX(X^\mathsf{T}V^{-1}X)^{-1}=0\); the transpose term is zero as well. Finally, \(DVD^\mathsf{T}\) is positive semidefinite because \(V\) is positive definite. This establishes the covariance comparison without requiring normality.

When \(V=\sigma^2I\), the GLS formula reduces to ordinary least squares, \((X^\mathsf{T}X)^{-1}X^\mathsf{T}y\). If errors are correlated, ordinary least squares can still be unbiased under the stated mean model, but its usual independent-error covariance formula is generally wrong and it does not exploit the covariance structure. Diagonal weighted least squares is appropriate when observations are uncorrelated but have different variances; it does not account for off-diagonal covariance. GLS handles both unequal variances and correlations.

A Mixed-Model Covariance Calculation

Diagnostics, Limitations, and a Common Pitfall

GLS is only as appropriate as its mean and covariance assumptions. In a mixed-model GWAS, the covariance depends on the relationship matrix and the variance components. If the specified \(V\) poorly describes the dependence or variability in the phenotypes, the GLS weights and the resulting coefficient covariance can be unreliable. Residual diagnostics from quantitative-trait GWAS can help identify patterns inconsistent with the fitted model; a covariance matrix should also be checked for symmetry and positive definiteness before it is used.

The formulas above condition on \(V\) being known. In practice, variance components are estimated, so replacing \(V\) by an estimate gives a feasible GLS calculation. The displayed covariance formula then treats that estimate as fixed; it does not, by itself, account for uncertainty from estimating the variance components. Estimating those components is a separate part of mixed-model fitting.

• Rank failure: If a column of \(X\) is redundant, then \(X^\mathsf{T}V^{-1}X\) is singular and the stated unique estimator does not exist. Remove or reparameterize redundant columns.
• Invalid covariance: A non-positive-definite \(V\) cannot be used in this GLS formula. Check the covariance construction and parameters rather than attempting to interpret an invalid inverse.
• Numerical instability: Explicitly forming \(V^{-1}\) is often less stable than solving linear systems using a factorization of \(V\). The mathematical estimator is unchanged by using stable linear algebra.
• Overstated certainty: A standard error computed with estimated variance components is not automatically a complete accounting of all uncertainty. State whether the covariance was treated as known or estimated.

Common mistake: “The phenotype correlation is positive, so GLS must increase the standard error of every coefficient.” That conclusion does not follow. GLS changes the weighting and the covariance of all fitted coefficients; the effect on a particular coefficient’s variance depends on \(X\) and the full matrix \(V\). Compare the relevant diagonal entry of \((X^\mathsf{T}V^{-1}X)^{-1}\), rather than inferring its direction from one correlation alone.