Mixed Models For Related Donors And Population Structure

Why Model Related Donors?

A linear expression QTL model treats the expression measurements from different donors as independent once the included covariates and genotype have been accounted for. That assumption can be implausible when donors are genetically related. Relatives, or donors with similar ancestry, may have more similar expression for many reasons beyond the particular variant being tested. If this residual similarity is ignored, the model can treat correlated observations as if they were independent and misstate the uncertainty in a QTL effect.

A linear mixed model represents this dependence through a covariance matrix derived from genetic relatedness. The genotype at the QTL remains a fixed effect of direct interest; a donor-level random effect represents residual similarity associated with genome-wide relatedness. The model therefore distinguishes the mean change in expression associated with the tested genotype from the covariance among donors’ expression measurements.

Model and Assumptions

Definition: For $n$ donors with one molecular phenotype measured per donor, let $y$ be the $n$-vector of phenotypes, $X$ the design matrix for an intercept and included covariates, and $g$ the $n$-vector of genotype dosages at the tested variant. A donor-level linear mixed model is $$ y = X\beta + g\alpha + u + \varepsilon, $$ where $\beta$ contains fixed covariate effects, $\alpha$ is the fixed QTL effect per dosage allele, $u\sim N(0,\tau K)$ is a donor-level random effect, and $\varepsilon\sim N(0,\sigma^2 I)$ is residual error. Assume $u$ and $\varepsilon$ are independent. Thus $y\mid X,g\sim N(X\beta+g\alpha,V)$, where $V=\tau K+\sigma^2 I$.

Here $K$ is a known, positive semidefinite genetic relatedness matrix, with one row and column per donor. Its entries encode similarity: larger off-diagonal values indicate greater modeled relatedness. A common construction is a genomic relationship matrix calculated from many genetic markers. Scaling conventions differ, so the numerical size of $\tau$ depends on the scale of $K$. The covariance $\tau K$ describes the part of phenotype covariance represented by genetic similarity; $\sigma^2 I$ describes independent residual variation. In particular, for two distinct donors $i$ and $j$, their modeled covariance is $\tau K_{ij}$.

The model assumes a suitable linear mean, a correctly specified covariance structure, and normally distributed random effects and errors for exact Gaussian likelihood-based inference. The residual errors have common variance $\sigma^2$ and are independent across donors. The genotype and covariates must be measured appropriately, and the fixed-effect design must contain enough independent information to estimate $\alpha$. In practice, variance components are estimated from the data, usually by maximum likelihood or restricted maximum likelihood (REML), rather than known in advance.

For a QTL analysis, $X$ commonly includes an intercept and prespecified technical or biological covariates. As in ordinary QTL regression, their selection and the interpretation of the genotype coefficient need scientific justification. A relatedness covariance is not a substitute for checking whether relevant fixed covariates, such as ancestry variables, are needed. The relatedness matrix may be built from genome-wide markers; analysts often exclude the tested region or nearby markers when constructing it, so that local signal at the tested locus does not unduly contribute to the background covariance.

Estimating the Fixed QTL Effect

Once the covariance parameters are specified, set $W=[X\ \ g]$ and $\theta=(\beta^\mathsf{T},\alpha)^\mathsf{T}$. If $V$ is positive definite and $W$ has full column rank, the generalized least-squares (GLS) estimator is

\hat{\theta}=(W^\mathsf{T}V^{-1}W)^{-1}W^\mathsf{T}V^{-1}y.

This estimator weights observations according to their covariance. It does not change the meaning of the fixed genotype coefficient: $\alpha$ is still the expected phenotype difference per additional dosage allele, conditional on the included fixed covariates. The model changes how that effect is estimated and how uncertain it is when donors are correlated.

Why GLS Has This Form

The estimator follows from the Gaussian model rather than from treating each donor as equally independent. For a fixed positive-definite $V$, the part of the negative log likelihood that depends on $\theta$ is one-half the weighted residual sum of squares, $(y-W\theta)^\mathsf{T}V^{-1}(y-W\theta)$. Differentiating it with respect to $\theta$ and setting the result to zero gives the GLS normal equations:

