Multivariable Mendelian Randomization for Direct Effects

Why Estimate Effects Jointly?

A genetic variant can influence more than one exposure. If a variant influences both an exposure of interest and a mediator, for example, its association with the outcome reflects both pathways. A univariable Mendelian randomization analysis that assigns the entire association to just one exposure can therefore give a misleading effect estimate. Multivariable Mendelian randomization (MVMR) models the genetic associations with multiple exposures together, so it can estimate each exposure’s effect conditional on the others included in the model.

The word “direct” needs care. In MVMR, a direct-effect estimate is a coefficient for one exposure while holding the other modeled exposures fixed. It is not automatically a total effect, a natural direct effect, or proof of a particular biological pathway. Its causal interpretation depends on the instrument and model assumptions.

Definition: Let $Z_j$ denote genetic variant $j$, and let $X_1,\ldots,X_K$ denote $K$ exposures. Write $\gamma_{jk}$ for the association of $Z_j$ with exposure $X_k$, and $\Gamma_j$ for its association with outcome $Y$. In a linear additive model, the summary-association relationship is $$ \Gamma_j = \sum_{k=1}^{K}\gamma_{jk}\beta_k + \alpha_j. $$ Here $\beta_k$ is the effect of exposure $X_k$ conditional on the other modeled exposures, and $\alpha_j$ represents any effect of $Z_j$ on $Y$ not operating through those exposures. If all instruments satisfy the exclusion restriction, then $\alpha_j=0$ for every $j$.

The Joint Estimation Model

For $J$ variants, arrange the exposure associations in a $J\times K$ matrix $G$, with entry $G_{jk}=\gamma_{jk}$. Let $\boldsymbol{\Gamma}$ be the vector of the $J$ SNP–outcome associations and $\boldsymbol{\beta}$ the vector of $K$ exposure effects. Under the exclusion restriction, the model is $\boldsymbol{\Gamma}=G\boldsymbol{\beta}$, apart from sampling error. Variants do not need to affect only one exposure: the model allows a row of $G$ to have several nonzero entries.

A common summary-data estimator is weighted regression through the origin. For independent variants, one may use weights $w_j=1/SE(\Gamma_j)^2$, where $SE(\Gamma_j)$ is the standard error of the SNP–outcome association. The estimator minimizes the weighted squared discrepancies between observed and model-predicted outcome associations:

\hat{\boldsymbol{\beta}} = \arg\min_{\boldsymbol{\beta}} \sum_{j=1}^{J}w_j\left(\Gamma_j-\sum_{k=1}^{K}\gamma_{jk}\beta_k\right)^2 = (G^\mathsf{T}WG)^{-1}G^\mathsf{T}W\boldsymbol{\Gamma},

where $W$ is the diagonal matrix with entries $w_j$. The inverse exists only when the columns of $G$ are linearly independent. With correlated variants, the covariance among their outcome-association estimates must be accounted for in the weighting; the independent-variant expression above is not appropriate without that qualification.

Identification requirements: The included variants must be relevant to the exposures jointly: $G$ must have sufficient rank to distinguish the $K$ exposure effects. In addition, each exposure must be conditionally well predicted by the genetic instruments, rather than being nearly redundant with the other included exposure associations. Rank is a mathematical minimum; weak conditional prediction can still make estimates imprecise and unstable.

Assumptions Behind a Direct-Effect Interpretation

The regression formula produces an estimate, but the formula alone does not make it causal. Interpretation as an effect of an exposure requires assumptions about the variants, the exposures, and the model. In particular:

Relevance: the variants jointly predict the included exposures, with enough conditional strength to distinguish their effects.
Independence: the variants are independent of factors that confound the modeled exposure–outcome relationships. Population structure, for example, can compromise this condition if it induces associations between genetic variants and confounders.
Exclusion: conditional on the included exposures, the variants affect the outcome only through those exposures. A variant’s effects through an included exposure are allowed; a direct effect or a pathway through an omitted exposure violates this assumption.
Model form: the additive linear model is an adequate representation on the scale of the associations being analyzed. Coefficients are conditional on the other included exposures.
Compatible association estimates: SNP–exposure and SNP–outcome associations must be harmonized to the same effect allele and must refer to compatible populations and phenotype definitions. Otherwise, the rows of $G$ and the entries of $\boldsymbol{\Gamma}$ may not describe the same genetic contrasts.

