Why Mixed-Model GWAS Errors Matter
A linear mixed model (LMM) can account for genetic similarity among sampled individuals, but it does not automatically make a genome-wide association study (GWAS) reliable. The model must represent the sample and its dependence structure well enough for its standard errors and test statistics to be useful. Errors in the relationship matrix, numerical fitting, or interpretation of diagnostics can undermine that goal even when the model formula looks appropriate.
This tutorial focuses on four practical failure modes: omitting relatedness, using an unsuitable relationship matrix, ignoring convergence warnings, and overinterpreting inflation statistics. The earlier material on reading QQ plots after mixed-model adjustment and distinguishing confounding from polygenic signal provides useful context: diagnostics are evidence to evaluate, not a substitute for checking the model and study design.
The Model and What Its Assumptions Buy
For a quantitative phenotype measured on \(n\) individuals, consider testing one SNP at a time with the model
Here \(y\) is the \(n\)-vector of phenotypes, \(C\) is the matrix of included covariates (typically including an intercept), \(x_j\) is the \(n\)-vector of dosages for SNP \(j\), and \(\beta_j\) is its conditional effect in phenotype units per additional dosage copy. The random effect \(u\) represents covariance associated with genetic similarity; \(K\) is an \(n\times n\) relationship matrix. The residuals \(\varepsilon\) are assumed independent with common variance, and independent of \(u\). Consequently,
This model relies on more than having a matrix called \(K\). Its rows and columns must correspond to the analyzed individuals in the same order as \(y\) and \(x_j\); it must be a valid covariance structure (in particular, positive semidefinite); and it should represent the relevant dependence adequately. The conditional mean and residual assumptions must also be reasonable for the trait and analysis. For example, this Gaussian LMM is not automatically an appropriate model for a binary outcome.
Under a correctly specified covariance \(V\), the generalized least-squares estimate of the fixed effects has covariance \((X^\mathsf{T}V^{-1}X)^{-1}\), where \(X\) contains the tested genotype and the covariates. Thus, an error in \(V\) can affect the estimated uncertainty even if the estimated SNP coefficient appears plausible. A Wald test, for example, compares \(\hat{\beta}_j\) with its estimated standard error; its p-value is only as trustworthy as the fitting and calibration assumptions supporting that standard error and reference distribution.
Worked Example: The Cost of Treating Related Observations as Independent
Worked Example: The Cost of Treating Related Observations as Independent
Consider a simplified illustration with two related individuals whose residual-like measurements each have variance \(1\) and correlation \(0.50\). Compare the variance of their average under the stated dependence with the variance calculated by incorrectly assuming independence. These numbers illustrate the mechanism; they are not a general formula for a GWAS SNP-effect standard error.
For two measurements \(Y_1,Y_2\), the average is \(\bar{Y}=(Y_1+Y_2)/2\). Because each variance is \(1\), their covariance is \(0.50\). Therefore,
The corresponding standard error is \(\sqrt{0.75}=0.8660\). Under the false independence assumption, the covariance term is set to zero, giving variance \((1+1)/4=0.50\) and standard error \(\sqrt{0.50}=0.7071\). The independence calculation is too small in this example because the two observations contain overlapping information.
For GWAS, the precise effect on a SNP test depends on how genotype patterns and relatedness align, so this calculation does not establish that every omitted-relatedness analysis will be anti-conservative. It does show why treating correlated observations as independent can miscalculate information and uncertainty. Depending on the data, test statistics may be poorly calibrated and apparent association evidence may be exaggerated or otherwise distorted.
Failure Mode 1: Omitting or Underrepresenting Relatedness
If related individuals are analyzed with an ordinary regression that assumes independent errors, the fitted model omits covariance that is present in the data. The consequences depend on the pattern of relatedness and its alignment with genotype and phenotype. One common risk is that a SNP correlated with family membership or another structured feature can appear associated because the model attributes structured phenotype variation to the SNP. Standard errors and null test statistics can also be miscalibrated.
The correction is not simply “always use an LMM.” Check the sample definition, relatedness structure, and intended analysis. If a mixed model is appropriate, verify that relatedness is represented in \(K\), then assess whether the fitted model and its diagnostics are plausible. A model with a random effect cannot repair every omitted covariate, phenotype problem, or sampling bias.
Failure Mode 2: Using an Unsuitable Relationship Matrix
A relationship matrix can be unsuitable even if its dimensions are correct. Common causes include using a matrix built for a different set of individuals, misordering its rows or columns, relying on a poorly matched genetic reference, or using a matrix that captures a dependence structure different from the one relevant to the analysis. A matrix that is numerically valid can still be scientifically inadequate.
To see why the entries matter, suppose four hypothetical individuals form two pairs. Imagine the relevant covariance structure has variance components \(\sigma_g^2=0.50\) and \(\sigma_e^2=0.50\), with \(K\) having diagonal entries \(1\), within-pair entries \(0.60\), and zero between-pair entries. The model gives marginal variance \(0.50(1)+0.50=1\) per individual and within-pair covariance \(0.50(0.60)=0.30\). If the fitted matrix is instead \(I\), it assigns within-pair covariance \(0\). It has erased a feature of the assumed dependence structure.
That mismatch does not imply a universal direction of error: the effect on a particular SNP’s estimate, standard error, and p-value depends on the data and genotype pattern. It does mean the fitted covariance model no longer reflects the specified structure. Check individual IDs and ordering, matrix dimensions, symmetry and eigenvalue validity, and whether the construction and sample make sense together. Where possible, compare diagnostics and results under defensible alternative specifications rather than selecting a matrix because it produces more discoveries or a more reassuring QQ plot.
Worked Example: Finding a Sample-Alignment Error
Worked Example: Finding a Sample-Alignment Error
Suppose a hypothetical analysis contains phenotypes for individuals A, B, C, and D in that order. The relationship matrix was assembled in the order B, A, C, D, but the software was given no indication of the difference. Its dimensions are still \(4\times4\), so a dimension check alone will not detect the problem.
If the intended relationship between A and B is \(0.60\), the first two rows and columns of the intended matrix should associate the \(0.60\) entry with A and B. In the supplied ordering, the first position refers to B and the second to A. If the software instead interprets position one as A and position two as B without reordering the matrix, other entries can be assigned to the wrong people. The matrix will then describe a different covariance pattern for the phenotype vector.
A concrete check is to align by identifiers before fitting: the phenotype, genotype, covariate, and matrix index for each position must all name the same person. Here, reorder the matrix rows and columns into A, B, C, D order and verify the resulting entries against the intended pair structure. The check concerns labels, not just matrix size. A successful numerical fit would not, by itself, reveal that the identities were mismatched.
Failure Mode 3: Ignoring Convergence Warnings
Variance components are usually estimated numerically. An optimizer warning may indicate that its stopping criterion was not met, a gradient remains substantial, an estimate is on a boundary, or the likelihood surface is poorly determined. A warning does not prove that the resulting association is wrong, but it does mean that the analyst should not silently treat the displayed estimates and p-values as automatically reliable.
Inspect the specific warning and fitting output. Check whether estimates are finite and within allowed bounds, whether the optimizer reached its stated convergence criterion, and whether the result changes materially when fitting is repeated with sensible starting values or a different numerical strategy. Check that the phenotype, covariates, and \(K\) are correctly scaled and aligned. A boundary estimate can be scientifically possible, but it can also signal weak information or model difficulty; its meaning depends on the setting and software.
Worked Example: A Warning Can Change an Association Conclusion
Suppose a hypothetical software run prints a SNP estimate of \(0.12\) phenotype units per dosage copy, a standard error of \(0.04\), and a two-sided normal-approximation Wald p-value. The run also reports that the variance-component optimizer did not converge. The statistic calculated from the printed values is
where \(Z\sim N(0,1)\). The arithmetic is correct, but the warning means the reported standard error and its calibration have not yet been established as trustworthy. For illustration, if an adequately checked refit gave the same coefficient but standard error \(0.06\), then \(z=0.12/0.06=2.00\), with two-sided normal-approximation \(p\approx0.0455\). These hypothetical outputs show why a failed fit cannot be dismissed as cosmetic; they do not imply that refitting must produce these particular values.
The appropriate conclusion from the first run is not “the SNP is significant” or “the SNP is not associated.” It is that the numerical result needs investigation before inference. If a stable, defensible fit cannot be obtained, report the analysis as unresolved rather than choosing whichever run gives the preferred p-value.
Failure Mode 4: Misreading Inflation Statistics
A QQ plot compares observed test statistics with their null expectations, and the genomic inflation factor \(\lambda\) is commonly calculated as the median observed 1-degree-of-freedom chi-square statistic divided by its null median, approximately \(0.455\). Inflation can reflect confounding or model misspecification, but it can also arise from widespread polygenic signal, sample size, or other features. Conversely, a near-null median does not establish that every locus is calibrated.
Worked Example: Interpreting a Genomic Inflation Factor
Suppose a hypothetical GWAS has median observed chi-square statistic \(0.620\). Using the null median \(0.455\), the estimated genomic inflation factor is
This value says that the median observed statistic is about \(36.3\%\) larger than the reference median under this calculation. It does not say that confounding is exactly \(36.3\%\), nor does it identify the cause. Interpret it alongside the QQ plot, the study’s sample structure, and other available diagnostics. As discussed in the earlier material on polygenicity and inflation diagnostics, a single scalar can combine signals from different sources.
A common mistake is to treat any \(\lambda>1\) as proof of confounding and to “fix” it by scaling all test statistics without investigating the model. The reverse mistake is to declare an analysis safe because \(\lambda\) is near \(1\). Neither conclusion follows from the statistic alone. The QQ plot’s shape and departures across the distribution can add context, but they too are diagnostic evidence rather than a unique causal explanation of the pattern.
How to Diagnose and Compare Remedies
A disciplined check asks which assumption is at issue and what evidence would reveal a problem. The following sequence is useful before interpreting association results.
Confirm that individuals align across phenotype, genotype, covariate, and relationship-matrix inputs; check missingness, coding, and matrix validity.
Read convergence messages, inspect variance estimates and optimizer diagnostics, and assess whether the fit is stable under reasonable checks.
Examine the QQ plot and inflation summaries without assigning a unique cause from either one alone.
When appropriate, compare a relatedness-aware model with simpler adjustments or alternative justified relationship matrices, and investigate important differences rather than choosing the most favorable result.
An ordinary regression is simpler but assumes independent errors unless its uncertainty calculation accounts for dependence. Adjusting for a small set of covariates can address measured structure, but does not necessarily represent residual genetic covariance. An LMM models covariance through \(K\), but depends on the adequacy of that matrix and the numerical fit. Inflation rescaling is a summary adjustment, not a replacement for diagnosing why the statistics depart from their null reference. These approaches are not interchangeable: the appropriate comparison depends on the study design and the source of dependence or miscalibration.
Diagnostics also have limits. QQ plots and \(\lambda\) summarize patterns across tests; they cannot validate each variant’s result or prove that a selected model is correct. Sensitivity to model choices is informative but does not establish which choice is true. And even a well-calibrated association test estimates evidence for association conditional on the model; it does not by itself show that a variant causes a phenotype.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In the two-observation example, what variance and standard error result from the independence assumption, and why are they wrong under the stated correlation?
- A relationship matrix has the correct dimensions but is ordered differently from the phenotype vector. Why is the analysis still at risk?
- A model reports a non-convergence warning and a small SNP p-value. What can be concluded immediately, and what should be checked?
- A GWAS has \(\lambda=1.36\). Does this establish that confounding caused the inflation? Explain.
Solutions. (1) Independence gives variance \((1+1)/4=0.50\) and standard error \(\sqrt{0.50}=0.7071\). With correlation \(0.50\), the covariance is \(0.50\), so the correct variance for this average is \((1+1+2(0.50))/4=0.75\), with standard error \(0.8660\); omitting the covariance understates uncertainty in this example. (2) The matrix entries can be assigned to the wrong individuals, so the fitted covariance does not describe the intended sample even though its size is correct. Align all inputs by individual identifiers and verify the reordered matrix. (3) The printed p-value is not yet established as reliable. Inspect the warning, estimates, gradients or other optimizer diagnostics, and stability under reasonable refitting checks before drawing an association conclusion. (4) No. The value summarizes median test-statistic inflation relative to a null reference; it does not identify its cause. Consider the QQ-plot pattern and study context as well as other diagnostics.