Multivariable Mendelian Randomization

Why Use Several Exposures in One MR Model?

A genetic variant can be associated with more than one exposure. If two exposures both relate to an outcome, analyzing each exposure separately may attribute to one exposure part of the effect carried through the other. Multivariable Mendelian randomization (MVMR) represents the variant–outcome association as a joint function of the variant associations with all included exposures. The resulting coefficient for one exposure is a conditional effect: its effect when the other modeled exposures are held fixed.

This interpretation is not a guarantee that the estimate is unbiased. It depends on causal assumptions, adequate information in the variant–exposure associations, and a correctly specified model. The purpose here is to define that model, derive its summary-data estimator, and distinguish its estimates from separate univariable MR estimates.

The Multivariable Summary-Data Model

Definition: Suppose $J$ genetic variants are used as instruments for $K$ exposures. Let $\gamma_{jk}$ be the association of variant $j$ with exposure $k$, and let $\Gamma_j$ be its association with the outcome. Let $\theta_k$ be the causal effect of exposure $k$ on the outcome, conditional on the other included exposures. The multivariable summary-data model is $$ \Gamma_j = \sum_{k=1}^{K}\gamma_{jk}\theta_k+\alpha_j+\varepsilon_j, \qquad j=1,\ldots,J. $$ Here, $\alpha_j$ is the variant's effect on the outcome through pathways not represented by the included exposures, and $\varepsilon_j$ is sampling error in the estimated variant–outcome association.

In the ideal model, $\alpha_j=0$ for every instrument. This is the joint exclusion restriction: conditional on the exposures in the model, the variants have no additional pathway to the outcome. In practice, the $\gamma_{jk}$ and $\Gamma_j$ are estimated associations, commonly taken from genome-wide association study summary data. They must refer to compatible alleles, exposure and outcome scales, and populations; allele harmonization is needed so that the associations have consistent directions.

For causal interpretation, the instruments must also be independent of confounders of the exposures and outcome, and the modeled exposures must be sufficiently well defined for their effects to have a meaningful interpretation. The model assumes a linear, additive relationship on the scales used, with no exposure interactions represented. If interactions or nonlinear effects are important, a single vector of coefficients may not describe the causal relationship. The instrument set must carry distinct information about the included exposures; having many variants alone does not ensure this.

Write $\boldsymbol{\Gamma}$ for the $J$-vector of variant–outcome associations and $\mathbf{G}$ for the $J$-by-$K$ matrix whose $j,k$ entry is $\gamma_{jk}$. Write $\boldsymbol{\theta}$ for the $K$-vector of causal effects. Under the joint exclusion restriction and the summary-data model,

\boldsymbol{\Gamma}=\mathbf{G}\boldsymbol{\theta}+\boldsymbol{\varepsilon}.

If an exposure is measured in one unit and the outcome in another, the corresponding $\theta_k$ is in outcome units per unit of that exposure, conditional on the remaining included exposures. For example, a coefficient could be change in outcome units per one-unit increase in exposure A, holding exposure B fixed. It is not the coefficient from a model containing exposure A alone.

Weighted Estimation and Its Derivation

Assume first that the variants are independent and that sampling errors in their outcome-association estimates are independent, with known or estimated variances $s_{\Gamma_j}^{\,2}$. Define $\mathbf{W}$ as the diagonal matrix with entries $1/s_{\Gamma_j}^{\,2}$. A weighted least-squares estimate minimizes the sum of squared outcome-association residuals, weighted by their precision. More generally, if the vector of errors has covariance matrix $\boldsymbol{\Sigma}$, use $\mathbf{W}=\boldsymbol{\Sigma}^{-1}$. A non-diagonal covariance may be necessary when variants are correlated through linkage disequilibrium or when sampling errors are correlated.

Estimator: If $\mathbf{G}$ has full column rank and $\mathbf{W}$ is positive definite, the weighted least-squares estimate is $$ \hat{\boldsymbol{\theta}} = \left(\mathbf{G}^{T}\mathbf{W}\mathbf{G}\right)^{-1} \mathbf{G}^{T}\mathbf{W}\boldsymbol{\Gamma}. $$

The full-column-rank condition means that no exposure column of $\mathbf{G}$ is an exact linear combination of the other columns. Without it, the data cannot distinguish all $K$ coefficients. In usual summary-data estimation, uncertainty in the exposure associations $\gamma_{jk}$ is not fully incorporated into this weighted least-squares formula. Treating those associations as fixed is an approximation, most defensible when their estimation error is small relative to the information they contain. The formula alone does not establish that this approximation is adequate.

