Fixed Covariates in Mixed-Model GWAS

Why Include Fixed Covariates?

A variant can appear associated with a phenotype because its dosage is related to another phenotype predictor, such as age, sex, ancestry, or study batch. A fixed covariate represents a measured predictor whose effect is estimated directly in the model. It is distinct from the random genetic effect, which represents covariance among individuals through a relationship matrix. Including both lets the model account for measured mean differences and residual correlation at the same time.

In a mixed-model GWAS, the coefficient for a tested variant is therefore conditional on the covariates in the model. It describes the modeled change in phenotype per dosage unit after accounting for those covariates and the covariance structure. It need not equal the coefficient from a regression that omits them.

Definition: For one tested variant, a linear mixed model with fixed covariates is $$ \mathbf y=\mathbf C\boldsymbol\gamma+\mathbf g_j\beta_j+\mathbf u+\boldsymbol\varepsilon, \qquad \mathbf u\sim N(\mathbf 0,\sigma_u^2\mathbf K), \qquad \boldsymbol\varepsilon\sim N(\mathbf 0,\sigma_e^2\mathbf I). $$ Here $\mathbf y$ is the $n$-vector of phenotypes, $\mathbf C$ is the fixed-effects design matrix (including an intercept), $\boldsymbol\gamma$ is its coefficient vector, $\mathbf g_j$ is the $n$-vector of dosages at variant $j$, and $\beta_j$ is the variant coefficient. The random genetic effect $\mathbf u$ and residual $\boldsymbol\varepsilon$ are independent. Thus $\operatorname{Var}(\mathbf y)=\mathbf V=\sigma_u^2\mathbf K+\sigma_e^2\mathbf I$.

Columns of $\mathbf C$ might encode age, sex, ancestry principal components, or batch indicators. A categorical variable with $k$ levels is typically represented by $k-1$ indicator columns when the model also has an intercept; the omitted level is the reference category. The entries of $\mathbf g_j$ commonly encode allele dosage as 0, 1, or 2, so $\beta_j$ is in phenotype units per additional copy of the coded allele.

The model assumes that the conditional mean has the specified linear form, that $\mathbf K$ is a valid positive-semidefinite covariance matrix, and that the random effects have the stated mean-zero covariance and independence properties. For the usual Gaussian likelihood-based inference, the random effects and residuals are also taken to be Gaussian. The fixed-effects design must have full column rank, and the variant must not be perfectly explained by the covariates. The phenotype, covariate, genotype, and rows and columns of $\mathbf K$ must refer to the same individuals in the same order.

How Covariates Enter the GLS Estimate

With $\mathbf V$ treated as known, generalized least squares (GLS) estimates fixed effects by minimizing the residual quadratic form weighted by $\mathbf V^{-1}$. For a candidate value $b$ of the variant coefficient, the covariate coefficients minimize $(\mathbf y-\mathbf C\boldsymbol\gamma-\mathbf g_jb)^\mathsf{T}\mathbf V^{-1}(\mathbf y-\mathbf C\boldsymbol\gamma-\mathbf g_jb)$. The resulting variant estimate is the part of the phenotype variation aligned with the variant after accounting for the covariates in this same weighted geometry.

Define the GLS covariate-removal matrix.
Let $\mathbf P_C=\mathbf V^{-1}-\mathbf V^{-1}\mathbf C(\mathbf C^\mathsf{T}\mathbf V^{-1}\mathbf C)^{-1}\mathbf C^\mathsf{T}\mathbf V^{-1}$. It removes the fitted contribution of $\mathbf C$ using the covariance-weighted inner product.

Profile out the covariate coefficients.
For a fixed $b$, the GLS estimate of $\boldsymbol\gamma$ is $(\mathbf C^\mathsf{T}\mathbf V^{-1}\mathbf C)^{-1}\mathbf C^\mathsf{T}\mathbf V^{-1}(\mathbf y-\mathbf g_jb)$. Substituting this into the quadratic form leaves $(\mathbf y-\mathbf g_jb)^\mathsf{T}\mathbf P_C(\mathbf y-\mathbf g_jb)$.

Minimize over the variant coefficient.
Differentiate the profiled quadratic form with respect to $b$, set the derivative to zero, and solve. The denominator is positive when the variant is not wholly explained by the covariates.

\frac{\partial}{\partial b}\left[(\mathbf y-\mathbf g_jb)^\mathsf{T}\mathbf P_C(\mathbf y-\mathbf g_jb)\right] =-2\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf y +2b\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf g_j=0, \qquad \hat\beta_j= \frac{\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf y} {\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf g_j}.

