Interpreting Relatedness Patterns in the Covariance Matrix

From Genomic Relatedness to Phenotype Covariance

A genomic relationship matrix summarizes how genetic similarity is distributed across the people in a study. In a mixed-model GWAS, that structure matters because people who are more related under the chosen matrix are modeled as having more similar random genetic effects, and therefore potentially more similar phenotypes. The relationship value itself is not a phenotype covariance: the variance component scales it, while residual variation also contributes to phenotype variance.

This tutorial interprets the patterns that arise in the covariance matrix. The focus is on what its diagonal and off-diagonal entries say about individual variances and pairwise dependence, and how groups of related individuals produce broader patterns. The interpretation depends on the matrix scaling convention, the variance components, and the assumptions of the model.

Model and assumptions: Let \(y\) be the vector of \(n\) phenotypes and \(X\) the design matrix of fixed covariates, including an intercept when appropriate. Write \(y=X\beta+u+\varepsilon\), where \(u\sim N(0,\sigma_g^2K)\) is the random genetic effect and \(\varepsilon\sim N(0,\sigma_e^2I)\) is the residual. Assume \(K\) is a known, symmetric, positive semidefinite genomic relationship matrix; \(\sigma_g^2\geq 0\) and \(\sigma_e^2>0\); and \(u\) and \(\varepsilon\) are independent. Conditional on \(X\), the mean is \(X\beta\) and the covariance is \(V=\sigma_g^2K+\sigma_e^2I\). Normality supports the usual Gaussian mixed-model likelihood and inference; the covariance identity itself follows from the stated covariance and independence assumptions.

The entries \(K_{ij}\) depend on the matrix's scaling convention, established in the earlier material on scaling a genomic relationship matrix. Under this model, \(K_{ij}\) contributes to the covariance between the genetic effects of individuals \(i\) and \(j\). It is not, by itself, a phenotype correlation. A larger positive \(K_{ij}\) means a larger modeled genetic covariance when \(\sigma_g^2\) is fixed, but the resulting phenotype correlation also depends on the diagonal entries and on residual variance.

Reading Diagonal and Off-Diagonal Entries

For an individual \(i\), the phenotype variance is \(V_{ii}=\sigma_g^2K_{ii}+\sigma_e^2\). For two distinct individuals \(i\) and \(j\), the phenotype covariance is \(V_{ij}=\sigma_g^2K_{ij}\), because the independent residuals contribute no off-diagonal covariance. Their phenotype correlation is the covariance divided by the product of their standard deviations.

\operatorname{Var}(y_i\mid X)=\sigma_g^2K_{ii}+\sigma_e^2,\qquad \operatorname{Cov}(y_i,y_j\mid X)=\sigma_g^2K_{ij}\quad(i\ne j),

\operatorname{Corr}(y_i,y_j\mid X) = \frac{\sigma_g^2K_{ij}} {\sqrt{(\sigma_g^2K_{ii}+\sigma_e^2)(\sigma_g^2K_{jj}+\sigma_e^2)}}.

If the diagonal entries are all 1, each phenotype has variance \(\sigma_g^2+\sigma_e^2\), and the correlation simplifies to \(\sigma_g^2K_{ij}/(\sigma_g^2+\sigma_e^2)\). Thus, even in this simple case, the phenotype correlation is generally smaller in magnitude than \(K_{ij}\) when residual variance is positive. If the diagonal entries differ, the same off-diagonal relationship value can correspond to different phenotype correlations for different pairs.

Worked Example: Comparing Two Pairs

Suppose a hypothetical study uses a relationship matrix with diagonal entries 1. Let \(\sigma_g^2=0.4\) phenotype-units squared and \(\sigma_e^2=0.6\) phenotype-units squared. One pair has \(K_{ij}=0.5\); another pair has \(K_{kl}=0.05\). Calculate each pair's phenotype covariance and correlation.

For either individual, the phenotype variance is \(0.4(1)+0.6=1.0\) phenotype-units squared. For the first pair, \(V_{ij}=0.4(0.5)=0.20\) phenotype-units squared, so its correlation is \(0.20/\sqrt{1.0(1.0)}=0.20\).