The estimator follows directly from minimizing the weighted residual sum of squares. With $\mathbf{r}(\boldsymbol{\theta})=\boldsymbol{\Gamma}-\mathbf{G}\boldsymbol{\theta}$, the objective is $Q(\boldsymbol{\theta})=\mathbf{r}(\boldsymbol{\theta})^T\mathbf{W}\mathbf{r}(\boldsymbol{\theta})$. Differentiating and setting the gradient to zero gives the normal equations:

\frac{\partial Q}{\partial\boldsymbol{\theta}} =-2\mathbf{G}^{T}\mathbf{W} (\boldsymbol{\Gamma}-\mathbf{G}\boldsymbol{\theta})=\mathbf{0}, \qquad \mathbf{G}^{T}\mathbf{W}\mathbf{G}\boldsymbol{\theta} =\mathbf{G}^{T}\mathbf{W}\boldsymbol{\Gamma}.

When $\mathbf{G}^{T}\mathbf{W}\mathbf{G}$ is invertible, multiplying the normal equations by its inverse yields the estimator above. This derivation explains both the weighted fit and why full column rank is required. It does not prove the causal assumptions: a precise fitted coefficient can still be biased if the joint exclusion restriction or other causal assumptions fail.

Worked Example: Calculating Two Conditional Effects

Suppose a hypothetical analysis includes two exposures, A and B, and has already harmonized its variant associations. For illustration, its weighted summary-data cross-products are $$ \mathbf{G}^{T}\mathbf{W}\mathbf{G} = \begin{pmatrix} 0.06 & 0.03\\ 0.03 & 0.03 \end{pmatrix}, \qquad \mathbf{G}^{T}\mathbf{W}\boldsymbol{\Gamma} = \begin{pmatrix} 0.024\\ 0.010 \end{pmatrix}. $$ These are sufficient statistics for the weighted estimate; the cross-product matrix is invertible, so the two coefficients can be estimated.

The determinant of the cross-product matrix is $0.06(0.03)-0.03(0.03)=0.0009$. Its inverse is therefore $$ \left(\mathbf{G}^{T}\mathbf{W}\mathbf{G}\right)^{-1} = \begin{pmatrix} 33.333 & -33.333\\ -33.333 & 66.667 \end{pmatrix}. $$ Multiplying by the cross-product vector gives $$ \hat{\boldsymbol{\theta}} = \begin{pmatrix} 33.333 & -33.333\\ -33.333 & 66.667 \end{pmatrix} \begin{pmatrix} 0.024\\ 0.010 \end{pmatrix} = \begin{pmatrix} 0.4667\\ -0.1333 \end{pmatrix}. $$

The first component is the estimated effect of A, conditional on B, in outcome units per unit of A. The second is the estimated conditional effect of B: $-0.1333$ outcome units per unit of B, holding A fixed. The negative sign follows from the second row calculation, $-33.333(0.024)+66.667(0.010)\approx-0.8000+0.6667=-0.1333$. These estimates are conditional coefficients; their signs and magnitudes do not by themselves validate the assumptions.

Why Conditional and Univariable Estimates Differ

A univariable MR analysis for exposure A uses the variant–A associations but does not jointly model the variants' associations with B. If the variants also affect B and B affects the outcome, the outcome associations can reflect both pathways. A multivariable model attempts to separate the modeled contributions by using each variant's association with every included exposure. Thus, the univariable and multivariable coefficients answer different questions and need not agree, even when both analyses use the same variants.

Worked Example: Separate Estimates Versus Joint Estimates

Consider hypothetical associations for two independent variants. Suppose the true summary relationships in this illustration follow $\Gamma_j=0.5\gamma_{jA}+0.2\gamma_{jB}$, with no additional pathway. The associations are:

Variant	$\gamma_{jA}$	$\gamma_{jB}$	$\Gamma_j$
1	0.10	0.08	0.066
2	0.04	0.10	0.040

For the joint model, the equations are $0.10\theta_A+0.08\theta_B=0.066$ and $0.04\theta_A+0.10\theta_B=0.040$. The determinant is $0.10(0.10)-0.08(0.04)=0.0068$. Solving, $$ \hat{\theta}_A =\frac{0.066(0.10)-0.08(0.040)}{0.0068} =\frac{0.0034}{0.0068}=0.5, \qquad \hat{\theta}_B =\frac{0.10(0.040)-0.04(0.066)}{0.0068} =\frac{0.00136}{0.0068}=0.2. $$

For comparison, fit each exposure separately through the origin with equal outcome-association weights. The univariable weighted slope for A is $$ \frac{0.10(0.066)+0.04(0.040)}{0.10^2+0.04^2} =\frac{0.0082}{0.0116}=0.7069. $$ For B it is $$ \frac{0.08(0.066)+0.10(0.040)}{0.08^2+0.10^2} =\frac{0.00928}{0.0164}=0.5659. $$

