Worked Analysis Of Individual-Level Genetic Correlation

Why Analyze Two Traits Together?

When the same people have measurements on two quantitative traits, an individual-level analysis can use both those measurements and the genomic relationships among participants. The target is often the genetic covariance between the traits: whether people who are more genetically similar also tend to be similar across the two traits, after accounting for modeled covariates and residual variation. Dividing that covariance by the genetic standard deviations gives the genetic correlation, a scale-free measure of the modeled shared additive genetic signal.

This is a model-based analysis. Its interpretation depends on the relationship matrix, the covariance structure, the fixed effects, and assumptions about random effects and residuals. A fitted genetic correlation describes covariance under that model; it does not establish that changing one trait would cause a change in the other.

The Bivariate Genomic Mixed Model

Suppose $n$ participants have measurements on traits 1 and 2. Let $y_1$ and $y_2$ be their $n$-vectors of measurements, $X_1$ and $X_2$ the fixed-effect design matrices, and $\beta_1$ and $\beta_2$ the corresponding fixed-effect coefficients. These may include an intercept and measured covariates such as age, sex, or ancestry variables. Let $K$ be the known $n\times n$ genomic relationship matrix, whose entries summarize genetic similarity between pairs of participants.

Write $u_1$ and $u_2$ for the vectors of additive genetic effects and $e_1$ and $e_2$ for residual effects. The model is

y_1=X_1\beta_1+u_1+e_1,\qquad y_2=X_2\beta_2+u_2+e_2.

The bivariate genetic effects have covariance structure $G\otimes K$, while the residual effects have covariance structure $R\otimes I_n$, where $I_n$ is the $n\times n$ identity matrix. In trait-block form, these assumptions mean

\operatorname{Cov}(u_1,u_1)=\sigma_{g1}^2K,\quad \operatorname{Cov}(u_2,u_2)=\sigma_{g2}^2K,\quad \operatorname{Cov}(u_1,u_2)=\sigma_{g12}K,

\operatorname{Cov}(e_1,e_1)=\sigma_{e1}^2I_n,\quad \operatorname{Cov}(e_2,e_2)=\sigma_{e2}^2I_n,\quad \operatorname{Cov}(e_1,e_2)=\sigma_{e12}I_n.

Thus $\sigma_{g1}^2$ and $\sigma_{g2}^2$ are genetic variances, and $\sigma_{g12}$ is the genetic covariance. The corresponding residual components describe variation not captured by the fixed effects and modeled genetic effects. Genetic and residual covariance are distinct model components: a nonzero residual covariance does not by itself imply a nonzero genetic covariance.

Definition: Provided both genetic variances are positive, the genetic correlation is $$ r_g=\frac{\sigma_{g12}}{\sqrt{\sigma_{g1}^2\sigma_{g2}^2}}. $$ It is dimensionless and lies between $-1$ and $1$ when the genetic covariance matrix is valid. The covariance $\sigma_{g12}$, in contrast, has units equal to the product of the two traits' units.

The usual model assumes that the genetic effects are jointly mean-zero Gaussian with covariance $G\otimes K$, and the residual effects are jointly mean-zero Gaussian with covariance $R\otimes I_n$. It also assumes genetic effects and residuals are independent, conditional on the model, and that the fixed effects adequately represent modeled mean differences. The genetic covariance matrix $G$ and residual covariance matrix $R$ must each be positive semidefinite. If a model includes missing trait measurements, estimation must use the observed-data covariance structure; the examples below use participants with both traits measured.

Worked Analysis From Fitted Components

Worked Example: Interpreting a Bivariate Model Fit

Suppose a hypothetical study has 600 participants with both traits measured. The traits are measured in standardized units, and the genomic relationship matrix and fixed-effect covariates have been specified. A bivariate REML fit returns the following estimates:

\widehat G= \begin{pmatrix} 0.40&0.18\\ 0.18&0.25 \end{pmatrix}, \qquad \widehat R= \begin{pmatrix} 0.60&0.12\\ 0.12&0.75 \end{pmatrix}.

The fitted genetic variances are 0.40 and 0.25 squared standardized units, and the genetic covariance is 0.18 squared standardized units. The residual variances are 0.60 and 0.75, with residual covariance 0.12. Before interpreting these estimates, check that the covariance matrices are admissible. For a two-by-two symmetric matrix, positive diagonal entries and a nonnegative determinant ensure positive semidefiniteness. Here,

\det(\widehat G)=0.40(0.25)-(0.18)^2=0.1000-0.0324=0.0676>0,

\det(\widehat R)=0.60(0.75)-(0.12)^2=0.4500-0.0144=0.4356>0.