W^\mathsf{T}V^{-1}W\hat{\theta}=W^\mathsf{T}V^{-1}y.

Full column rank makes $W^\mathsf{T}V^{-1}W$ invertible, giving the displayed estimator. Under the stated Gaussian model with known $V$, its covariance is

\operatorname{Var}(\hat{\theta}\mid W)= (W^\mathsf{T}V^{-1}W)^{-1}.

In an actual analysis, $V$ depends on estimated variance components, so reported standard errors generally account for that estimation according to the fitting method. The known-$V$ expression is useful for understanding the role of covariance and for small calculations; it should not be mistaken for a claim that variance components are known in practice.

Worked Example: Constructing the Covariance Matrix

Suppose three hypothetical donors have a relatedness matrix

K= \begin{pmatrix} 1&0.4&0\\ 0.4&1&0\\ 0&0&1 \end{pmatrix}.

For illustration, take $\tau=2$ and $\sigma^2=1$. The diagonal variance for each donor is $\tau K_{ii}+\sigma^2=2(1)+1=3$. The covariance between donors 1 and 2 is $\tau K_{12}=2(0.4)=0.8$; the covariances involving donor 3 are zero. Therefore

V=2K+I= \begin{pmatrix} 3&0.8&0\\ 0.8&3&0\\ 0&0&3 \end{pmatrix}.

This matrix says that donors 1 and 2 have correlated phenotypes after accounting for the fixed effects, while donor 3 has no modeled residual covariance with them. It does not say that the tested QTL effect is larger for either related donor. The genotype effect is a separate term in the mean model.

Worked Example: Estimating a QTL Effect With Two Related Pairs

Consider four hypothetical donors, with phenotype values in normalized expression units $y=(4,6,8,10)^\mathsf{T}$ and genotype dosages $g=(0,0,2,2)^\mathsf{T}$. There is an intercept and no other covariate in this illustration. Suppose donors 1 and 2 form one related pair, donors 3 and 4 form another, and

K= \begin{pmatrix} 1&0.5&0&0\\ 0.5&1&0&0\\ 0&0&1&0.5\\ 0&0&0.5&1 \end{pmatrix}.

For a transparent calculation, take the variance components as known at $\tau=4$ and $\sigma^2=1$. Each pair’s covariance block is

4 \begin{pmatrix} 1&0.5\\ 0.5&1 \end{pmatrix} + \begin{pmatrix} 1&0\\ 0&1 \end{pmatrix} = \begin{pmatrix} 5&2\\ 2&5 \end{pmatrix}.

Within each pair the two donors have the same genotype. With equal diagonal variances and equal within-pair covariances, GLS estimates each genotype group’s mean by its pair average. The average phenotype at dosage 0 is $(4+6)/2=5$ expression units; at dosage 2 it is $(8+10)/2=9$ expression units. The fitted change across two dosage alleles is $9-5=4$ expression units, so the estimated effect per allele is

\hat{\alpha}=\frac{9-5}{2}=2 \quad\text{normalized expression units per dosage allele}.

The intercept estimate is 5 expression units, the estimated mean at dosage 0. For each pair, the variance of its average is $(5+5+2(2))/4=3.5$ expression units squared. The two pair averages are independent under this block-diagonal $V$; hence the variance of their difference is $3.5+3.5=7$. Dividing the difference by 2 to express it per allele gives

SE(\hat{\alpha})=\frac{\sqrt{7}}{2}=1.323 \quad\text{normalized expression units per dosage allele}.

This calculation treats the covariance parameters as known solely to make the GLS arithmetic explicit. In a fitted study model they would ordinarily be estimated. The result is an estimated conditional association, not evidence that changing the genotype would cause the expression difference.

Population Structure and the Meaning of Adjustment

A genetic relatedness covariance can account for residual correlation among donors, but covariance adjustment and confounding adjustment are not interchangeable. A fixed QTL effect compares expected phenotypes across genotype values under the specified mean model. If ancestry is associated both with genotype and with the phenotype mean, omitting an ancestry effect can leave the genotype coefficient confounded. Representing some ancestry-related similarity in $K$ does not, by itself, guarantee that the fixed-effect mean is correctly specified.

