Sandwich Standard Errors for Genetic Association

Why Use Sandwich Standard Errors?

A genetic association regression estimates the relationship between a phenotype and genotype, often while adjusting for covariates. The fitted coefficient describes the association under the specified regression model. Its standard error, however, also depends on assumptions about how outcomes vary around that model. For a quantitative trait, the usual linear-regression standard error assumes a common residual variance. If residual variability instead changes with genotype or covariates, that standard error can be inaccurate even when the conditional mean is modeled correctly.

A sandwich standard error estimates uncertainty using two ingredients: how responsive the estimating equation is to changes in the coefficients, and how variable the individual observations’ contributions to that equation are. This construction applies to both linear and logistic association models. Its validity still depends on important conditions, particularly independent observations and a suitable mean model.

Definition and General Form

The two matrices have distinct roles. The bread describes how much the estimating equation changes when the coefficient changes; the meat records the variability of the individual contributions. In ordinary regression examples, the bread is symmetric, so the expression is often written with the bread on both sides without showing a transpose. The general formula includes the transpose: a covariance matrix must be symmetric, and an arbitrary sensitivity matrix need not be symmetric.

In finite-sample notation, let \(\hat{A}=-\sum_i \partial\psi_i(\hat{\beta})/\partial\beta^\mathsf{T}\) and \(\hat{B}=\sum_i\psi_i(\hat{\beta})\psi_i(\hat{\beta})^\mathsf{T}\). The empirical covariance estimate is \(\hat{A}^{-1}\hat{B}(\hat{A}^{-1})^\mathsf{T}\). This is commonly called the HC0 sandwich estimate in regression settings. It is a large-sample variance estimator, not a guarantee of accurate inference in a small sample.

Linear Association Regression

For individual \(i\), let \(Y_i\) be a quantitative trait and let \(x_i\) be a \(p\)-dimensional column vector containing an intercept, genotype coding, and any included covariates. The linear mean model is \(E(Y_i\mid x_i)=x_i^\mathsf{T}\beta\). The genotype coefficient is interpreted in the units of \(Y_i\) per unit of the chosen genotype coding, conditional on the included covariates.

The OLS estimator solves \(\sum_i x_i(Y_i-x_i^\mathsf{T}\beta)=0\). Its individual contribution is \(\psi_i(\beta)=x_i(Y_i-x_i^\mathsf{T}\beta)\). Differentiating gives \(-x_ix_i^\mathsf{T}\), so the empirical bread is \(X^\mathsf{T}X\), where \(X\) has rows \(x_i^\mathsf{T}\). The meat is \(\sum_i x_ix_i^\mathsf{T}\hat{e}_i^2\), where \(\hat{e}_i=Y_i-x_i^\mathsf{T}\hat{\beta}\). Thus:

Here the bread is symmetric, which is why the final factor can be written without a transpose. The diagonal element for the genotype coefficient estimates its variance; its square root is the sandwich standard error. The formula requires a full-rank design matrix. For consistency, the linear conditional mean must be appropriate, observations must be independent, and the data must have suitable finite moments. The key advantage over the conventional variance is that the error variance need not be constant across \(x_i\).

Under the stronger homoskedastic model, \(\operatorname{Var}(Y_i\mid x_i)=\sigma^2\), the model-based covariance is estimated by \(s^2(X^\mathsf{T}X)^{-1}\), where \(s^2=\sum_i\hat{e}_i^2/(n-p)\). If the conditional variance is constant, the expected meat is proportional to the bread, and the robust and model-based estimates target the same large-sample covariance. If conditional variances differ, the model-based formula generally misses that pattern; the sandwich estimates it through the squared residuals and design vectors.

Logistic Association Regression

