Gaussian Mixed Models for Genotype-Based Covariance

Why Genotype-Based Covariance Matters

In a genetic study, individuals with similar genotypes may also have correlated phenotypes because their genetic backgrounds share effects across many markers. A model that treats every individual's phenotype as independent can overlook this covariance. The Gaussian mixed model represents it through a random genetic effect whose covariance is specified by a genomic relationship matrix.

The central idea is precise: the relationship matrix does not directly predict a phenotype or establish that one individual's genotype causes another's outcome. Instead, it describes a pattern of covariance among genetic effects. Combining that covariance with independent residual variation gives the marginal covariance of the observed phenotypes.

Model and Assumptions

Definition: Let $y$ be an $n$-vector of quantitative phenotypes, $X$ an $n \times p$ fixed-effect design matrix, and $\beta$ a $p$-vector of fixed-effect coefficients. The Gaussian genotype-based mixed model is $$ y = X\beta + g + e, $$ where $g$ is an $n$-vector of polygenic effects with $g \sim N(0,\sigma_g^2 K)$, and $e$ is an $n$-vector of residual effects with $e \sim N(0,\sigma_e^2 I)$. Here, $K$ is a specified $n \times n$ genomic relationship matrix, $I$ is the $n \times n$ identity matrix, and $\sigma_g^2,\sigma_e^2 \ge 0$ are variance components.

The model uses $K$ as an input matrix summarizing genotype-based relationships among the $n$ individuals. It must be symmetric and positive semidefinite so that $\sigma_g^2 K$ can be a covariance matrix. The particular convention used to construct $K$ determines the scale and interpretation of its entries; its construction and scaling are separate choices from the covariance model itself.

The standard Gaussian specification assumes that $X$ is treated as fixed, that $\beta$ is a fixed but unknown parameter, and that $g$ and $e$ are independent. The residual assumption $e \sim N(0,\sigma_e^2 I)$ means residuals have mean zero, common variance $\sigma_e^2$, and zero covariance between distinct individuals. The polygenic effects have mean zero and covariance $\sigma_g^2 K$. These assumptions describe the model; they are not guaranteed merely because a relationship matrix has been supplied.

It is useful first to view the model conditional on $g$. Given the genetic effects, the phenotype vector has mean $X\beta+g$ and residual covariance $\sigma_e^2 I$. Once $g$ is treated as random and averaged over, the phenotypes have a different covariance: residual variation remains, and covariance from the shared genetic structure is added.

Deriving the Marginal Covariance

The marginal distribution is the distribution of $y$ after accounting for the randomness in both $g$ and $e$, rather than conditioning on a particular value of $g$. Since the sum of independent Gaussian vectors is Gaussian, the marginal distribution is Gaussian as well.

\begin{aligned} E(y) &= E(X\beta+g+e) \\ &= X\beta+E(g)+E(e) \\ &= X\beta, \\ \operatorname{Cov}(y) &= \operatorname{Cov}(g+e) \\ &= \operatorname{Cov}(g)+\operatorname{Cov}(e) \\ &\quad +\operatorname{Cov}(g,e)+\operatorname{Cov}(e,g) \\ &= \sigma_g^2 K+\sigma_e^2 I. \end{aligned}

The cross-covariance terms are zero because $g$ and $e$ are independent. The fixed part $X\beta$ contributes to the mean, not to the covariance, because it is nonrandom under this model. Writing $V=\sigma_g^2 K+\sigma_e^2 I$, the full marginal statement is

y \sim N(X\beta,V).

Thus the $i$th diagonal entry of $V$ is the marginal variance of individual $i$'s phenotype, while the $(i,j)$th off-diagonal entry is the covariance between individuals $i$ and $j$. In particular, for $i \ne j$, the residual covariance is zero under the stated model, so that $V_{ij}=\sigma_g^2K_{ij}$. The diagonal is $V_{ii}=\sigma_g^2K_{ii}+\sigma_e^2$. This distinction is central: the relationship matrix supplies genetic covariance, while the residual term supplies individual-specific noise.