This is the GLS version of partial regression: both phenotype and genotype are compared after the covariates have been accounted for. If $\mathbf V$ is known, the variance of this estimate under the model is $(\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf g_j)^{-1}$. In practice, variance components are usually estimated, so standard errors using the fitted covariance are plug-in quantities. The formula also shows why a variant that is nearly predictable from the covariates contains little independent information: its residualized contribution, and hence the denominator, is small.

Worked Example: A Covariate Changes the Variant Estimate

Consider four hypothetical individuals with phenotype $y$, variant dosage $g$, and a binary covariate $c$:

Individual	$y$	$g$	$c$
1	1	0	0
2	4	0	1
3	6	1	1
4	6	1	1

For clarity, suppose the covariance used in this calculation is $\mathbf V=\mathbf I$. First fit the unadjusted regression with an intercept and $g$. The mean phenotype among $g=0$ individuals is $(1+4)/2=2.5$; among $g=1$ individuals it is $(6+6)/2=6$. Thus the unadjusted dosage coefficient is $6-2.5=3.5$ phenotype units per dosage unit.

Now fit the model with an intercept, $g$, and $c$. The data satisfy $y=1+2g+3c$ for all four individuals. For example, individual 2 has fitted value $1+2(0)+3(1)=4$, and individual 3 has fitted value $1+2(1)+3(1)=6$. The adjusted estimate is therefore $\hat\beta=2$, with covariate coefficient 3. The genotype coefficient changes from 3.5 to 2 because the $g=1$ individuals also have a higher prevalence of $c=1$. In this hypothetical data set, the adjusted model fits exactly; that feature is for arithmetic clarity, not a realistic claim about genomic data.

Covariate Adjustment with a Non-Identity Covariance

In a mixed model, covariate adjustment is not generally ordinary mean-centering or unweighted regression. It uses $\mathbf V^{-1}$, so observations with different covariance patterns contribute differently. The following special case shows how a random-effect covariance can be present even when the coefficient calculation has a simple result.

Worked Example: A Random Intercept Covariance and an Fixed Intercept

Use the four observations in the preceding example and suppose $\mathbf K=\mathbf 1\mathbf 1^\mathsf{T}$, $\sigma_u^2=0.5$, and $\sigma_e^2=1$. This is a valid, though deliberately simple, covariance structure: the random effect is shared by all individuals. Then $\mathbf V=\mathbf I+0.5\mathbf 1\mathbf 1^\mathsf{T}$, and direct multiplication gives $\mathbf V^{-1}=\mathbf I-\frac{1}{6}\mathbf 1\mathbf 1^\mathsf{T}$.

Because $\mathbf C$ contains an intercept, GLS contrasts that distinguish the variant and covariate effects are orthogonal to the shared $\mathbf 1$ direction. On that contrast space, $\mathbf V^{-1}$ acts as the identity: if $\mathbf 1^\mathsf{T}\mathbf v=0$, then $\mathbf V^{-1}\mathbf v=\mathbf v-\frac{1}{6}\mathbf 1(\mathbf 1^\mathsf{T}\mathbf v)=\mathbf v$. Consequently this particular covariance does not change the fixed-effect estimates from the preceding example: the adjusted variant estimate remains 2. This is a special property of a shared random intercept and an included fixed intercept, not a general rule that covariance modeling leaves coefficients unchanged.

For a general relationship matrix, $\mathbf V^{-1}$ can weight different directions in the data differently. Then the adjusted estimate can differ from the ordinary least-squares estimate even when the same fixed covariates are included. The essential distinction is that $\mathbf C$ models systematic mean differences, whereas $\mathbf K$ contributes to the covariance of observations around that mean.

Interpreting the Conditional Effect

The adjusted coefficient answers a specific model-based question: among individuals compared at the same modeled covariate values, what phenotype difference is associated with a one-unit difference in dosage, accounting for the fitted covariance? For a quantitative phenotype measured in millimetres, for example, $\hat\beta_j=0.4$ means an estimated 0.4 mm higher phenotype per additional copy of the coded allele, conditional on the included covariates and covariance model. It is not automatically a causal effect.

Covariate choice therefore matters. An omitted variable related to both dosage and phenotype may leave confounding. Conversely, adding an irrelevant or poorly measured predictor can reduce precision, and adding highly collinear predictors can make the adjusted variant estimate unstable. A covariate should have a defensible role in the analysis rather than being selected only because it changes a result.