Some assumptions can be assessed only indirectly. A full-rank association matrix does not establish adequate conditional strength, and good fit of the regression does not prove that all exclusion restrictions hold. Robust methods and sensitivity analyses can help explore departures from some assumptions, but do not make invalid instruments valid by definition.

Worked Example: Estimating Two Conditional Effects

Suppose a hypothetical analysis uses three independent variants and two continuous exposures, $X_1$ and $X_2$. All association estimates use the same effect-allele orientation. The SNP–outcome associations have equal standard errors, so the weighted regression reduces to ordinary least squares through the origin.

Variant	$\gamma_{j1}$	$\gamma_{j2}$	$\Gamma_j$
1	0.10	0.00	0.06
2	0.00	0.20	0.04
3	0.10	0.20	0.11

For this example, assume the variants satisfy independence and exclusion, the exposures and outcome follow the stated additive linear model, and their conditional genetic prediction is adequate. These are hypothetical assumptions, not facts established by the table. The design matrix and outcome vector are $$ G= \begin{pmatrix} 0.10&0\\ 0&0.20\\ 0.10&0.20 \end{pmatrix}, \qquad \boldsymbol{\Gamma}= \begin{pmatrix} 0.06\\ 0.04\\ 0.11 \end{pmatrix}. $$ The normal-equation components are $$ G^\mathsf{T}G= \begin{pmatrix} 0.02&0.02\\ 0.02&0.08 \end{pmatrix}, \qquad G^\mathsf{T}\boldsymbol{\Gamma}= \begin{pmatrix} 0.017\\ 0.030 \end{pmatrix}. $$ The determinant of $G^\mathsf{T}G$ is $0.02(0.08)-0.02(0.02)=0.0012$, which is positive. Solving the two normal equations gives $$ \hat{\beta}_1 = \frac{0.017(0.08)-0.02(0.030)}{0.0012} = 0.6333, \qquad \hat{\beta}_2 = \frac{0.02(0.030)-0.02(0.017)}{0.0012} = 0.2167. $$ If $X_1$ is measured in years, $X_2$ in kilograms, and $Y$ in millimeters, these estimates have units millimeters per year and millimeters per kilogram, respectively. Under the stated assumptions, the estimate for $X_1$ describes its outcome effect conditional on $X_2$, and the estimate for $X_2$ describes its effect conditional on $X_1$. They are not effects for changing the exposures together in a particular way.

Why the Joint Model Works

The estimator follows from minimizing the weighted residual sum of squares. Let $Q(\boldsymbol{\beta})=(\boldsymbol{\Gamma}-G\boldsymbol{\beta})^\mathsf{T}W(\boldsymbol{\Gamma}-G\boldsymbol{\beta})$. Differentiating with respect to $\boldsymbol{\beta}$ and setting the derivative equal to zero gives the normal equations:

\frac{\partial Q}{\partial\boldsymbol{\beta}} = -2G^\mathsf{T}W(\boldsymbol{\Gamma}-G\boldsymbol{\beta})=0 \quad\Longrightarrow\quad G^\mathsf{T}WG\boldsymbol{\beta}=G^\mathsf{T}W\boldsymbol{\Gamma}.

If $G^\mathsf{T}WG$ is invertible, multiplying both sides by its inverse yields $\hat{\boldsymbol{\beta}}=(G^\mathsf{T}WG)^{-1}G^\mathsf{T}W\boldsymbol{\Gamma}$. If it is not invertible, the data do not uniquely determine all the included effects. This is an identification problem, not something that can be repaired by simply reporting more decimal places.

Worked Example: When Two Effects Cannot Be Separated

Consider three hypothetical variants whose exposure associations are $(0.10,0.20)$, $(0.20,0.40)$, and $(0.15,0.30)$. The second exposure association is twice the first for every variant, so the columns of $G$ are proportional. Suppose the outcome associations are $0.05$, $0.10$, and $0.075$, respectively, with no sampling error for this illustration.

For each variant, the model requires $\Gamma_j=\gamma_{j1}\beta_1+\gamma_{j2}\beta_2$. The first row therefore gives $$ 0.05=0.10\beta_1+0.20\beta_2 \quad\Longrightarrow\quad \beta_1+2\beta_2=0.5. $$ The other two rows give the same equation, not new information. For example, both $(\beta_1,\beta_2)=(0.5,0)$ and $(0.1,0.2)$ satisfy it. Thus the joint genetic associations identify only the combination $\beta_1+2\beta_2$, not the two effects separately. The matrix has rank one, below the two needed columns. A separate causal interpretation for each exposure is not identified by these data.

