Heteroskedasticity And Standard Errors In Linear GWAS

Why Residual Variance Matters in a Linear GWAS

In a quantitative-trait GWAS, a linear regression estimates how a trait changes with SNP genotype, usually while adjusting for covariates. The coefficient describes the estimated change in the trait per unit of the genotype coding; its standard error describes the estimated sampling variability of that coefficient. If individuals with different genotypes or covariate values have different residual variability, the coefficient can still be useful, but the usual standard-error formula may no longer be valid.

Here the focus is the uncertainty calculation, not a change to the mean model. In particular, using a heteroskedasticity-robust standard error does not remove confounding, correct an incorrectly specified conditional mean, or establish that a SNP causes a trait difference. It changes how uncertainty is estimated when the residual variance is not constant.

The key assumptions for the standard robust calculation are that the conditional mean is correctly specified, \(E(\varepsilon_i\mid G_i,C_i)=0\), and the observations are independent across individuals. The conditional variance may differ between individuals, but should be finite; for the usual large-sample inference, no small set of observations should dominate the fit. Genotype coding and covariate adjustment must also match the scientific model. These assumptions are distinct: allowing unequal variances does not protect against a nonzero conditional mean error or dependence between observations.

With an additive genotype dosage coded 0, 1, or 2, \(\beta_G\) is the fitted trait difference per additional effect allele, conditional on the included covariates. For example, if the trait is measured in milligrams per deciliter, \(\beta_G\) is in milligrams per deciliter per allele. This interpretation follows the quantitative-trait effects covered in the earlier material on interpreting quantitative-trait effects.

Conventional and Robust Standard Errors

Write the regression in matrix form as \(Y=X\beta+\varepsilon\), where \(X\) includes an intercept, the SNP dosage, and the covariates. Ordinary least squares estimates \(\hat\beta=(X^\mathsf{T}X)^{-1}X^\mathsf{T}Y\), provided the columns of \(X\) are linearly independent. If the conditional mean is correct, then \(E(\hat\beta\mid X)=\beta\). Under homoskedasticity and independent errors, the conventional covariance estimator is \(s^2(X^\mathsf{T}X)^{-1}\), where \(s^2\) is the residual sum of squares divided by \(n-p\), and \(p\) is the number of fitted regression coefficients.

When residual variance differs by individual, the actual conditional covariance of the estimator is instead \(\operatorname{Var}(\hat\beta\mid X)=(X^\mathsf{T}X)^{-1}X^\mathsf{T}\Omega X(X^\mathsf{T}X)^{-1}\), where \(\Omega\) is diagonal, with entry \(i\) equal to \(\operatorname{Var}(\varepsilon_i\mid X)\). The conventional formula replaces all these different variances with one pooled estimate. Its standard errors can therefore be too small or too large.

A heteroskedasticity-robust sandwich estimator replaces the unknown errors and their variances with fitted residuals. One widely used finite-sample adjustment is HC1. Let \(\hat\varepsilon_i=Y_i-X_i^\mathsf{T}\hat\beta\), where \(X_i^\mathsf{T}\) is row \(i\) of the design matrix. The HC1 covariance estimator is:

The standard error for the SNP coefficient is the square root of the diagonal entry corresponding to \(\hat\beta_G\). HC1 multiplies the HC0 sandwich estimate by \(n/(n-p)\). This correction does not make the estimator exact in small samples; it is a common adjustment to the asymptotic sandwich estimator. Other robust variants, such as HC2 and HC3, adjust residual contributions differently, particularly to account for leverage.

Why the Sandwich Formula Works

The robust calculation follows from the way least squares combines individual outcomes. Conditional on \(X\), the estimation error is \(\hat\beta-\beta=(X^\mathsf{T}X)^{-1}X^\mathsf{T}\varepsilon\). If errors are independent but have variances \(\sigma_i^2\), their covariance matrix is diagonal, \(\Omega=\operatorname{diag}(\sigma_1^2,\ldots,\sigma_n^2)\). Applying the covariance rule for a linear transformation gives the true covariance below. Replacing each unknown \(\sigma_i^2\) by \(\hat\varepsilon_i^2\) yields HC0; HC1 then applies its degrees-of-freedom multiplier.

This derivation also clarifies what the robust estimator does not do. It estimates uncertainty for the specified least-squares coefficients under a potentially nonconstant diagonal error covariance. It does not alter the fitted coefficients, and it does not accommodate correlated errors between individuals merely by using HC1.

Worked Calculations

Diagnostics, Alternatives, and Limitations

Residual diagnostics can reveal patterns that make a constant-variance assumption implausible. Plot residuals against fitted values, genotype dosage, and important covariates; a changing vertical spread or a funnel shape suggests that residual variance may depend on the fitted mean or predictors. In a GWAS, inspect plots and summaries rather than relying only on a formal test: with large samples, small departures can be statistically detectable without materially changing a standard error. A residual plot cannot prove homoskedasticity, and the fitted residuals are estimates, not the unobserved errors themselves.

The conventional formula is more efficient when its homoskedasticity assumptions hold, but can misestimate uncertainty otherwise. HC1 is a practical alternative under independent observations and a correctly specified conditional mean; with a large enough sample, it estimates the sandwich covariance without requiring equal variances. If observations are dependent, such as through family or other grouping structure, the individual-level HC1 formula is not a general solution. Correlation requires an approach that accounts for that dependence, rather than simply substituting HC1.

Robust standard errors also have limitations. They do not repair omitted-variable bias, an unsuitable genotype model, influential observations, or a wrong conditional mean. They can be unstable in small samples or when a few individuals have high leverage. Very large residuals can strongly affect the squared-residual terms, while HC1's simple multiplier does not specifically correct high leverage. Check sample size, leverage, residual patterns, genotype counts, and model specification; consider whether a different variance model or a dependence-aware analysis is scientifically warranted.

A frequent mistake is to assume that “robust” means “larger” or “safer in every respect.” The examples show that HC1 can be larger or smaller than the conventional standard error, depending on how residual magnitudes align with the design. The correction is to report which estimator was used, inspect whether its assumptions are reasonable, and interpret it as an uncertainty estimate for the fitted model—not as a remedy for every possible GWAS problem.