Worked LMM GWAS Analysis with Related Samples

Why Relatedness Changes the Association Analysis

In a GWAS of unrelated individuals, a regression model often treats the remaining trait variation as independent across people. That assumption can be untenable when the sample includes relatives: family members may resemble one another in the trait for genetic and shared-background reasons. If their dependence is ignored, standard errors and test statistics can be miscalibrated. A linear mixed model (LMM) represents this dependence in the covariance of the trait while testing the association of each variant.

The analysis below concerns a quantitative trait and an additive, single-variant test. It proceeds from a null model, which includes covariates and relatedness but not the tested variant, to a score test and then to genome-wide calibration checks. All numbers are hypothetical. A small calculation can illustrate the mechanics, but it cannot establish that a real GWAS has adequate power or valid calibration.

Definition: For $n$ individuals, let $y$ be the $n$-vector of quantitative trait values, $X$ the $n \times q$ matrix of fixed covariates, and $g$ the $n$-vector of dosages for one tested variant (usually coded 0, 1, or 2 copies of an allele). The additive LMM is $$ y = X\beta + g\beta_g + u + \varepsilon, $$ where $u \sim N(0,\sigma_g^2 K)$ is a random effect for genetic relatedness and $\varepsilon \sim N(0,\sigma_e^2 I)$ is residual error, independent of $u$. The known matrix $K$ describes pairwise relatedness, and $I$ is the identity matrix. Under the variant null hypothesis, $\beta_g=0$. The null-model covariance is $V=\sigma_g^2K+\sigma_e^2I$.

The fixed covariates commonly include an intercept and selected technical or population-structure covariates. The covariance matrix $K$ is usually constructed from genome-wide genotype data and is treated as fixed for the analysis. The model assumes that the conditional mean is linear in the included fixed effects, the random effects have the stated covariance, and residuals have the stated covariance. The normality assumption supports the usual likelihood-based fit and test calibration; the score test is generally interpreted using a large-sample reference distribution. A mixed model does not automatically correct every form of confounding or model misspecification.

Fit the Null Model Before Testing Variants

For each analysis, estimate the null-model variance components using the trait, covariates, and relatedness matrix, without adding the variant currently being tested. Restricted maximum likelihood (REML) is a common way to estimate the variance components while accounting for estimation of the fixed effects. With estimates $\widehat{\sigma}_g^2$ and $\widehat{\sigma}_e^2$, use $\widehat V=\widehat{\sigma}_g^2K+\widehat{\sigma}_e^2I$ for the subsequent test. In practice, the software performs the optimization; the analyst still needs to check that the fitted covariance and design are usable.

The algebra below requires a positive-definite covariance matrix $V$, so $V^{-1}$ exists, and a full-column-rank covariate matrix $X$, so $X^\mathsf{T}V^{-1}X$ is invertible. These are conditions, not automatic consequences of writing down an LMM. For example, a positive residual variance and positive-semidefinite $K$ imply $V$ is positive definite; if a fitted covariance is singular, the displayed inverse-based formulas do not apply without a specified alternative treatment.

Definition: When $V$ is positive definite and $X$ has full column rank, the generalized least-squares estimate under the null and the corresponding residual projection matrix are $$ \widehat{\beta}_0=(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}y, \qquad P_V=V^{-1}-V^{-1}X(X^\mathsf{T}V^{-1}X)^{-1}X^\mathsf{T}V^{-1}. $$ The matrix $P_V$ removes the fitted covariate contribution in the covariance-weighted geometry. In particular, $P_VX=0$. In applications using estimated variance components, these expressions are evaluated at the fitted null covariance.

A variance-component fit should be checked for convergence, plausible estimates, and a covariance matrix that is numerically usable. A boundary estimate is not automatically invalid, but may call for particular care with standard asymptotic approximations. Also verify that the trait and covariate rows align with the individuals and order used in $K$ and the genotype data. A mathematically valid fit to mismatched rows is not a valid analysis.

Worked Example: Building the Fitted Covariance

Suppose a hypothetical cohort contains two pairs of relatives. For a four-person illustration, let $K$ have a 1 for each person's diagonal and for the pairwise entries within each pair, and 0 for entries between pairs: $$ K= \begin{pmatrix} 1&1&0&0\\ 1&1&0&0\\ 0&0&1&1\\ 0&0&1&1 \end{pmatrix}. $$ Assume a null-model REML fit reports $\widehat{\sigma}_g^2=0.5$ trait-units squared and $\widehat{\sigma}_e^2=1.0$ trait-units squared. Take $X$ to contain only an intercept.