For the second pair, \(V_{kl}=0.4(0.05)=0.020\) phenotype-units squared, and its correlation is \(0.020/\sqrt{1.0(1.0)}=0.020\). The model therefore assigns the first pair ten times the phenotype covariance of the second pair. It does not assign a phenotype correlation of 0.5 to the first pair: 0.5 is the relationship value, whereas the modeled phenotype correlation is 0.20 after accounting for the variance components.

These are marginal covariances conditional on the fixed design \(X\). Including covariates in the mean model does not change the formula for \(V\). It does, however, change the mean being modeled, so covariance should not be confused with similarity in observed phenotype values: two people can have similar observations partly because their fixed covariates predict similar means.

How Group Structure Appears in the Matrix

A group of related individuals produces a recognizable block pattern: off-diagonal entries within a group may be larger than entries between groups. A covariance matrix translates that pattern into phenotype units through \(\sigma_g^2\), while the residual variance adds to the diagonal only. The matrix as a whole must also be a valid covariance matrix. In particular, \(K\) must be positive semidefinite, which means \(a^\mathsf{T}Ka\geq 0\) for every vector \(a\). With nonnegative variance components, this implies \(V\) is positive semidefinite; the condition \(\sigma_e^2>0\) makes \(V\) positive definite.

Worked Example: Two Related Pairs in One Matrix

Consider four hypothetical individuals arranged as two pairs. Let their relationship matrix have diagonal entries 1, relationship 0.6 within each pair, and relationship 0.1 between the pairs:

K= \begin{pmatrix} 1&0.6&0.1&0.1\\ 0.6&1&0.1&0.1\\ 0.1&0.1&1&0.6\\ 0.1&0.1&0.6&1 \end{pmatrix}.

To check positive semidefiniteness, first consider within-pair contrast vectors \((1,-1,0,0)^\mathsf{T}\) and \((0,0,1,-1)^\mathsf{T}\). Multiplication by \(K\) gives 0.4 times each vector, so both have eigenvalue 0.4. On vectors constant within each pair, of the form \((a,a,b,b)^\mathsf{T}\), the action reduces to the two-dimensional matrix \(\begin{pmatrix}1.6&0.2\\0.2&1.6\end{pmatrix}\). Its eigenvalues are \(1.6+0.2=1.8\) and \(1.6-0.2=1.4\). Thus the eigenvalues of \(K\) are 1.8, 1.4, 0.4, and 0.4, all nonnegative.

Now let \(\sigma_g^2=0.5\) and \(\sigma_e^2=0.5\), in phenotype-units squared. Each diagonal entry of \(V\) is \(0.5(1)+0.5=1.0\). Within-pair phenotype covariance is \(0.5(0.6)=0.30\); between-pair covariance is \(0.5(0.1)=0.050\). Since each phenotype variance is 1.0, those covariances are also the corresponding correlations, 0.30 and 0.050. The eigenvalues of \(V=0.5K+0.5I\) are obtained by transforming each eigenvalue \(\lambda\) of \(K\) to \(0.5\lambda+0.5\): they are 1.4, 1.2, 0.7, and 0.7. Every eigenvalue is positive, as required for a nonsingular covariance matrix.

The two larger eigenvalues correspond to patterns that vary at the pair level: a shared shift across all four individuals, or a contrast between the two pairs. The smaller eigenvalues correspond to contrasts within pairs. This is one way to see that covariance structure applies to combinations of observations as well as to individual pairs. It does not mean that a particular realized group must have similar observed phenotypes; it describes variation across hypothetical repetitions under the model.

Negative Entries and Valid Covariance

A negative off-diagonal value is not automatically an error. Under some centering and scaling conventions, a relationship value can be negative relative to the reference population. It then implies negative modeled genetic covariance for that pair, not a negative variance. The full matrix still must be positive semidefinite, and the resulting phenotype covariance must be compatible with valid variances. An isolated negative entry cannot be judged without considering the matrix as a whole.

Worked Example: A Negative Relationship Value

For two hypothetical individuals, suppose \(K=\begin{pmatrix}1&-0.2\\-0.2&1\end{pmatrix}\), \(\sigma_g^2=0.3\), and \(\sigma_e^2=0.7\), with variance components measured in phenotype-units squared. The eigenvalues of \(K\) are \(1-0.2=0.8\) and \(1+0.2=1.2\), so \(K\) is positive definite. The phenotype covariance matrix is

