Estimating Genetic Correlation And Its Uncertainty

From Genetic Covariance to Genetic Correlation

A genetic covariance describes how additive genetic effects for two traits vary together, but its numerical size depends on the measurement scales of both traits. Genetic correlation standardizes that covariance, making it easier to describe the direction and strength of shared additive genetic variation. Estimating it also requires uncertainty calculations: the estimate is a function of several estimated variance components, not a directly observed quantity.

Suppose a bivariate genetic model has additive genetic variances \(\sigma^2_{g1}\) and \(\sigma^2_{g2}\), and additive genetic covariance \(\sigma_{g12}\). The variances are nonnegative, and the covariance matrix must be positive semidefinite. The corresponding population genetic correlation is defined as follows.

Definition: The genetic correlation between traits 1 and 2 is \(r_g=\sigma_{g12}/\sqrt{\sigma^2_{g1}\sigma^2_{g2}}\), provided both genetic variances are positive. It is the genetic covariance divided by the product of the genetic standard deviations. It is dimensionless and lies between \(-1\) and \(1\).

The estimate is usually the plug-in estimate: substitute estimates of the three genetic components into the definition. Write the estimated variances as \(\widehat{\sigma}^2_{g1}\) and \(\widehat{\sigma}^2_{g2}\), and the estimated covariance as \(\widehat{\sigma}_{g12}\). Then

\widehat{r}_g = \frac{\widehat{\sigma}_{g12}} {\sqrt{\widehat{\sigma}^2_{g1}\widehat{\sigma}^2_{g2}}}.

This calculation presumes that the component estimates came from an appropriate fitted model and that both estimated genetic variances are positive. If the model estimates a valid positive-semidefinite genetic covariance matrix, its components satisfy \(\widehat{\sigma}_{g12}^2\leq \widehat{\sigma}^2_{g1}\widehat{\sigma}^2_{g2}\), so the plug-in correlation is in the allowed range. If a numerical procedure yields a value outside that range, investigate the model, constraints, or computation rather than reporting it as a valid correlation.

Worked Example: Calculating a Genetic Correlation

In hypothetical units, a fitted bivariate model estimates genetic variance \(0.40\) for trait 1, genetic variance \(0.25\) for trait 2, and genetic covariance \(0.10\). Calculate the estimated genetic correlation.

Substitute the estimates. \[ \widehat{r}_g=\frac{0.10}{\sqrt{(0.40)(0.25)}}. \] The variance product is \(0.40\times0.25=0.10\), so its square root is \(\sqrt{0.10}=0.3162\). Therefore \[ \widehat{r}_g=\frac{0.10}{0.3162}=0.3162. \]

The estimated genetic correlation is about \(0.316\), a positive, dimensionless estimate. The covariance-matrix constraint is satisfied: \(0.10^2=0.0100\), while \(0.40\times0.25=0.100\). This calculation describes estimated shared genetic variation; by itself it does not establish a causal relationship between the traits.

Delta-Method Uncertainty

The delta method approximates the sampling variance of a smooth function of estimated parameters using a first-order Taylor expansion. For genetic correlation, let \(\boldsymbol{\theta}=(v_1,v_2,c)^\mathsf{T}\), where \(v_1=\sigma^2_{g1}\), \(v_2=\sigma^2_{g2}\), and \(c=\sigma_{g12}\). Let \(\widehat{\boldsymbol{\theta}}\) be their estimated values, and let \(\widehat{\mathbf V}\) estimate the covariance matrix of those component estimates. This matrix must include their covariances as well as their variances when those covariances are nonzero.

For \(g(v_1,v_2,c)=c(v_1v_2)^{-1/2}\), the gradient is

\nabla g(v_1,v_2,c)= \begin{pmatrix} -\dfrac{g(v_1,v_2,c)}{2v_1}\\[4pt] -\dfrac{g(v_1,v_2,c)}{2v_2}\\[4pt] \dfrac{1}{\sqrt{v_1v_2}} \end{pmatrix}.

The estimated delta-method variance is the gradient evaluated at the component estimates, multiplied on both sides of the estimated component covariance matrix:

\widehat{\operatorname{Var}}(\widehat r_g) \approx \nabla g(\widehat{\boldsymbol{\theta}})^\mathsf{T} \widehat{\mathbf V} \nabla g(\widehat{\boldsymbol{\theta}}), \qquad SE(\widehat r_g)= \sqrt{\widehat{\operatorname{Var}}(\widehat r_g)}.

The approximation assumes the component estimates have a sufficiently regular, approximately joint sampling distribution and that the function is locally well approximated by its linear expansion. In particular, both genetic variances should be positive and not so close to zero that the denominator makes the function unstable. A commonly used approximate 95% Wald interval is \(\widehat r_g\pm1.96\,SE(\widehat r_g)\). Its coverage is approximate, not guaranteed; it can extend beyond the possible correlation range.

Why the Delta Formula Works

A first-order expansion around the true component vector \(\boldsymbol{\theta}\) gives \(g(\widehat{\boldsymbol{\theta}})\approx g(\boldsymbol{\theta})+ \nabla g(\boldsymbol{\theta})^\mathsf{T} (\widehat{\boldsymbol{\theta}}-\boldsymbol{\theta})\). If the component estimator has covariance matrix \(\mathbf V\), the variance of this linear term is \(\nabla g^\mathsf{T}\mathbf V\nabla g\). Replacing the unknown parameters and covariance matrix with their estimates yields the displayed standard-error formula. The approximation omits higher-order terms, which explains why it can be unreliable near boundaries or when component estimates are highly variable.

The derivatives follow directly by differentiating \(g=c(v_1v_2)^{-1/2}\). With respect to \(c\), the derivative is \((v_1v_2)^{-1/2}\). With respect to \(v_1\), it is \(-\tfrac12c v_1^{-3/2}v_2^{-1/2}=-g/(2v_1)\); with respect to \(v_2\), it is \(-\tfrac12c v_1^{-1/2}v_2^{-3/2}=-g/(2v_2)\). These are the three entries used in the variance calculation.

Worked Example: A Delta-Method Standard Error and Interval

Continue with hypothetical component estimates \(v_1=0.40\), \(v_2=0.25\), and \(c=0.10\), giving \(\widehat r_g=0.3162\). Suppose the estimated covariance matrix of these three component estimates, in the order \((v_1,v_2,c)\), is diagonal with entries \(0.0025\), \(0.0016\), and \(0.0009\). This assumption of zero estimated covariances is part of this example, not a general rule.

The gradient at the estimates is \[ \nabla g= \begin{pmatrix} -0.3162/(2\times0.40)\\ -0.3162/(2\times0.25)\\ 1/\sqrt{(0.40)(0.25)} \end{pmatrix} = \begin{pmatrix} -0.3953\\ -0.6325\\ 3.1623 \end{pmatrix}. \] Because the estimated covariance matrix is diagonal, the variance estimate is the sum of each diagonal component variance multiplied by the squared corresponding gradient entry: \[ (-0.3953)^2(0.0025)+(-0.6325)^2(0.0016)+(3.1623)^2(0.0009) \approx0.000391+0.000640+0.009000=0.01003. \] Thus \(SE(\widehat r_g)=\sqrt{0.01003}=0.1001\). The approximate 95% Wald interval is \[ 0.3162\pm1.96(0.1001)=0.3162\pm0.1962=(0.1200,\,0.5124). \]

Under the model and approximation used, this interval gives a range of values reasonably compatible with the estimate and its sampling uncertainty. Its nominal 95% coverage refers to the long-run behavior of intervals constructed by this method under its assumptions; it does not mean there is a 95% probability that this fixed interval contains the fixed parameter.

Profile Likelihood and Bootstrap Alternatives

A symmetric Wald interval is convenient, but it may be a poor summary when the likelihood is asymmetric or the estimate is near a boundary. Two alternatives retain more of the shape or sampling behavior of the estimator.

Profile likelihood. For a candidate value \(r\) of genetic correlation, maximize the model likelihood over the nuisance parameters, including the component scales, subject to the correlation being \(r\). This gives the profile log likelihood \(\ell_p(r)\). Let \(\widehat r_g\) maximize it. Under regular likelihood conditions, an approximate 95% profile-likelihood interval consists of values satisfying \(2\{\ell_p(\widehat r_g)-\ell_p(r)\}\leq3.84\), using the 95th percentile of a chi-square distribution with one degree of freedom. This is an asymptotic likelihood interval; its calibration can be poor near boundaries or with limited information. The optimization must respect the model's valid covariance constraints.