Worked Example: Two Individuals With a Shared Genetic Covariance

Suppose a hypothetical study has two individuals and uses $K=\begin{pmatrix}1&0.25\\0.25&1\end{pmatrix}$. Let $\sigma_g^2=0.4$ phenotype-units squared and $\sigma_e^2=0.6$ phenotype-units squared. Assume the model's Gaussian and independence conditions hold.

The marginal covariance is $V=0.4K+0.6I$. Substituting the matrices gives

V = 0.4\begin{pmatrix}1&0.25\\0.25&1\end{pmatrix} +0.6\begin{pmatrix}1&0\\0&1\end{pmatrix} = \begin{pmatrix}0.4&0.1\\0.1&0.4\end{pmatrix} + \begin{pmatrix}0.6&0\\0&0.6\end{pmatrix} = \begin{pmatrix}1&0.1\\0.1&1\end{pmatrix}.

Each phenotype has marginal variance $1$ phenotype-units squared. Their covariance is $0.1$ phenotype-units squared, coming from $0.4 \times 0.25$; the residual component adds nothing to their off-diagonal covariance. Because both variances are $1$, their marginal correlation is $0.1/\sqrt{1 \times 1}=0.1$. This is the correlation implied by this complete model, not a claim that the matrix entry alone is a phenotype correlation.

Reading the Entries of the Covariance Matrix

The formula for $V$ turns the structure of $K$ into a direct statement about phenotype covariance. When $\sigma_g^2>0$, a positive $K_{ij}$ implies positive genetic covariance between individuals $i$ and $j$; a negative $K_{ij}$ implies negative genetic covariance. If $K_{ij}=0$, the model assigns no genetic covariance to that pair. These statements concern covariance under the chosen model and matrix convention. They do not mean that every entry is a correlation, nor that a pairwise relationship causes a phenotype.

The diagonal entries matter as well. If $K_{ii}$ differs across individuals, then their modeled genetic variances differ: the genetic contribution for individual $i$ is $\sigma_g^2K_{ii}$. Since the residual variance is common in this model, the total marginal variance can also differ across individuals. A unit diagonal is not required by the covariance equation itself.

Worked Example: A Three-Individual Covariance Pattern

Consider a hypothetical phenotype measured on three individuals, with $K=\begin{pmatrix}1&0.5&0\\0.5&1&0.2\\0&0.2&1\end{pmatrix}$, $\sigma_g^2=0.3$ phenotype-units squared, and $\sigma_e^2=0.7$ phenotype-units squared. Assume $K$ is the specified symmetric positive-semidefinite matrix and the residuals are independent of the genetic effects.

The marginal covariance is

\begin{aligned} V &= 0.3K+0.7I \\ &= \begin{pmatrix} 0.3&0.15&0\\ 0.15&0.3&0.06\\ 0&0.06&0.3 \end{pmatrix} + \begin{pmatrix} 0.7&0&0\\ 0&0.7&0\\ 0&0&0.7 \end{pmatrix} \\ &= \begin{pmatrix} 1&0.15&0\\ 0.15&1&0.06\\ 0&0.06&1 \end{pmatrix}. \end{aligned}

The model gives each phenotype variance $1$ phenotype-units squared. Individuals 1 and 2 have covariance $0.3 \times 0.5=0.15$ phenotype-units squared, and individuals 2 and 3 have covariance $0.3 \times 0.2=0.06$ phenotype-units squared. Individuals 1 and 3 have zero covariance under this model because their $K$ entry is zero and residuals are independent. The zero is a model implication, not evidence that the two phenotypes could never be associated for other reasons omitted from this model.

Contrasts and the Role of the Diagonal