Worked Example: Adjusting for Study Batch

Suppose a hypothetical study records phenotype $y$, dosage $g$, and batch for five individuals. Batch A has observations $(g,y)=(0,10),(1,12)$. Batch B has $(g,y)=(0,20),(1,22),(1,24)$. Fit a linear model with an intercept, a batch indicator, and dosage, with a common dosage coefficient for both batches.

The unadjusted coefficient is the difference in mean phenotype between dosage groups. For $g=0$, the mean is $(10+20)/2=15$. For $g=1$, the mean is $(12+22+24)/3=58/3=19.333$. The unadjusted coefficient is $58/3-15=13/3=4.333$ phenotype units per dosage unit.

With batch in the model, the common dosage coefficient is the within-batch association pooled across batches. For batch A, the dosage slope is $12-10=2$, and the within-batch sum of squared centered dosages is $0.5$. For batch B, the dosage-1 mean is $(22+24)/2=23$, so its slope relative to dosage 0 is $23-20=3$; its within-batch sum of squared centered dosages is $2/3$. The pooled adjusted slope is therefore $(2(0.5)+3(2/3))/(0.5+2/3)=3/(7/6)=18/7=2.571$.

The adjusted estimate is 2.571 phenotype units per dosage unit, rather than 4.333. The difference arises because the unadjusted comparison mixes within-batch dosage differences with the higher baseline phenotype in batch B. This example assumes the dosage effect is common across batches; a model with a dosage-by-batch interaction would ask a different question.

Diagnostics, Alternatives, and Limitations

Before interpreting an adjusted association, check that the design is estimable and represents the intended comparison. Useful checks include:

Confirm that each covariate is coded as intended and that categorical variables have an appropriate reference level.
Check for rank deficiency or near-collinearity among the variant dosage and columns of $\mathbf C$. A nearly zero value of $\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf g_j$ indicates little genotype information beyond the covariates.
Verify that phenotype, genotype, covariates, and relationship matrix use matching individuals and ordering. Misalignment can create a plausible-looking but invalid fit.
Inspect missingness and sample inclusion. If the individuals analyzed differ across variants, estimates may not be directly comparable across the scan.
Assess whether the linear mean model and covariance specification are scientifically plausible. Fixed covariate adjustment does not repair an incorrectly specified relationship matrix or guarantee removal of unmeasured confounding.

The main alternative is an unadjusted regression or association test, which omits $\mathbf C$. It may be simpler, but it targets a marginal association and can mix covariate differences with the variant association. Including the covariates in the mixed model targets the conditional coefficient described above. A further alternative is to change the covariance model; that addresses dependence and weighting, not the choice of fixed mean predictors. These approaches are not interchangeable, and the appropriate specification depends on the study design and scientific question.

Common pitfall: “The variant coefficient got smaller after adding ancestry covariates, so the variant has no effect.” This conclusion does not follow. Adjustment changes the estimand from a marginal association to a conditional one; the change may reflect confounding, collinearity, or reduced information. Examine the adjusted estimate, its uncertainty, the covariate relationships, and the model assumptions before drawing a conclusion.

Key takeaway: Fixed covariates enter the mean of the mixed model, while the random genetic effect enters its covariance. The tested variant coefficient is estimated from the genotype information remaining after covariance-weighted adjustment for the fixed covariates.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

In the model, what is the distinction between $\mathbf C\boldsymbol\gamma$ and $\mathbf u$?
Why can the coefficient from a regression omitting ancestry covariates differ from the coefficient after their inclusion?
What does a small value of $\mathbf g_j^\mathsf{T}\mathbf P_C\mathbf g_j$ indicate?
For the batch example, state the adjusted dosage coefficient and its units in words.

Solutions. (1) $\mathbf C\boldsymbol\gamma$ is the fixed, modeled mean contribution of measured covariates; $\mathbf u$ is a mean-zero random effect whose covariance is specified through $\mathbf K$. (2) The unadjusted estimate can combine within-covariate dosage differences with phenotype differences associated with those covariates; adjustment estimates the conditional association instead. (3) It indicates little genotype variation remains after covariate adjustment in the GLS metric, so the variant coefficient is weakly informed and may be unstable. (4) The estimate is $18/7=2.571$ phenotype units per dosage unit, conditional on batch under the assumed common dosage effect.

Individual	\(y\)	\(g\)	\(c\)
1	1	0	0
2	4	0	1
3	6	1	1
4	6	1	1