Covariance-Aware Tests Across Omics Layers

Why Test Molecular Outcomes Jointly?

A genetic variant may be associated with several molecular outcomes measured in the same individuals, such as gene expression and protein abundance. Testing each outcome separately answers whether there is evidence for an association with that particular outcome. A joint test instead asks whether the vector of genetic effects across the outcomes is zero. Because the outcomes are correlated, their effect estimates are correlated too; a joint test must account for that covariance rather than treating the estimates as independent.

The joint test is an omnibus test: rejection indicates evidence that at least one component of the effect vector differs from zero. It does not by itself identify which outcome is associated, establish a direction of effect, or show that one molecular outcome causes another.

Model and Joint Hypothesis

Suppose \(n\) individuals have measurements on \(K\) molecular outcomes. Let \(Y_i\) be the \(K\)-dimensional outcome vector for individual \(i\), \(G_i\) the genotype dosage for a variant, and \(C_i\) a vector of covariates. A dosage might be coded as 0, 1, or 2 copies of an effect allele. Consider the multivariate linear model

Y_i = \alpha + \beta_G G_i + \Gamma C_i + \varepsilon_i,\qquad \varepsilon_i \sim N_K(0,\Sigma).

Here \(\beta_G\) is a \(K\)-component vector: its \(k\)th element is the expected change in outcome \(k\) per additional effect allele, conditional on the covariates. The matrix \(\Sigma\) describes residual variances and covariances among outcomes after accounting for genotype and covariates. For example, a positive off-diagonal element means that individuals with higher-than-predicted values for one outcome tend also to have higher-than-predicted values for another.

Definition: The joint null hypothesis for one variant across \(K\) outcomes is \(H_0:\beta_G=0\), meaning that all \(K\) genotype effects are zero. The alternative is \(H_a:\beta_G\ne 0\), meaning that at least one component is nonzero.

The standard multivariate regression assumptions are that the conditional mean is correctly represented by the model, observations are independent across individuals, the outcomes have a common residual covariance matrix \(\Sigma\), and the design matrix has enough independent information to estimate the coefficients. For the Gaussian likelihood-ratio test below, the residual vectors are additionally assumed to be multivariate normal. Using the same individuals and covariates for all outcomes gives a common design and makes the covariance-aware calculations especially direct. If observations are missing in different outcome layers, these simple formulas do not automatically apply.

The Joint Wald Test

Let \(\hat{\beta}_G\) be the estimated vector of genetic effects, and let \(\widehat V\) be its estimated \(K\) by \(K\) covariance matrix. Its diagonal entries are the estimated variances of the outcome-specific effect estimates; its off-diagonal entries quantify their estimated covariance. The Wald statistic measures the estimated effect vector relative to its uncertainty, using the full covariance matrix.

Definition: The joint Wald statistic is \(W=\hat{\beta}_G^{\,T}\widehat V^{-1}\hat{\beta}_G\). Under \(H_0\), with a correctly specified model and standard large-sample regularity conditions, \(W\) has an approximate chi-square distribution with \(K\) degrees of freedom.

For a common design, the structure of \(\widehat V\) follows from the residual covariance matrix. After adjusting for covariates, let \(g_c\) be the genotype vector residualized against those covariates, and let \(S_G=g_c^Tg_c\). If \(\widehat\Sigma\) estimates the residual covariance across outcomes, then the estimated covariance matrix of the genetic effect estimates is \(\widehat V=\widehat\Sigma/S_G\). Thus the Wald statistic can also be written as \(S_G\hat{\beta}_G^{\,T}\widehat\Sigma^{-1}\hat{\beta}_G\). In practice, it is important to use the covariance of the effect estimates, not merely the outcome variances while setting their covariance to zero.

Worked Example: Correlated Expression and Protein Effects

Suppose a hypothetical study estimates a per-allele expression effect of \(0.12\) standard deviations and a protein effect of \(0.12\) standard deviations. The estimated covariance matrix of these two effect estimates, in squared standard-deviation-per-allele units, is