Worked Example: When Genotype and Ancestry Cannot Be Separated

Suppose four hypothetical donors have genotype dosages $g=(0,0,2,2)^\mathsf{T}$, ancestry scores $a=(0,0,1,1)^\mathsf{T}$, and expression values $y=(5,5,9,9)^\mathsf{T}$, in normalized expression units. Here $g=2a$ for every donor. A model with an intercept and genotype alone estimates a difference of $9-5=4$ expression units between dosage 0 and dosage 2, or $4/2=2$ units per allele.

But the same observations are exactly described by $y=5+4a$, with no separate genotype contribution. If both genotype and ancestry are included, their columns are linearly dependent: $g=2a$. The design matrix is not full rank, so the data cannot identify separate fixed effects for genotype and ancestry. Adding a relatedness covariance can model correlated residuals, but it cannot create information that distinguishes two perfectly collinear predictors. The reported genotype-only slope therefore cannot be interpreted as an ancestry-adjusted QTL effect in this example.

Inference and Practical Interpretation

A common test is $H_0:\alpha=0$ against an alternative in which $\alpha\ne 0$, using a Wald statistic or a likelihood-based comparison after fitting the covariance structure. A Wald statistic divides the estimated effect by its standard error; under an appropriate large-sample approximation it is compared with a standard normal distribution. The p-value describes how incompatible the observed statistic is with the null model under the specified assumptions. It is not the probability that the null hypothesis is true.

Report the effect in the phenotype’s units per dosage allele, alongside its standard error or confidence interval, and describe the fixed covariates and covariance model. A positive $\hat{\alpha}$ indicates higher expected phenotype values with each additional dosage allele, conditional on the modeled covariates. It does not establish a causal effect. The claim is also conditional on the chosen phenotype processing, sample, fixed-effect design, and assumptions about covariance.

The matrix $K$ is not a generic correction that guarantees valid inference. It must correspond to the donors in the same order as $y$, be positive semidefinite, and be constructed and scaled consistently. The variance component $\tau$ must be estimable from the available data, and the fixed design must allow the genotype effect to be distinguished from other fixed effects. If most genotype variation occurs between ancestry groups, the data may contain little information about a within-group genotype association even when the model can be fitted.

Common Mistakes

Wrong reasoning: “The random effect accounts for population structure, so no ancestry covariates are needed.” Correction: The random effect models covariance; it does not automatically remove a mean association between ancestry and phenotype. Consider the fixed-effect mean model and possible confounding separately.
Wrong reasoning: “A significant mixed-model coefficient is causal because relatedness was included.” Correction: The coefficient is a conditional association under the model. Covariance adjustment alone does not establish causation or eliminate every source of confounding.
Wrong reasoning: “The random effect is the QTL effect.” Correction: The tested genotype’s fixed coefficient $\alpha$ describes the conditional mean change per allele. The random effect $u$ represents residual donor-level variation with covariance $\tau K$.
Wrong reasoning: “Including genotype and ancestry always adjusts for both.” Correction: If their columns are collinear or nearly collinear, the effects cannot be separated precisely or at all. Check design rank and the amount of genotype variation within ancestry groups.

Key takeaway: A donor-level mixed model uses $V=\tau K+\sigma^2 I$ to represent correlated phenotypes while retaining the tested genotype as a fixed effect. Its QTL coefficient is a conditional per-allele association; valid interpretation still depends on an identifiable fixed-effect design and a suitable mean and covariance model.

Check Your Understanding

Try these questions before moving on.

In $V=\tau K+\sigma^2 I$, what determines the covariance between two distinct donors?
Why does a genetic relatedness random effect not automatically control confounding by ancestry?
For a genotype coded as 0, 1, or 2, what does a fixed coefficient of $-0.3$ mean if the phenotype is measured in normalized expression units?
What happens to separate genotype and ancestry estimates if the genotype column is exactly twice the ancestry column?
Why might an analysis exclude markers near the tested QTL when constructing the genetic relatedness matrix?