Direct Effects in a Mediation Question

In a mediation application, the modeled exposures might be an upstream exposure $X$ and a proposed mediator $M$, with outcome $Y$. The joint model is $\Gamma_j=\gamma_{jX}\beta_X+\gamma_{jM}\beta_M$ when the exclusion restriction holds. Here $\beta_X$ is the effect of $X$ conditional on $M$, while $\beta_M$ is the effect of $M$ conditional on $X$. MVMR can use variants that influence both $X$ and $M$, provided their joint association patterns contain enough independent information and the validity assumptions hold.

Under a suitable linear additive causal model, $\beta_X$ can be interpreted as a direct effect with the mediator held fixed. That interpretation is not automatic: it relies on valid instruments for the included exposures and a correctly specified model. In particular, a genetic effect on $Y$ through an unmodeled pathway can violate exclusion. Also, the conditional coefficient alone does not give the total effect or establish the amount mediated. Those are distinct quantities with their own definitions and assumptions.

Worked Example: Interpreting an Exposure and Mediator Coefficient

Suppose a hypothetical MVMR analysis models $X$, measured in years, and mediator $M$, measured in kilograms, against outcome $Y$, measured in millimeters. For three independent variants, the exposure-association rows are $(0.20,0)$, $(0,0.10)$, and $(0.20,0.10)$. The corresponding outcome associations are $0.12$, $0.03$, and $0.15$. Assume equal outcome-association standard errors, adequate conditional strength, and the independence, exclusion, and linear-model assumptions.

The normal-equation quantities are $$ G^\mathsf{T}G= \begin{pmatrix} 0.08&0.02\\ 0.02&0.02 \end{pmatrix}, \qquad G^\mathsf{T}\boldsymbol{\Gamma} = \begin{pmatrix} 0.054\\ 0.018 \end{pmatrix}. $$ The determinant is $0.08(0.02)-0.02(0.02)=0.0012$, so the effects are algebraically identifiable. Solving gives $$ \hat{\beta}_X = \frac{0.054(0.02)-0.02(0.018)}{0.0012} = 0.60, \qquad \hat{\beta}_M = \frac{0.08(0.018)-0.02(0.054)}{0.0012} = 0.30. $$ With the stated units, the estimates are 0.60 millimeters per year for $X$, conditional on $M$, and 0.30 millimeters per kilogram for $M$, conditional on $X$. If the assumptions support causal interpretation, the first is a direct, conditional effect estimate. It should not be relabeled the total effect of $X$; the joint regression by itself has not calculated that quantity.

Common Mistakes and Practical Checks

Assigning a shared genetic association to only one exposure. A variant associated with both included exposures contributes through both terms in the model. Omitting one term can distort the remaining coefficient.
Assuming full rank guarantees strong instruments. Full rank only says the effects can be distinguished in principle. Near-proportional columns can still make conditional estimates highly uncertain.
Calling every conditional coefficient a natural direct effect. The MVMR coefficient is conditional on modeled exposures. Further assumptions are needed for stronger mediation interpretations.
Treating good regression fit as proof of instrument validity. A fit statistic cannot establish independence or rule out every unmodeled pathway from a variant to the outcome.

In practice, inspect the exposure-association matrix for rank and near-collinearity, assess conditional instrument strength, and state which exposures are included in the model. Interpret each coefficient with its conditioning set and units. Standard regression uncertainty calculations can also be inadequate when SNP–exposure associations are imprecise; weak-instrument methods and sensitivity analyses address related concerns, but do not remove the need to state the assumptions.

Key takeaway: Multivariable Mendelian randomization estimates exposure effects jointly by relating each variant’s outcome association to its associations with all included exposures. A coefficient is a conditional direct-effect estimate only when the genetic instruments, conditional relevance, and model support that causal interpretation.

Check Your Understanding

Try these questions using the definitions and model above.

In the MVMR summary-association model, what does $\alpha_j=0$ assert for variant $j$?
Why can a variant associated with both an exposure and a mediator still be used in MVMR, under the stated assumptions?
What does it imply if two columns of the exposure-association matrix are proportional across all variants?
How should the coefficient for an upstream exposure be described when the mediator is also included in the model?
Why does a full-rank matrix not, by itself, establish that the instruments are sufficiently strong for stable estimation?