Multivariate Regression for Joint Molecular Phenotypes

Why Model Molecular Phenotypes Together?

A genetic variant may be associated with several molecular measurements from the same people: for example, gene expression, protein abundance, and a metabolite concentration. Analyzing these outcomes together makes the shared study design explicit. It also gives a direct way to ask whether the estimated genotype association is larger for one outcome than another. The model does not, by itself, establish that one molecular layer causes a change in another.

The essential feature is that each individual contributes a vector of outcomes, not three independent records. The outcomes can have different units and can be correlated within an individual. The regression coefficients describe associations conditional on the included covariates; the residual covariance describes how the outcomes vary together after that adjustment.

Definition: For individual \(i\), let \(Y_i=(Y_{i1},\ldots,Y_{iq})\) be a vector of \(q\) molecular phenotypes, \(G_i\) a genotype dosage, and \(C_i\) a vector of covariates. A multivariate linear regression models the conditional mean of the whole outcome vector as a linear function of the predictors, while allowing the outcome residuals for an individual to have a covariance matrix \(\Sigma\).

Model and Assumptions

For three aligned molecular outcomes, such as expression, protein, and metabolite measurements, a useful model is

Y_i = \alpha + G_i\beta + C_i^{\mathsf T}\Gamma + \varepsilon_i, \qquad i=1,\ldots,n.

Here \(Y_i\), \(\alpha\), \(\beta\), and \(\varepsilon_i\) are three-component row vectors; \(\beta=(\beta_1,\beta_2,\beta_3)\) contains the genotype coefficients for the three outcomes. The covariate coefficient matrix \(\Gamma\) allows each covariate to have a different association with each outcome. With an intercept included, the conditional mean is \(E(Y_i\mid G_i,C_i)=\alpha+G_i\beta+C_i^{\mathsf T}\Gamma\).

Usually \(G_i\) is coded as the count of copies of a specified allele, taking values 0, 1, or 2. Then \(\beta_j\) is the expected difference in outcome \(j\) associated with one additional allele, holding the included covariates fixed. A dosage from imputation can instead be a value between 0 and 2, with the same per-allele interpretation under the linear model. These are associations, not automatically causal effects.

For the standard linear-model analysis, assume that the conditional mean is correctly specified, that individuals are independent, and that the predictors have adequate variation and no exact linear dependencies. The residual vectors have mean zero and a common covariance matrix \(\Sigma\), which is positive definite. For exact small-sample \(t\)- and \(F\)-based inference, additionally assume the residual vectors are multivariate normal conditional on the predictors. The outcomes need not be independent within an individual: the off-diagonal elements of \(\Sigma\) represent their residual covariance.

The same individuals and the same design matrix should be used for every outcome in this formulation. If a phenotype is missing for some individuals, fitting each outcome on a different subset no longer gives one common complete multivariate regression. The analysis must then be framed to account for the missingness; simply treating absent measurements as observed values is not valid.

Fitting the Regression

Let \(X\) be the \(n\)-by-\(p\) design matrix containing an intercept, genotype, and the chosen covariates, and let \(Y\) be the \(n\)-by-\(q\) matrix whose rows are the individuals’ outcome vectors. The coefficient matrix \(B\) has one column of regression coefficients per outcome. The fitted model is \(Y=XB+E\), where \(E\) is the residual matrix. Least squares chooses \(B\) to minimize the sum of squared residuals over all outcomes.

\widehat B=(X^{\mathsf T}X)^{-1}X^{\mathsf T}Y.

This formula applies when \(X\) has full column rank. It shows an important practical fact: with the same design matrix for all outcomes, the multivariate least-squares coefficient for each outcome is its usual single-outcome least-squares coefficient. Fitting the outcomes together does not change those point estimates. It does, however, preserve the association between outcome residuals, which matters when making comparisons across outcome coefficients.

Derivation: The least-squares criterion is \(L(B)=\operatorname{tr}\{(Y-XB)^{\mathsf T}(Y-XB)\}\), the sum of squared residuals across individuals and outcomes. Differentiating with respect to \(B\) gives \(-2X^{\mathsf T}(Y-XB)\). Setting this derivative to zero yields \(X^{\mathsf T}XB=X^{\mathsf T}Y\). Since \(X^{\mathsf T}X\) is invertible under the full-rank assumption, solving gives \(\widehat B=(X^{\mathsf T}X)^{-1}X^{\mathsf T}Y\).