Bootstrap. Repeatedly generate or resample data under a procedure that represents the study's sampling structure, refit the same model, and recalculate \(\widehat r_g\) for each replicate. A percentile interval uses the empirical quantiles of the replicate estimates. Resampling must preserve relevant dependence: for example, individuals cannot be treated as independent sampling units when the design requires resampling families or clusters. A parametric bootstrap instead simulates data from the fitted model and is conditional on that model being an adequate description of the data. Replicates that fail to fit or produce boundary estimates should be examined and handled transparently, not silently discarded.

Worked Example: Reading a Profile-Likelihood Interval

In a hypothetical fitted model, the profile log likelihood is maximized at \(\widehat r_g=0.32\), with maximum value \(-418.60\). The values \(r=0.07\) and \(r=0.55\) each give profile log likelihood \(-420.52\).

For either endpoint, the likelihood-ratio distance from the maximum is \[ 2\{(-418.60)-(-420.52)\}=2(1.92)=3.84. \] Thus the approximate 95% profile-likelihood interval is \((0.07,\,0.55)\), using the one-degree-of-freedom chi-square cutoff. It is asymmetric around \(0.32\): the lower distance is \(0.25\), and the upper distance is \(0.23\). The interval reflects the profile likelihood's shape rather than forcing symmetry. The endpoint calculation is meaningful only if the likelihood optimization and its asymptotic calibration are appropriate for the fitted model.

Worked Example: A Bootstrap Percentile Interval

Suppose a hypothetical analysis produces \(1{,}000\) valid bootstrap estimates by resampling independent study units while preserving the paired trait measurements within each unit and refitting the same bivariate model. The original-data estimate is \(0.316\). After sorting the bootstrap estimates, the empirical 2.5th and 97.5th percentiles are \(0.08\) and \(0.52\).

The percentile interval is therefore \((0.08,\,0.52)\). It is formed directly from the bootstrap distribution, rather than by adding and subtracting a standard error. Its interpretation relies on the bootstrap procedure representing the relevant sampling process and on the model being refitted consistently. If units are related or clustered, the resampling scheme must account for that dependence; resampling individual records independently would generally give an inappropriate uncertainty estimate.

Choosing and Reporting an Uncertainty Estimate

The delta method is efficient when component estimates are regular, their joint covariance matrix is available, and the correlation estimate is comfortably away from problematic boundaries. Profile likelihood can better reflect asymmetry when the likelihood is informative and model fitting under fixed correlations is feasible. Bootstrap intervals can represent complicated estimator behavior, but only when the resampling or simulation scheme matches the design and the model is refitted reliably. No method repairs weak identifiability or misspecified data structure.

Report the point estimate, interval method, interval endpoints, and enough model context to identify what genetic components were used. State that the genetic correlation is dimensionless; unlike the covariance, it has no trait-unit product. An interval is an uncertainty summary for this parameter under the stated model, not a measure of the phenotypic correlation and not evidence of causation.

Common pitfall: Treating an estimated genetic covariance as though it were already a correlation, or calculating a standard error from the covariance estimate alone. The correction is to standardize by both estimated genetic standard deviations and propagate uncertainty in all three components, including their covariances. A second warning is to report a symmetric Wald interval automatically: near a boundary or with an unstable variance estimate, inspect whether profile likelihood or a design-aware bootstrap is more defensible.

Key takeaway: Genetic correlation is estimated covariance divided by the product of genetic standard deviations. Its uncertainty depends jointly on both variance estimates and the covariance estimate; choose an interval method whose approximation and resampling assumptions fit the model and study design.

Check Your Understanding

Use the definitions and methods above to answer these questions.

If the estimated genetic covariance is negative while both estimated genetic variances are positive, what is the sign of the plug-in genetic correlation?
Why must the delta-method variance calculation account for covariance among the three estimated genetic components when those estimates are correlated?
What assumption about the component covariance estimates was made in the delta-method worked example, and why is it not a general rule?
What feature of a profile likelihood can make its interval different from a symmetric Wald interval?
Why might resampling individual records independently produce an invalid bootstrap interval for a study with related or clustered participants?