V=0.3K+0.7I = \begin{pmatrix} 1.00&-0.060\\ -0.060&1.00 \end{pmatrix}.

Each phenotype variance is \(0.3+0.7=1.00\) phenotype-units squared. The phenotype covariance is \(-0.060\) phenotype-units squared, and the correlation is \(-0.060/\sqrt{1.00(1.00)}=-0.060\). This valid model assigns a small negative covariance to the pair. The negative value does not mean either phenotype has a negative variance, and it does not establish a causal antagonism between the individuals.

Diagnostics, Alternatives, and Limits

Interpretation should be checked against both the matrix and the model. Useful checks include the following:

Confirm that \(K\) is symmetric and that its diagonal and off-diagonal values follow the intended scaling convention.
Check that \(K\) is positive semidefinite up to numerical tolerance. A materially negative eigenvalue signals a problem with the matrix or its construction, not a meaningful negative variance.
Compare the observed pattern of \(K\) with the implied pattern of \(V\). Under this model, off-diagonal covariances must be proportional to the corresponding entries of \(K\), with proportionality factor \(\sigma_g^2\).
Check whether the assumed residual structure is plausible. The model uses independent residuals with common variance; residual correlation or unequal residual variances require a different covariance specification.
Interpret variance components in the context of the phenotype scale. A covariance has squared phenotype units, whereas a correlation and a relationship value are dimensionless.

A frequent mistake is to call \(K_{ij}\) the correlation between phenotypes. The correction is to calculate \(V_{ij}=\sigma_g^2K_{ij}\), then divide by the square root of the two phenotype variances to obtain the correlation. Another mistake is to read a large positive covariance as evidence that relatedness caused the observed phenotype values. The covariance is a model-based description of joint variation, not a causal conclusion about a particular pair.

The full genomic relationship matrix is more flexible than a simple family-block model, because it can represent graded relationships across all pairs rather than assigning only “same family” or “different family” values. An independent-residual model is simpler still, but has zero off-diagonal phenotype covariance and cannot represent the dependence encoded by \(K\). Each alternative makes different assumptions: a block model imposes a chosen group structure, while the genomic model uses the specified matrix. No covariance pattern alone proves that the matrix captures every relevant source of dependence. Misspecified relatedness, omitted covariates, or residual dependence can make the fitted model inadequate.

Key takeaway: In the genomic mixed model, the relationship matrix supplies the pattern of genetic covariance, and the genetic variance component sets its scale. Residual variance adds to phenotype variances but, under the independent-residual assumption, not to covariances between distinct individuals.

Check Your Understanding

Work through each question, then compare your reasoning with the solutions.

With unit diagonal entries, \(\sigma_g^2=0.25\), and \(\sigma_e^2=0.75\), what are the phenotype covariance and correlation for a pair with \(K_{ij}=0.4\)?
For a pair with \(K_{ii}=1.2\), \(K_{jj}=0.8\), \(K_{ij}=0.3\), \(\sigma_g^2=0.5\), and \(\sigma_e^2=0.5\), calculate both phenotype variances, the covariance, and the correlation.
Why does a negative off-diagonal entry in \(K\) not, by itself, invalidate the covariance model?
In the two-pair matrix from the worked example, what are the within-pair and between-pair phenotype covariances when \(\sigma_g^2=0.2\), regardless of the residual variance?

Solutions. For question 1, the covariance is \(0.25(0.4)=0.10\) phenotype-units squared. Each phenotype variance is \(0.25+0.75=1.00\), so the correlation is \(0.10/\sqrt{1.00(1.00)}=0.10\). For question 2, the variances are \(0.5(1.2)+0.5=1.10\) and \(0.5(0.8)+0.5=0.90\) phenotype-units squared. The covariance is \(0.5(0.3)=0.15\) phenotype-units squared, and the correlation is \(0.15/\sqrt{1.10(0.90)}\), approximately 0.151. For question 3, validity depends on the positive-semidefiniteness of the full \(K\), not on every entry being nonnegative; with nonnegative variance components and the stated residual model, the resulting \(V\) is a valid covariance matrix. For question 4, within-pair covariance is \(0.2(0.6)=0.12\) phenotype-units squared and between-pair covariance is \(0.2(0.1)=0.020\) phenotype-units squared. Independent residual variance affects the diagonal but not either off-diagonal covariance.