The separate estimates, 0.7069 for A and 0.5659 for B, differ from the joint estimates of 0.5 and 0.2 because each separate fit leaves the other exposure's contribution in the outcome associations. The example uses an exact joint model, so the difference is not caused by random sampling error. In real data, sampling error and violations of assumptions can add further differences.

Identification, Diagnostics, and Limitations

A central diagnostic is whether the exposure-association matrix provides independent information about the exposures. If two columns are identical or nearly proportional, the model cannot reliably separate their effects. Large standard errors or unstable estimates can result when columns are nearly dependent, even if software returns a numerical answer. Conditional instrument strength is a related issue, but its formal assessment requires tools beyond the estimator derived here.

Worked Example: When Two Effects Cannot Be Separated

Suppose three hypothetical variants have associations $(0.02,0.04)$, $(0.05,0.10)$, and $(0.08,0.16)$ with exposures A and B. In every row, $\gamma_{jB}=2\gamma_{jA}$. The modeled outcome association is therefore $$ \Gamma_j =\gamma_{jA}\theta_A+\gamma_{jB}\theta_B =\gamma_{jA}(\theta_A+2\theta_B). $$

Only the combination $\theta_A+2\theta_B$ is determined by these associations. For example, $(\theta_A,\theta_B)=(0.4,0.1)$ and $(0.2,0.2)$ both give $\theta_A+2\theta_B=0.6$, so both predict the same variant–outcome associations. No increase in the number of variants with this exact proportionality can distinguish the two effects. In matrix terms, the columns are linearly dependent, $\mathbf{G}^{T}\mathbf{W}\mathbf{G}$ is singular, and the separate conditional coefficients are not identified.

Other checks concern whether the data and model are coherent. Verify that exposure and outcome associations refer to aligned effect alleles and comparable scales, that variant dependence is handled in the covariance or by an appropriate selection strategy, and that the weighting reflects the uncertainty structure being assumed. Large residuals can indicate incompatibility with the fitted joint linear relationship, but they do not identify the cause: residual variation may arise from sampling error, model misspecification, or pathways outside the included exposures. Conversely, small residuals do not establish the causal assumptions.

MVMR is not simply a way to “adjust away” every unwanted pathway. Its interpretation is conditional on the exposures actually included. Omitting a relevant exposure can leave a pathway unmodeled; including highly related exposures can make their effects difficult to distinguish. Measurement error in exposure associations, weak conditional information, and violations of the joint exclusion restriction can also undermine estimates. Robust methods and pleiotropy diagnostics address some forms of model failure, but do not remove the need to state the target exposures and assumptions clearly.

Common Mistakes and Practical Interpretation

Calling a multivariable coefficient an unadjusted effect. It is conditional on the other included exposures under the model.
Assuming separate univariable estimates must match. They use different models and can retain contributions from other exposures.
Reading a returned estimate as proof of identification. If exposure-association columns are dependent, distinct effects cannot be separated; near-dependence can make estimates unstable.
Treating the regression equation as the causal argument. The equation produces an estimate, while causal interpretation additionally requires valid instruments, the joint exclusion restriction, and an adequate functional form.
Ignoring allele alignment or covariance. Inconsistent allele directions or unmodeled dependence can invalidate the calculation and its uncertainty assessment.

The key distinction is between what the model estimates and what assumptions allow that estimate to mean. The weighted regression estimates the coefficients that best relate the variant–outcome associations to all included variant–exposure associations. Under the causal assumptions, those coefficients represent conditional effects. Without adequate identification or valid assumptions, they remain model-based associations rather than established causal effects.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In a two-exposure MVMR model, what does the coefficient for exposure A represent?
Suppose two exposure-association columns in $\mathbf{G}$ are exactly proportional. What can and cannot be identified?
For the hypothetical associations in the second worked example, calculate the univariable slope for A and explain why it is not the joint coefficient for A.
Name two assumptions or data conditions required before interpreting the weighted estimate causally.

Solutions. (1) It is the effect of A on the outcome conditional on the other included exposure, under the model's causal assumptions. (2) The two separate effects cannot be identified; only a particular combination of them is informed by the proportional columns. (3) The equal-weight slope is $[0.10(0.066)+0.04(0.040)]/[0.10^2+0.04^2]=0.0082/0.0116=0.7069$. It is a separate fit that does not model B's contribution, unlike the joint coefficient of 0.5. (4) Examples include instrument independence from exposure–outcome confounders, the joint exclusion restriction, adequate rank and information in $\mathbf{G}$, correct allele alignment, and an adequate linear additive model.