Worked Example: Estimating Three Genotype Slopes

Suppose a hypothetical study measures expression, protein, and metabolite outcomes on the same 10 individuals. For a regression with an intercept and genotype dosage, genotype is centered at its sample mean. The centered genotype sum of squares is \(S_{GG}=5\), and the cross-products of centered genotype with each outcome are \(1.5\), \(1.0\), and \(-0.5\), respectively. These are sufficient summaries for the three slopes in this model.

With a centered genotype and an intercept, the slope for outcome \(j\) is the cross-product divided by \(S_{GG}\). Thus the estimates are

\widehat\beta_{\mathrm{expression}}=\frac{1.5}{5}=0.30,\qquad \widehat\beta_{\mathrm{protein}}=\frac{1.0}{5}=0.20,\qquad \widehat\beta_{\mathrm{metabolite}}=\frac{-0.5}{5}=-0.10.

If expression was measured in normalized expression units, protein in abundance units, and metabolite in concentration units, the respective interpretations are those units per additional allele. The negative metabolite slope describes a lower fitted metabolite concentration per additional allele in this hypothetical sample. The numerical sizes cannot be ranked as if they shared a scale.

Covariate Adjustment and Outcome Scale

Covariates can include variables such as age, sex, study batch, and ancestry principal components when appropriate to the study design. They should be selected and coded before interpreting the genotype coefficients. If all outcomes use the same individuals and covariates, adjustment has a convenient equivalent description: remove the linear associations of genotype and each outcome with the covariates, then regress each residualized outcome on residualized genotype. This is the partial-regression form of the same least-squares coefficient.

Comparing coefficients across molecular layers also requires attention to units. A difference such as “0.30 expression units versus 0.20 protein units” is not a meaningful numerical comparison if the units differ. One option is to express each phenotype in standard-deviation units, using a clearly defined scale and the same scaling convention throughout. Then a coefficient is an outcome standard deviation per allele. Scaling makes a numerical contrast interpretable on a common scale, but it does not make the biological measurements interchangeable.

Worked Example: Genotype Associations Adjusted for Covariates

Consider a hypothetical complete-case analysis with 120 individuals, three molecular outcomes, genotype, an intercept, and covariates. After linearly removing the intercept and covariates from genotype and each outcome, suppose the residualized genotype sum of squares is 48. The cross-products between residualized genotype and residualized expression, protein, and metabolite are 7.2, 3.6, and \(-2.4\). Assume the design matrix is full rank and the linear conditional-mean model is suitable.

The adjusted genotype slopes are the residualized cross-products divided by the residualized genotype sum of squares:

\widehat\beta_{\mathrm{expression}}=\frac{7.2}{48}=0.15,\qquad \widehat\beta_{\mathrm{protein}}=\frac{3.6}{48}=0.075,\qquad \widehat\beta_{\mathrm{metabolite}}=\frac{-2.4}{48}=-0.050.

For example, the expression estimate is an adjusted increase of 0.15 expression units per additional allele, conditional on the specified covariates. This calculation describes the fitted slopes; it does not show whether any slope differs statistically from zero. That requires an uncertainty estimate, and the conclusion depends on the model assumptions and the residual variation.

Testing Whether Two Layer-Specific Associations Differ

A direct comparison is a test of a coefficient contrast. For expression and protein, the null hypothesis of equal genotype associations is \(H_0:\beta_{\mathrm{expression}}-\beta_{\mathrm{protein}}=0\), against a specified alternative such as a two-sided difference. The contrast must be defined on scales that make its value meaningful. For example, comparing standardized outcomes makes the difference a contrast in standard deviations per allele.

The uncertainty of this difference depends on the covariance between the two fitted slopes. Because both outcomes were measured on the same individuals, their residuals may be correlated. For a model with an intercept and genotype only, let \(S_{GG}=\sum_i(G_i-\bar G)^2\), and let \(\widehat\Sigma\) be the estimated residual covariance matrix, calculated as \(E^{\mathsf T}E/(n-2)\). The estimated variance of the difference between two slopes is