For a binary trait \(Y_i\in\{0,1\}\), the logistic mean model specifies \(P(Y_i=1\mid x_i)=p_i\), with \(\operatorname{logit}(p_i)=x_i^\mathsf{T}\beta\). The coefficient for genotype is a conditional log odds ratio per unit of the genotype coding, when the other covariates are held fixed. The individual score contribution is \(\psi_i(\beta)=x_i(Y_i-p_i)\). Its negative derivative is \(p_i(1-p_i)x_ix_i^\mathsf{T}\), giving bread \(\sum_i p_i(1-p_i)x_ix_i^\mathsf{T}\) at the fitted values. The empirical meat is \(\sum_i x_ix_i^\mathsf{T}(Y_i-\hat{p}_i)^2\).

For an independent Bernoulli outcome with a correctly specified conditional probability, the conditional variance is \(p_i(1-p_i)\). Consequently, \(E\{(Y_i-p_i)^2\mid x_i\}=p_i(1-p_i)\), so the expected meat equals the expected bread. This explains why the robust and usual model-based logistic covariance estimates are asymptotically equivalent under the correctly specified independent Bernoulli model. The robust form is useful when the variance or likelihood assumptions used for model-based inference are not trusted, but it does not make a wrong conditional mean correct. If that mean is misspecified, the coefficient may target a best-fitting, model-dependent value rather than the intended association parameter.

Why the Transpose Matters

The sandwich form is not merely a convenient way to write a variance. Its matrix order follows from propagating the variability of the estimating equation through the coefficient estimator. Linearizing the estimating equation around its target gives a coefficient error approximately equal to the inverse sensitivity times the average contribution. Taking the covariance of that expression places the transpose on the second inverse:

For example, suppose an estimating system has \(A=\begin{pmatrix}1&1\\0&1\end{pmatrix}\) and \(B=I\), the two-dimensional identity matrix. Then \(A^{-1}=\begin{pmatrix}1&-1\\0&1\end{pmatrix}\), and the covariance calculation is \(A^{-1}B(A^{-1})^\mathsf{T}=\begin{pmatrix}2&-1\\-1&1\end{pmatrix}\). It is symmetric, as a covariance matrix must be. Omitting the transpose would give \(\begin{pmatrix}1&-2\\0&1\end{pmatrix}\), which is not symmetric and cannot be a covariance matrix. In the OLS and logistic examples above, the bread is symmetric, so the transpose does not change the displayed regression formula.

Assumptions, Limitations, and Common Mistakes

A sandwich standard error concerns uncertainty under the chosen estimating equation; it is not a general repair for an association analysis. In the linear model, the robust covariance protects against nonconstant conditional residual variance when the conditional mean is correctly specified. In logistic regression with independent Bernoulli observations, the model-based and sandwich forms agree asymptotically if the probability model is correct. In either setting, the usual individual-level sandwich assumes independent observations and adequate sample information for large-sample approximation. In particular, it does not account for dependence among relatives; that requires a method that treats the dependence appropriately.

Robust standard errors also do not remove confounding, correct an inappropriate genotype coding, fix selection bias, or resolve bias from missing data. Nor do they compensate for a misspecified mean relationship between phenotype, genotype, and covariates. A robust standard error can quantify uncertainty around the coefficient estimated by the fitted model while that coefficient remains a poor estimate of the scientific target.

• Mistake: “Robust” means robust to any model failure. The correction is that sandwich variance estimates address specified forms of variance misspecification under their sampling assumptions. They do not fix a wrong mean model or dependence between observations.
• Mistake: Dropping the transpose in the general formula. The correction is to use \(A^{-1}B(A^{-1})^\mathsf{T}\). Writing \(A^{-1}BA^{-1}\) is justified in these regression formulas because their bread matrix is symmetric, not as a general rule.
• Mistake: Treating HC0 as a small-sample guarantee. The correction is that HC0 is a large-sample estimator; it can perform poorly with few observations or influential design points. Variants such as HC1 multiply the HC0 covariance by \(n/(n-p)\) in OLS, a finite-sample adjustment that does not make it exact.