The fitted covariance is $\widehat V=0.5K+1.0I$, giving $$ \widehat V= \begin{pmatrix} 1.5&0.5&0&0\\ 0.5&1.5&0&0\\ 0&0&1.5&0.5\\ 0&0&0.5&1.5 \end{pmatrix}. $$ Each two-by-two block has eigenvalues 2.0 and 1.0, both positive, so this $\widehat V$ is positive definite. The intercept column has rank 1, so the inverse needed for the GLS intercept exists. The example checks the fitted covariance and design; the reported variance-component estimates are hypothetical fit output, not estimates derived from these four entries alone. Such a tiny sample would not support reliable GWAS inference.

Test One Variant With a Score Statistic

Under the null model, the score for the additive variant effect is $U=g^\mathsf{T}P_Vy$, and its information is $I_g=g^\mathsf{T}P_Vg$. These quantities use the null-model fit: there is no need to refit the variance components under the alternative for each variant. The test is available only if $I_g>0$. If $I_g=0$, the variant has no distinguishable dosage variation after the covariate adjustment in this model, and the formula cannot test its effect.

Theorem: If $V$ is positive definite, $X$ has full column rank, and $g^\mathsf{T}P_Vg>0$, the one-degree-of-freedom score statistic is $$ T=\frac{(g^\mathsf{T}P_Vy)^2}{g^\mathsf{T}P_Vg}. $$ Under the null and standard large-sample regularity conditions, $T$ is approximately distributed as a chi-square random variable with 1 degree of freedom. The score-based effect estimate and its standard error are $\widehat{\beta}_{g,\mathrm{score}}=U/I_g$ and $SE=1/\sqrt{I_g}$, respectively. With estimated variance components, this is a plug-in test; the chi-square calibration is approximate, not an exact finite-sample guarantee.

The justification follows from the Gaussian likelihood with covariance held fixed at its null value. After accounting for $X$, the derivative of the log likelihood with respect to $\beta_g$ at zero is $U=g^\mathsf{T}P_Vy$, and the corresponding information is $I_g=g^\mathsf{T}P_Vg$. For a one-parameter score test, the squared score divided by its information gives $T$. The inverse conditions above ensure these expressions are defined; the additional positive-information condition ensures the denominator is nonzero.

Worked Example: A Related-Sample Variant Test

Use the fitted covariance from the previous example and hypothetical trait values $y=(2,4,5,7)^\mathsf{T}$, in trait units. Test a hypothetical variant with dosages $g=(0,1,0,1)^\mathsf{T}$, using an intercept in $X$. The inverse of each covariance block is $$ \begin{pmatrix}1.5&0.5\\0.5&1.5\end{pmatrix}^{-1} = \begin{pmatrix}0.75&-0.25\\-0.25&0.75\end{pmatrix}. $$

The GLS intercept is $\widehat{\beta}_0=9/2=4.5$ trait units: $V^{-1}\mathbf{1}$ has four entries equal to 0.5, so $\mathbf{1}^\mathsf{T}V^{-1}y=9$ and $\mathbf{1}^\mathsf{T}V^{-1}\mathbf{1}=2$. Applying $P_V$ to $y$ gives $P_Vy=(-1.75,0.25,-0.25,1.75)^\mathsf{T}$. Thus $U=g^\mathsf{T}P_Vy=0.25+1.75=2.0$.

The GLS intercept for $g$ is 0.5, and $P_Vg=(-0.5,0.5,-0.5,0.5)^\mathsf{T}$. Therefore $I_g=g^\mathsf{T}P_Vg=1.0>0$. The score estimate is $2.0/1.0=2.0$ trait units per dosage unit, its standard error is $1/\sqrt{1.0}=1.0$ trait unit per dosage unit, and $T=2.0^2/1.0=4.0$. Comparing with the 1-degree-of-freedom chi-square distribution gives an approximate two-sided $p$-value of 0.0455. This very small illustration is not a reliable stand-alone GWAS; it demonstrates the calculation, not the adequacy of its asymptotic calibration.

In a larger study, the effect estimate describes the modeled difference in mean trait value per additional allele copy, conditional on the fixed covariates and under the specified LMM. The p-value is the probability, under the null model and its reference approximation, of a statistic at least as extreme as the observed one. It is not the probability that the null hypothesis is true.

Genome-Wide Calibration Checks

A valid-looking test at one variant does not establish that the genome-wide analysis is calibrated. Inspect the distribution of test statistics, including a quantile-quantile (QQ) plot, and examine the genomic inflation factor $\lambda$. For one-degree-of-freedom chi-square statistics, $\lambda$ is commonly calculated as the median observed statistic divided by the median of the null chi-square distribution, approximately 0.4549. The earlier material on QQ plots and genomic inflation explains how to read these summaries; here the key point is that they should be checked after applying the LMM.

Worked Example: Interpreting a Genome-Wide Inflation Summary