\widehat{\operatorname{Var}}(\widehat\beta_j-\widehat\beta_k) = \frac{\widehat\Sigma_{jj}+\widehat\Sigma_{kk}-2\widehat\Sigma_{jk}}{S_{GG}}.

The minus-two-covariance term follows from the variance rule \(\operatorname{Var}(A-B)=\operatorname{Var}(A)+\operatorname{Var}(B)-2\operatorname{Cov}(A,B)\). It is not generally correct to treat two layer-specific slope estimates as independent. With covariates, the same contrast principle applies, using the genotype coefficient’s appropriate design-matrix variance factor in place of \(1/S_{GG}\).

Worked Example: Testing an Expression–Protein Difference

Suppose a hypothetical study has 12 individuals, with all three outcomes standardized before fitting an intercept-and-genotype model. The estimated expression and protein slopes are 0.30 and 0.20 standard deviations per allele. The centered genotype sum of squares is 6. The estimated residual variances are 0.36 for expression and 0.25 for protein, and their estimated residual covariance is 0.12. Assume independent individuals, a correctly specified linear mean, and multivariate normal residuals for this small-sample \(t\) test.

The estimated difference is \(0.30-0.20=0.10\) standard deviations per allele. Its estimated variance and standard error are

\widehat{\operatorname{Var}}(\widehat\beta_{\mathrm{expression}}-\widehat\beta_{\mathrm{protein}}) = \frac{0.36+0.25-2(0.12)}{6} = \frac{0.37}{6} = 0.06167, \qquad SE=\sqrt{0.06167}=0.2483.

The test statistic is \(t=0.10/0.2483=0.403\), with \(df=12-2=10\). The two-sided \(p\)-value is approximately 0.696. A 95% confidence interval is \(0.10\pm2.228(0.2483)\), or approximately \((-0.453,0.653)\) standard deviations per allele. The data do not provide evidence of a difference at the 0.05 significance level. This is not proof that the associations are equal: the interval also reflects substantial uncertainty about their difference.

What the Joint Model Does—and Does Not—Say

A multivariate regression is useful when the scientific question concerns several outcomes measured on the same individuals. It gives a consistent set of layer-specific estimates under a shared design, and it retains the residual covariance needed for valid comparisons. A targeted contrast can ask whether two specified genotype associations differ. The answer is about the model’s conditional associations, not whether genotype acts through one molecular layer to affect another.

A contrast between two outcomes is only one possible question. Broader tests involving multiple coefficients or several outcomes require their own clearly stated null hypotheses and appropriate accounting for the covariance structure. In all cases, report the outcome scales, the covariates, the genotype coding, and the coefficient units so readers can tell what a reported difference represents.

Common Pitfalls

Treating the same people’s outcomes as independent. Separate regressions can give the same least-squares point estimates, but an across-layer coefficient comparison must account for residual covariance. Ignoring it can give an incorrect standard error.
Comparing numbers measured in incompatible units. A larger numerical slope does not necessarily mean a stronger association when expression and protein are on different scales. Define a common scale if a numerical difference is the target.
Interpreting a coefficient as causal. A genotype coefficient is a conditional association under the model. Covariate adjustment alone does not establish a causal effect or a pathway from one molecular outcome to another.
Using different analysis samples without noticing. The shared-response formulation assumes aligned outcome rows. If the outcomes are observed on different sets of people, the resulting estimates do not come from one common complete multivariate fit.

Key takeaway: Fit one shared-design regression to the vector of aligned molecular outcomes. Interpret each genotype coefficient in its own outcome units, and use the residual covariance—not an independence assumption—when testing a difference between layer-specific associations.

Check Your Understanding

Answer each question by referring to the model and assumptions in this tutorial.

In a multivariate regression, what does the genotype coefficient for the protein outcome represent when genotype is coded as allele dosage?
Why are the least-squares coefficient estimates the same as fitting each outcome separately when every outcome uses the same design matrix?
Why does a test of the difference between expression and protein slopes involve their residual covariance?
What additional interpretation problem arises when the expression and protein outcomes are measured in different units?
Does finding different genotype coefficients for two molecular outcomes establish that one outcome mediates the genotype association with the other? Explain.