Mixed Models for Population Structure Control

Why Model Genetic Similarity?

A genetic association analysis can be distorted when individuals with similar ancestry also tend to share phenotype-related genetic or environmental factors. Their outcomes may therefore be correlated, while their allele frequencies may differ from those of other ancestry groups. Treating all observations as independent can make the uncertainty in a SNP association estimate too small and can produce misleading evidence of association.

A linear mixed model addresses this dependence by including a random effect whose covariance across individuals is specified using a genomic relationship matrix. The tested SNP remains a fixed effect: its coefficient describes the estimated phenotype difference associated with one additional effect allele, conditional on the other terms in the model. The relationship matrix instead describes how much background similarity is expected between each pair of individuals.

Model and Interpretation

For \(n\) individuals, let \(y\) be the vector of continuous phenotypes, \(g\) the vector of dosages for the SNP being tested, and \(C\) a matrix of measured covariates such as age or sex. Dosage is commonly coded as the number of copies of a specified allele, from 0 to 2. Write the fixed-effect design matrix as \(X=[\mathbf{1}, C, g]\), with corresponding coefficient vector \(\beta\). The model is

Here \(u\) is the random background genetic effect, \(e\) is residual error, and \(I\) is the \(n \times n\) identity matrix. The standard model assumes \(u\) and \(e\) are independent, both have mean zero, and the residual errors have equal variance and are mutually independent. The relationship matrix \(K\) is treated as known for the association fit and must be positive semidefinite. Its scale depends on how it was constructed; often its diagonal entries are near 1, but this is not guaranteed.

The fixed effects describe the conditional mean, while the random effect describes covariance around that mean. Specifically, the model implies

For different individuals \(i\) and \(j\), the modeled covariance is \(V_{ij}=\sigma_g^2 K_{ij}\). A larger \(K_{ij}\) therefore implies greater covariance in their background effects, when \(\sigma_g^2\) is positive. The diagonal elements of \(V\) give the modeled phenotype variances. The tested SNP's coefficient, denoted \(\beta_g\), is interpreted in the phenotype's units per one additional copy of the effect allele, conditional on the included covariates and the modeled background structure. It is an association estimate, not by itself evidence of a causal effect.

Estimating and Testing the SNP Effect

The variance components \(\sigma_g^2\) and \(\sigma_e^2\) are generally estimated from the data, commonly by maximum likelihood or restricted maximum likelihood. Once estimates are available, they determine an estimated covariance matrix \(\hat V\). The fixed effects can then be estimated by generalized least squares (GLS). GLS gives less independent weight to observations that are modeled as strongly correlated, rather than treating every phenotype as an independent contribution.

For a specified positive-definite covariance matrix \(V\), the GLS estimate and its covariance are

In practice, variance components are estimated, so association software substitutes \(\hat V\) and typically reports an estimated standard error. To test the SNP association, the hypotheses are \(H_0:\beta_g=0\) and \(H_a:\beta_g\ne 0\). A common large-sample Wald statistic is \(z=\hat\beta_g/SE(\hat\beta_g)\); under the null and the model assumptions, it is approximately standard normal. A two-sided p-value is \(2P(Z\ge |z|)\) for \(Z\sim N(0,1)\). A corresponding approximate 95% confidence interval is \(\hat\beta_g \pm 1.96SE(\hat\beta_g)\).

This large-sample calibration is not an automatic guarantee of valid inference. It relies on an adequate mean and covariance model, suitable estimation of variance components, and sufficient information for the approximation being used. The model is normally used for quantitative phenotypes; binary phenotypes require a model appropriate to a binary outcome rather than simply assuming a normally distributed continuous response.

Why the Covariance Term Helps With Structure

Ancestry can create genetic similarity across many markers, not just at the SNP under test. When this background similarity is also associated with phenotype, the random effect \(u\) gives the model a way to represent correlated phenotype deviations among genetically similar individuals. The SNP coefficient is then estimated using the covariance-weighted information in the sample, rather than using the independent-observation assumption.

This adjustment is not equivalent to asserting that ancestry has no relationship with phenotype, nor does it guarantee that every source of population confounding has been removed. It relies on \(K\) representing relevant background similarity and on the random-effect covariance being an adequate description of how that similarity relates to phenotype. A poorly constructed or poorly scaled matrix, an inadequate mean model, or a mismatch between the assumed and actual covariance can leave residual confounding or give unreliable standard errors. The relationship-matrix preparation and diagnostics matter as much as the choice to use a mixed model.

Deriving the GLS Estimate

The GLS formula follows by treating the phenotype vector as multivariate normal with mean \(X\beta\) and covariance \(V\), for fixed \(V\). Up to a constant, its log-likelihood as a function of \(\beta\) is

Differentiating with respect to \(\beta\) and setting the derivative to zero gives

If \(X^\mathsf{T}V^{-1}X\) is invertible, solving these normal equations yields the GLS estimate stated above. Its covariance is \((X^\mathsf{T}V^{-1}X)^{-1}\), conditional on \(X\) and known \(V\); in applications, replacing \(V\) with an estimate makes the reported uncertainty an approximation that also depends on variance-component estimation.

Common Pitfall and Practical Interpretation

• Wrong reasoning: “The model includes a relationship matrix, so ancestry confounding is automatically eliminated.” Correction: The matrix specifies a covariance pattern. Control depends on whether that pattern and the fixed-effect mean model adequately represent the relevant structure.
• Wrong reasoning: “A large random-effect variance means the SNP has a large effect.” Correction: \(\sigma_g^2\) describes variation in the background random effect under the chosen \(K\) scale; \(\beta_g\) is the tested SNP's fixed-effect association coefficient. They are different quantities.
• Wrong reasoning: “The random effect is another covariate, so its coefficient is an ancestry effect.” Correction: The random effect has a covariance model across individuals; it is not a single fixed ancestry coefficient with a direct per-unit interpretation.
• Wrong reasoning: “A significant Wald test proves that the tested allele changes the phenotype.” Correction: The test assesses evidence for a nonzero conditional association under the specified model. It does not by itself establish a causal effect.