\widehat V = \begin{pmatrix} 0.010 & 0.008\\ 0.008 & 0.010 \end{pmatrix}, \qquad \hat{\beta}_G = \begin{pmatrix} 0.12\\ 0.12 \end{pmatrix}.

The correlation between the effect estimates is \(0.008/\sqrt{0.010(0.010)}=0.80\). The eigenvalue of \(\widehat V\) in the common-effect direction \((1,1)^T\) is \(0.018\). Since \(\hat{\beta}_G\) lies in that direction, its squared length is \(0.12^2+0.12^2=0.0288\), giving

W = \frac{0.0288}{0.018}=1.600.

With \(K=2\), the reference distribution is \(\chi^2_2\), so the approximate p-value is \(P(\chi^2_2\ge 1.600)=\exp(-1.600/2)=0.4493\). The joint test does not provide strong evidence against the null. Each separate univariate Wald statistic is \(0.12^2/0.010=1.440\), with a two-sided standard-normal p-value of about \(0.230\). The joint result accounts for the fact that these similar-direction estimates are highly correlated; it does not count them as two independent pieces of evidence.

Worked Example: Opposite Effects Along a Low-Variance Direction

Now suppose the same hypothetical covariance matrix applies, but the estimated expression effect is \(0.12\) standard deviations per allele and the protein effect is \(-0.12\). The separate univariate statistics are still \(1.440\) each. The effect vector now lies in the contrast direction \((1,-1)^T\), whose covariance eigenvalue is \(0.010-0.008=0.002\). Therefore,

W = \frac{0.12^2+(-0.12)^2}{0.002} =\frac{0.0288}{0.002} =14.40.

For two degrees of freedom, the approximate p-value is \(\exp(-14.40/2)=0.0007466\). The joint test gives strong evidence that the vector of effects is not zero, even though neither separate univariate test is conventionally significant. This occurs because, under the estimated positive covariance, an expression effect and protein effect of opposite signs point in a direction with relatively little sampling variation. The result is evidence about the joint association pattern, not evidence that the variant causes one outcome to change the other.

The Multivariate Likelihood-Ratio Test

A second approach compares the fit of two nested Gaussian multivariate regression models. The reduced model includes the covariates but omits the genotype. The full model includes both the covariates and the genotype, with a separate genetic effect for each outcome. Let \(E_0\) be the residual sum-of-products matrix from the reduced model and \(E_1\) the corresponding matrix from the full model. Both are \(K\) by \(K\); their diagonal entries are residual sums of squares, and their off-diagonal entries are residual cross-products.

Definition: With the residual covariance estimated freely in each model, the likelihood-ratio statistic is \(LR=n\log\bigl(|E_0|/|E_1|\bigr)\), where \(|E|\) denotes a matrix determinant. Under the joint null and standard regularity conditions, \(LR\) is approximately chi-square with \(K\) degrees of freedom for one genotype term tested across \(K\) outcomes.

The determinants measure the overall residual variation across outcomes, including their covariance. A smaller determinant for the full model indicates an improvement in joint fit. This statistic is not generally the sum of the \(K\) separate univariate likelihood-ratio statistics.

Why the Joint Wald and Likelihood-Ratio Tests Are Related

For one added genotype term, the full model reduces the residual sum-of-products matrix by a rank-one matrix. If \(\hat{\beta}_G\) is the full-model effect estimate and \(S_G\) is the residualized genotype sum of squares, then \(E_1=E_0-S_G\hat{\beta}_G\hat{\beta}_G^T\). The matrix determinant lemma, \(|A+uv^T|=|A|(1+v^TA^{-1}u)\), gives a compact expression for the likelihood change.

t=S_G\hat{\beta}_G^TE_0^{-1}\hat{\beta}_G,\qquad \frac{|E_1|}{|E_0|}=1-t,\qquad LR=-n\log(1-t).

