Worked LMM Analysis of a Related Cohort

From Related Samples to an Analysis Plan

When a quantitative phenotype is measured in a related cohort, participants’ outcomes are not generally independent: relatives may resemble one another because they share genetic background. A linear mixed model represents this dependence through a covariance matrix built from genomic relationships. The same model can include a candidate variant as a fixed effect, so its estimated association is adjusted for the modeled relatedness and other covariates.

A complete analysis requires more than a coefficient and a p-value. The analyst must state the mean model, explain the covariance model and its assumptions, report the variance-component estimates, and check whether the fitted model is an adequate description of the data. The examples below use explicitly hypothetical data and model output.

What the Covariance Model Says

The mean model \(X\boldsymbol{\beta}\) describes systematic differences in expected phenotype associated with the included fixed effects. It commonly includes an intercept, relevant measured covariates, and—in a candidate-variant analysis—the candidate genotype dosage. The random effect represents a shared polygenic contribution whose covariance across participants is specified by \(K\). The residual term represents remaining variation under the model.

For participants \(i\) and \(j\), the model implies

Thus, an off-diagonal entry \(K_{ij}\) contributes modeled covariance between two different participants. The diagonal of \(K\) also matters: when \(K_{ii}=1\), the modeled variance of participant \(i\)'s phenotype around its fixed-effect mean is \(\sigma_g^2+\sigma_e^2\). The interpretation of the two variance parameters depends on the scale and construction of \(K\); their values should therefore be reported together with how that matrix was formed.

A standard analysis assumes that the fixed-effect mean is adequately specified, \(K\) is appropriate for the analyzed participants and is positive semidefinite, and the residuals have the stated common variance and independence structure. Normality supports the usual likelihood-based estimation and approximate Wald inference. These are modeling assumptions, not facts guaranteed by fitting software. In a candidate-variant analysis, the variant’s dosage is included as a fixed effect; a leave-one-chromosome-out relationship matrix is one common way to reduce overlap between the candidate’s chromosome and the polygenic covariance model.

Estimating the Variance Components

In practice, the variance parameters are estimated from the phenotype data and the specified \(K\), typically by restricted maximum likelihood (REML). REML accounts for the estimation of fixed effects when estimating covariance parameters. Estimates are constrained to be nonnegative because they represent variances. An estimate at zero is a boundary estimate, not evidence that the corresponding source of variation is known to be absent.

The fitted covariance matrix is obtained by substituting the estimates into \(V=\sigma_g^2K+\sigma_e^2I\). Inference about fixed effects then uses this estimated covariance structure. Variance estimates are in squared phenotype units; their standard errors, when reported, are also in squared phenotype units. The estimates describe this model fit and should not be interpreted beyond the assumptions and scale used to construct \(K\).

Testing a Candidate Genetic Effect

Let \(G_i\) denote participant \(i\)'s dosage of the candidate allele, coded 0, 1, or 2. The candidate-variant model includes \(G\) as a column of \(X\), alongside the intercept and other prespecified covariates. Its coefficient \(\beta_G\) is the expected phenotype difference associated with a one-allele increase in dosage, conditional on the included covariates and the mixed-model covariance structure. It is an association estimate, not by itself a causal effect.

A Wald test evaluates \(H_0:\beta_G=0\) against \(H_a:\beta_G\ne0\). If the fitted coefficient is \(\hat{\beta}_G\) and its estimated standard error is \(SE(\hat{\beta}_G)\), the Wald statistic is their ratio. With a large-sample normal reference distribution, the two-sided p-value is the probability under a standard normal distribution of a statistic at least as extreme in absolute value. The normal approximation, correct mean and covariance specifications, and sufficiently reliable variance estimates matter to the calibration of this test.

A confidence interval using the same normal approximation is \(\hat{\beta}_G\pm1.96SE(\hat{\beta}_G)\) for a 95% interval. It describes values compatible with the data and model under the procedure’s repeated-sampling coverage properties; it is not a probability statement about a fixed parameter being inside this particular interval.

Checks, Reporting, and Limitations

A fitted model should be examined before its association result is treated as dependable. The checks should be chosen to probe the assumptions and the data features relevant to the analysis, rather than to produce a single pass-or-fail label. For example, residual plots can reveal trends in residual spread or mean, and a normal quantile plot can reveal marked departures from the assumed residual distribution. Because participants are related, it is also useful to examine whether residual similarity remains patterned by genomic relatedness after fitting the model. Calibration checks across null or control variants can help reveal systematic test-statistic inflation, though they do not identify every possible source of misspecification.

Model checks do not prove that the covariance structure is correct. Potential limitations include a poorly specified fixed-effect mean, population structure not adequately represented by covariates or \(K\), an unsuitable relationship matrix, heterogeneity in residual variance, influential observations, and genotype or phenotype measurement error. If the candidate variant contributes to the relationship matrix, its signal may partly enter the random effect; a leave-one-chromosome-out matrix can reduce this concern. The choice of matrix and covariates should be stated transparently.

A useful justification for the covariance formula follows directly from the model. Since the fixed effects are treated as nonrandom, they contribute no covariance. Independence of \(\boldsymbol{u}\) and \(\boldsymbol{\varepsilon}\) makes their cross-covariance zero. Therefore,

This identity explains why the same genomic matrix enters the covariance between participants and the marginal variance for each participant: the diagonal and off-diagonal entries encode different parts of the same covariance structure.

Common Mistake

• Wrong reasoning: “The candidate-variant p-value is small, so the variant causes the phenotype, and the model has accounted for all relatedness.”
• Correction: The test evaluates association conditional on the specified fixed effects and covariance model. A small p-value does not prove causality or guarantee that \(K\), the covariates, or residual assumptions are adequate. Report the estimate and uncertainty, state the model, and describe diagnostics and limitations.