Both matrices have positive diagonal entries and determinants, so both are positive definite. The genetic correlation estimate is

\widehat r_g= \frac{0.18}{\sqrt{0.40(0.25)}}= \frac{0.18}{\sqrt{0.1000}}= \frac{0.18}{0.3162}=0.5692.

Under this fitted model, the traits have a moderately positive genetic correlation. The estimated genetic covariance is positive, and the correlation expresses its magnitude relative to the two genetic standard deviations. It is not the correlation between the observed trait values: the residual covariance also contributes to their observed covariation.

Uncertainty for the Genetic Correlation

A point estimate alone does not show how precisely the genetic correlation is estimated. A common large-sample approach is the delta method: approximate the variance of a smooth function of estimated covariance components by propagating their estimated sampling variance through the function's derivatives. Define $A=\sigma_{g1}^2$, $B=\sigma_{g2}^2$, and $C=\sigma_{g12}$, so that $r_g=C(AB)^{-1/2}$. Its derivatives are

\frac{\partial r_g}{\partial A}=-\frac{r_g}{2A},\qquad \frac{\partial r_g}{\partial B}=-\frac{r_g}{2B},\qquad \frac{\partial r_g}{\partial C}=\frac{1}{\sqrt{AB}}.

If $\widehat\theta=(\widehat A,\widehat B,\widehat C)^\mathsf{T}$ has estimated covariance matrix $\widehat V$, let $d$ be the vector of these derivatives evaluated at the fitted components. The delta-method variance estimate is $d^\mathsf{T}\widehat Vd$, and its square root is the approximate standard error of $\widehat r_g$. A 95% Wald interval is $\widehat r_g\pm1.96\,SE(\widehat r_g)$. This approximation is most useful when the component estimates are away from boundary values and the sampling distribution is reasonably close to normal.

Worked Example: Calculating a Delta-Method Interval

Continue with the hypothetical estimates $A=0.40$, $B=0.25$, and $C=0.18$. Suppose the fitted model reports the following illustrative estimated sampling covariance matrix for these three component estimates. The zero off-diagonal entries here are a simplifying hypothetical feature of this example, not a general property of REML estimates.

\widehat V= \begin{pmatrix} 0.001600&0&0\\ 0&0.001225&0\\ 0&0&0.001225 \end{pmatrix}.

At $\widehat r_g=0.5692$, the gradient components are

\frac{\partial r_g}{\partial A}=-\frac{0.5692}{2(0.40)}=-0.7115,\quad \frac{\partial r_g}{\partial B}=-\frac{0.5692}{2(0.25)}=-1.1384,\quad \frac{\partial r_g}{\partial C}=\frac{1}{\sqrt{0.40(0.25)}}=3.1623.

Because the illustrative component-estimate covariances are zero, the propagated variance is the sum of the three derivative-squared contributions:

\widehat{\operatorname{Var}}(\widehat r_g) =(-0.7115)^2(0.001600)+(-1.1384)^2(0.001225)+(3.1623)^2(0.001225) =0.01465.

Therefore $SE(\widehat r_g)=\sqrt{0.01465}=0.1210$, and the approximate 95% interval is

0.5692\pm1.96(0.1210) =0.5692\pm0.2372 =(0.3320,\ 0.8064).

The interval describes uncertainty in the model-based genetic correlation: under the usual repeated-sampling interpretation and the approximation used, the procedure would cover the underlying parameter in about 95% of comparable studies. It does not assign a 95% probability to the fixed parameter being inside this particular interval.

Testing Zero Genetic Covariance

To assess evidence that the modeled genetic covariance differs from zero, test $H_0:\sigma_{g12}=0$ against $H_a:\sigma_{g12}\ne0$. When the two genetic variances are positive, zero covariance is equivalent to zero genetic correlation. A Wald test divides the estimated covariance by its standard error and compares the result with a standard normal reference distribution. This is an approximate large-sample test and requires a valid standard error from the fitted model.

At the null value, the covariance parameter is locally unconstrained in either direction when the two variances are strictly positive: sufficiently small positive or negative covariances still give a positive definite $G$. Thus a two-sided one-degree-of-freedom Wald test is appropriate under the usual regularity conditions. If a variance is at a boundary or the fit is otherwise nonregular, ordinary Wald calibration may be unreliable.

Worked Example: Testing Whether Genetic Covariance Is Zero

For the hypothetical fitted model above, suppose the estimated standard error of $\widehat{\sigma}_{g12}=0.18$ is 0.035. At significance level $\alpha=0.05$, calculate the Wald statistic:

z=\frac{\widehat{\sigma}_{g12}-0}{SE(\widehat{\sigma}_{g12})} =\frac{0.18}{0.035}=5.143.

The two-sided normal-reference p-value is $2P(Z\geq5.143)$, approximately $2.7\times10^{-7}$, where $Z$ has a standard normal distribution. Since this is below 0.05, reject $H_0$. The data provide strong evidence against zero genetic covariance under the specified model. This test does not establish that the covariance is causal, nor does it assess whether the model's assumptions are correct.

Residual Diagnostics and Model Checks

Inference from a fitted covariance structure is only as credible as the model's adequacy. Diagnostics should examine whether the fitted mean and covariance structures leave systematic patterns. Useful checks include plotting residuals against fitted values and covariates, looking for outliers or changes in spread, and checking whether residual cross-trait association remains after accounting for the fitted genetic effects. Because the model includes related participants, residual patterns should also be considered in relation to $K$, not only as if all observations were independent.

Residuals used for screening are often conditional residuals, calculated after estimating fixed effects and predicting genetic effects. Those predicted effects are shrunk toward zero, so conditional residuals are not independent observations with exactly the raw residual covariance $R$. A residual plot or sample residual correlation is therefore a diagnostic signal, not a replacement estimator for $\sigma_{e12}$. Model-based or appropriately standardized residual checks should account for the fitted covariance structure.

Worked Example: Screening Paired Residuals

Suppose a hypothetical diagnostic extract gives six paired conditional residuals, in standardized trait units:

Participant	Trait 1 residual	Trait 2 residual
1	-1	-0.8
2	-1	0.8
3	0	-0.8
4	0	0.8
5	1	-0.8
6	1	0.8

Both residual means are zero. The sum of cross-products is

\sum_{i=1}^{6} e_{i1}e_{i2} =(-1)(-0.8)+(-1)(0.8)+(0)(-0.8)+(0)(0.8)+(1)(-0.8)+(1)(0.8)=0.

Consequently, the sample residual correlation is zero for these six pairs. This small calculation does not prove the residual model is correct: the sample is tiny, and the residuals are conditional estimates. A complete check would also inspect residuals against fitted values and covariates, investigate unusually large standardized residuals, and assess whether residual structure is related to genomic similarity or participant grouping.

A Practical Analysis Sequence

Specify the paired-trait model.
Define outcomes, covariates, the genomic relationship matrix, and the genetic and residual covariance structures. State the assumptions before interpreting fitted components.

Fit and validate covariance components.
Estimate the model, then check that the fitted genetic and residual covariance matrices are positive semidefinite and that the fit has not encountered a boundary or convergence problem.

Compute the genetic correlation and uncertainty.
Divide the estimated genetic covariance by the product of the genetic standard deviations, and use the model's estimated component covariance matrix to obtain an interval.

Test the covariance and inspect diagnostics.
Report the test, p-value, and contextual conclusion, then check whether residual patterns raise concerns about the fitted mean or covariance assumptions.

Common Pitfalls

Calling every cross-trait association genetic. The observed traits can covary through genetic effects, residual effects, or both. Interpret $\sigma_{g12}$ specifically as the model's genetic covariance component.
Reporting a correlation without its components. A value of $\widehat r_g$ is easier to interpret alongside the estimated genetic variances and covariance, with units stated for the covariance.
Ignoring covariance-matrix validity. A fitted covariance outside the positive-semidefinite region is not a valid covariance structure. Check the matrix, convergence, and any boundary estimates before drawing conclusions.
Reading a small p-value as causal evidence. Rejecting zero covariance is evidence against that null under the model and test assumptions; it does not show that one trait causes the other.
Treating a residual check as a formal proof. A clean plot cannot confirm every assumption, and a residual correlation from fitted residuals is not automatically an unbiased estimate of the residual covariance.

Key takeaway: Individual-level genetic correlation analysis jointly models both traits' genetic effects through the genomic relationship matrix. Interpret the fitted covariance components, correlation, uncertainty, and zero-covariance test together, and treat residual diagnostics as an essential check on the model rather than as a formality.

Check Your Understanding

Use the model and distinctions developed above to answer these questions.

Why does the genetic covariance $\sigma_{g12}$ have different units from the genetic correlation $r_g$?
What matrix checks would you make before interpreting fitted genetic and residual covariance estimates?
If the residual covariance is nonzero but the estimated genetic covariance is near zero, what does each component describe?
What does a two-sided Wald test of zero genetic covariance conclude when its p-value is below the chosen significance level?
Why is the sample correlation of conditional residuals not automatically an estimate of the model's residual correlation?