A contrast such as the difference between two phenotypes depends on both their marginal variances and their covariance. For any fixed vector $a$, the variance of the linear combination $a^\mathsf{T}y$ is $a^\mathsf{T}Va$. This follows from the covariance identity for a linear transformation: $\operatorname{Cov}(Ay)=A\operatorname{Cov}(y)A^\mathsf{T}$, here with $A=a^\mathsf{T}$. For the difference $y_i-y_j$, this becomes $V_{ii}+V_{jj}-2V_{ij}$.

Worked Example: Variance of a Pairwise Phenotype Difference

For a hypothetical pair of individuals, suppose $K=\begin{pmatrix}1.2&0.3\\0.3&0.8\end{pmatrix}$, $\sigma_g^2=0.5$ phenotype-units squared, and $\sigma_e^2=0.25$ phenotype-units squared. The matrix is symmetric and positive definite: its diagonal entries are positive and its determinant is $1.2 \times 0.8-0.3^2=0.87>0$. Assume the genetic effects and residuals are independent as specified by the model.

First calculate the marginal covariance:

V = 0.5\begin{pmatrix}1.2&0.3\\0.3&0.8\end{pmatrix} +0.25\begin{pmatrix}1&0\\0&1\end{pmatrix} = \begin{pmatrix}0.85&0.15\\0.15&0.65\end{pmatrix}.

Therefore, the variance of the difference is

\operatorname{Var}(y_1-y_2)=0.85+0.65-2(0.15)=1.20

phenotype-units squared. Equivalently, the genetic contribution to this variance is $0.5(1.2+0.8-2(0.3))=0.70$ phenotype-units squared, and the residual contribution is $0.25+0.25=0.50$ phenotype-units squared, for a total of $0.70+0.50=1.20$. This calculation shows why the diagonal entries cannot be discarded when interpreting genetic covariance.

Why This Model Matters—and a Common Mistake

Genotype-based covariance is useful because it represents polygenic similarity without assigning the entire genetic contribution to one marker at a time. It can account for dependence among phenotypes associated with related genetic backgrounds while retaining fixed effects such as an intercept or measured covariates in $X\beta$. The covariance matrix $V$ is the model's summary of how variable the phenotype vector is around its mean. This tutorial specifies that structure; estimating its variance components or using it for inference requires additional procedures.

Common mistake: A reader sees a positive off-diagonal entry in $K$ and concludes that the two individuals' phenotypes must have that same correlation. This is incorrect. The off-diagonal covariance is $\sigma_g^2K_{ij}$, not $K_{ij}$ by itself. Correlation also divides by the square root of the two marginal variances, which include both genetic and residual contributions. For example, in the two-individual example, the entry $K_{12}=0.25$ yielded covariance $0.1$ and correlation $0.1$, not $0.25$.

A second interpretive boundary is equally important: covariance is not causation. The model describes dependence under its assumptions; it does not show that one individual's genotype affects another individual's phenotype. If important sources of dependence or variation are omitted, the specified covariance may not adequately describe the data.

Key takeaway: Under $y=X\beta+g+e$, with $g \sim N(0,\sigma_g^2K)$, $e \sim N(0,\sigma_e^2I)$, and independent genetic and residual effects, the marginal phenotype vector is $y \sim N(X\beta,\sigma_g^2K+\sigma_e^2I)$. The relationship matrix shapes genetic covariance; the residual term adds individual-specific variance.

Check Your Understanding

Use the model assumptions stated above to answer the following questions.

Why does the fixed-effect term $X\beta$ contribute to the marginal mean but not to the covariance matrix $V$?
If $K_{ij}=0.4$ and $\sigma_g^2=0.2$ phenotype-units squared, what is the modeled covariance between distinct individuals $i$ and $j$, assuming residuals are independent?
What assumption makes the cross-covariance terms between $g$ and $e$ equal to zero?
Explain why an off-diagonal entry of $K$ is not automatically the correlation between two phenotypes.
For two individuals, write the formula for the variance of their phenotype difference using entries of $V$.