Suppose a hypothetical scan reports a median one-degree-of-freedom chi-square statistic of 0.500 across the tested variants. Then $$ \lambda=\frac{0.500}{0.4549}=1.099. $$ The median statistic is about 9.9% above its null median. This is a signal to inspect the QQ plot and analysis details, not proof that relatedness remains uncontrolled. A departure can reflect residual confounding or model problems, but widespread polygenic association can also elevate test statistics. Conversely, a near-one $\lambda$ does not prove that every variant's p-value is valid or that the model is correctly specified.

Check whether departures appear broadly or mainly in the most significant tail, and consider whether the pattern is consistent across relevant subsets and technical groupings. These are diagnostic clues, not mechanical rules for changing p-values. A mixed model can reduce dependence-induced miscalibration, but it does not repair sample mismatches, an inappropriate $K$, unmodeled structure, phenotype problems, or an unsuitable mean model.

When a Variant Cannot Be Tested

The condition $g^\mathsf{T}P_Vg>0$ has a direct interpretation. The covariate-adjusted genotype must contain information not already explained by $X$. Since $P_V$ is positive semidefinite, this quantity cannot be negative; it is zero when the tested dosage lies in the covariate span, in the covariance-weighted sense relevant to the model. A zero value is not evidence for a zero genetic effect: the effect is not separately identifiable from the included covariates in that test.

Worked Example: No Residual Genotype Information

Suppose $X$ contains only an intercept and all four hypothetical individuals have dosage 1, so $g=(1,1,1,1)^\mathsf{T}=X$. Because $P_VX=0$, it follows that $P_Vg=0$ and $g^\mathsf{T}P_Vg=0$. The score statistic would require division by zero and is undefined. The correct conclusion is that this sample cannot distinguish a dosage effect from the intercept; it is not appropriate to report a score-test p-value by substituting zero information into the formula.

Alternatives, Limitations, and a Common Pitfall

Ordinary linear regression is simpler, but its usual independent-error standard errors do not represent the relatedness covariance specified by $K$. A mixed model incorporates that covariance and is the relevant comparison when dependence among individuals matters. Within the mixed-model framework, a Wald test estimates the variant effect under the alternative and divides it by its estimated standard error; a likelihood-ratio test compares fitted likelihoods under the null and alternative. A score test evaluates the null fit and is convenient for scanning many variants. The methods can differ in finite samples, especially when variance components are estimated or the model is near a boundary, so their results should not be treated as interchangeable in every setting.

The LMM also has limits. It depends on an appropriate relatedness matrix and covariance assumptions, and the additive single-variant model does not test every possible genetic effect. Imperfect genotype data, unmodeled technical effects, or a mismatched trait distribution may undermine inference. Calibration plots and $\lambda$ summarize patterns across tests; they cannot certify the validity of any particular association.

Pitfall: assuming the written model guarantees invertibility. A researcher might use the GLS and score formulas whenever a matrix is called a covariance matrix. But a covariance matrix can be singular, and covariate columns can be linearly dependent. Verify $V$ is positive definite and $X$ has full column rank before using the inverse formulas.
Pitfall: treating zero information as a null result. If $g^\mathsf{T}P_Vg=0$, the test is undefined, not nonsignificant. Check dosage variation after covariate adjustment and report that the effect is not estimable in this design.
Pitfall: reading $\lambda>1$ as proof of confounding. Polygenic signal can contribute to inflation. Interpret the QQ plot and inflation summary together with the model and study design rather than assigning a single cause from $\lambda$ alone.

Key takeaway: A related-sample GWAS uses a null LMM to model trait covariance, then tests each variant only when the covariance and covariate design support the required calculations and the variant has positive residual information. Genome-wide diagnostics help assess calibration, but they do not replace model checks.

Check Your Understanding

Answer each question, then compare with the solutions.

What role does $K$ play in the null-model covariance $V$?
State the matrix conditions required for the GLS estimate and the additional condition required for a variant score test.
If $U=3$ and $I_g=4$, calculate the score statistic, score estimate, and standard error.
A scan has $\lambda=1.00$. Does that establish that every variant test is calibrated? Explain.

Solutions. (1) $K$ encodes pairwise relatedness, and $\sigma_g^2K$ contributes the covariance attributed to the random genetic effect. (2) $V$ must be positive definite and $X$ must have full column rank; additionally, $g^\mathsf{T}P_Vg>0$. (3) $T=U^2/I_g=9/4=2.25$; the score estimate is $U/I_g=3/4=0.75$ trait units per dosage unit, and $SE=1/\sqrt{I_g}=0.5$ trait units per dosage unit. (4) No. A genome-wide median summary cannot certify individual tests or rule out other model problems; inspect the distribution and the analysis assumptions as well.