The determinant ratio follows by setting \(A=E_0\), \(u=-S_G\hat{\beta}_G\), and \(v=\hat{\beta}_G\). Since the full-model residual matrix is positive definite in the usual nonsingular case, \(1-t\) is positive. For small \(t\), the expansion \(-\log(1-t)\approx t\) shows that \(LR\approx nt\). If the full-model maximum-likelihood estimate of the residual covariance is \(\widehat\Sigma=E_1/n\), the corresponding Wald statistic is \(W=S_G\hat{\beta}_G^T\widehat\Sigma^{-1}\hat{\beta}_G=nt/(1-t)\). Thus both statistics use the same fitted effects and covariance information, and are close when the evidence is modest; they need not be numerically identical. In large samples, both are compared with the same chi-square reference distribution under the null.

Worked Example: Comparing Joint and Separate Likelihood-Ratio Tests

Consider hypothetical expression and protein data from \(n=100\) individuals, with the same covariates and complete outcome measurements for everyone. Suppose the reduced-model and full-model residual sum-of-products matrices are

E_0= \begin{pmatrix} 100&60\\ 60&100 \end{pmatrix}, \qquad E_1= \begin{pmatrix} 96&56\\ 56&96 \end{pmatrix}.

The determinants are \(|E_0|=100(100)-60(60)=6400\) and \(|E_1|=96(96)-56(56)=6080\). The joint likelihood-ratio statistic is therefore

LR=100\log(6400/6080) =100\log(1.05263) =5.129.

Against \(\chi^2_2\), the approximate p-value is \(\exp(-5.129/2)=0.0769\). At a significance level of \(0.05\), this joint test does not reject the null. Separately, each outcome's residual sum of squares falls from \(100\) to \(96\), producing a univariate likelihood-ratio statistic \(100\log(100/96)=4.082\). Each has an approximate \(\chi^2_1\) p-value of about \(0.043\). These marginal results do not change the joint conclusion: separate tests answer different questions, and the joint determinant reflects the covariance between the outcomes. The p-values are approximate because the chi-square reference distributions are large-sample results.

Choosing and Interpreting the Test

The joint Wald test is often convenient when the effect estimates and their covariance matrix are already available. The likelihood-ratio test is natural when both nested multivariate models can be fitted and compared. Under the stated model and large-sample conditions, the tests target the same joint null, although their finite-sample values can differ. Neither test removes the need to consider how many variants, outcome sets, or other hypotheses were tested; the stated p-value is for this particular joint test and does not itself adjust for a larger testing family.

A joint rejection says that the data are inconsistent with all genetic effects being zero, under the model and test assumptions. It does not establish that every outcome has a nonzero effect. To identify which outcomes contribute, inspect the estimated effects and uncertainty and use appropriately planned follow-up tests. The covariance describes residual co-variation after model adjustment; it is not, by itself, evidence of a causal relationship between molecular layers.

Common Mistake

Wrong reasoning: “The expression and protein estimates are correlated, so I can add their two univariate test statistics as if they were independent.”
Correction: The sampling covariance of the effect estimates changes the uncertainty in different directions of the effect vector. Use the full covariance matrix in the Wald statistic, or use the joint likelihood from the multivariate model. Ignoring the off-diagonal terms can either overstate or understate the evidence.
Wrong reasoning: “A significant joint test means the variant affects both outcomes.”
Correction: The joint alternative requires only one or more components to be nonzero. Outcome-specific estimates and follow-up inference are needed to describe which effects are supported.

Key takeaway: A covariance-aware joint test evaluates the entire genetic effect vector against its estimated uncertainty. Correlation between molecular outcomes is part of the test, not a reason to treat separate results as independent.

Check Your Understanding

Answer these questions using the model and test definitions above.

For three molecular outcomes, how many degrees of freedom does the large-sample joint Wald test use when testing one variant’s effects?
What does an off-diagonal entry of the covariance matrix of the effect estimates represent?
Why can two individually nonsignificant univariate tests still yield a significant joint Wald test?
In the likelihood-ratio test, what do \(E_0\) and \(E_1\) represent, and why are their determinants used?
What conclusion is justified by rejecting the joint null, and what conclusion is not justified without further analysis?