Robust Estimation in Multivariable MR

Why Robust Estimation Is More Demanding With Multiple Exposures

In multivariable Mendelian randomization (MVMR), genetic variants can be associated with more than one exposure. The aim is to estimate the effect of each exposure on an outcome while accounting for the others. This can distinguish, for example, the effect of one correlated risk factor from the effect of another. It also creates a specific vulnerability: a variant may be strongly associated with an exposure marginally but provide little information about that exposure conditional on the other exposures.

Robust estimation addresses departures from a chosen model, such as direct effects of variants on the outcome that are not mediated by the included exposures. It does not remove the need for valid identification assumptions. In particular, neither a robust fit nor a diagnostic statistic can establish that the included genetic variants are independent of confounders or that their direct effects have the structure required by the estimator.

Definition: For $J$ genetic variants and $K$ exposures, let $\hat{\boldsymbol{\gamma}}_j=(\hat{\gamma}_{j1},\ldots,\hat{\gamma}_{jK})^\mathsf{T}$ denote variant $j$'s estimated associations with the exposures, and let $\hat{\Gamma}_{Yj}$ denote its estimated association with the outcome. A summary-association model is $$ \hat{\Gamma}_{Yj}=\hat{\boldsymbol{\gamma}}_j^\mathsf{T}\boldsymbol{\theta}+\alpha_j+\varepsilon_j, $$ where $\boldsymbol{\theta}$ contains the exposure effects, $\alpha_j$ is the variant's direct effect on the outcome not through the included exposures, and $\varepsilon_j$ is sampling error. An estimator's pleiotropy assumption specifies what patterns of $\alpha_j$ it can tolerate.

This model presumes that the included exposures are defined on scales for which the coefficients have an interpretable meaning, that variant associations are appropriately harmonized, and that any dependence between variants is handled in the analysis. For the simple weighted formulas below, assume variants are independent, outcome-association standard errors are known or estimated accurately, and uncertainty in the exposure-association estimates is small enough to be ignored. When that last approximation is poor, weak conditional instruments can distort estimates; methods and diagnostics must address that problem explicitly.

Multivariable IVW and MVMR-Egger

Multivariable inverse-variance weighted (IVW) estimation fits the model with no intercept and no direct effects, using the inverse variance of each outcome association as its weight. In matrix notation, let $\hat{\boldsymbol{\Gamma}}_Y$ be the vector of outcome associations, let $\hat{\boldsymbol{\Gamma}}$ be the $J$-by-$K$ matrix whose rows are $\hat{\boldsymbol{\gamma}}_j^\mathsf{T}$, and let $W$ be the diagonal matrix with entries $1/\operatorname{SE}(\hat{\Gamma}_{Yj})^2$. If the weighted design has full column rank, then

\hat{\boldsymbol{\theta}}_{\mathrm{IVW}} = (\hat{\boldsymbol{\Gamma}}^\mathsf{T}W\hat{\boldsymbol{\Gamma}})^{-1} \hat{\boldsymbol{\Gamma}}^\mathsf{T}W\hat{\boldsymbol{\Gamma}}_Y.

The strongest version of the IVW assumption is that all included variants satisfy the exclusion restriction, so $\alpha_j=0$. Under a weaker balanced-pleiotropy argument, the direct effects must average to zero in the relevant weighted relationship with the exposure associations; importantly, they must not create a systematic association between instrument strength and direct effects that biases the slopes. Merely having positive and negative direct effects is not enough to ensure that condition. IVW also requires adequate conditional instrument strength: the exposure columns must distinguish the effects being estimated.

MVMR-Egger adds an intercept, fitting a weighted regression of outcome associations on all exposure associations. The intercept allows a nonzero average direct effect under its assumptions; it does not give each variant an unrestricted direct-effect parameter. The multivariable InSIDE assumption requires the direct effects to be independent of the variants’ instrument strengths for the exposures, considered jointly in the regression. The design matrix containing the intercept and all exposure columns must have full column rank. As in univariable MR-Egger, uncertainty in estimated variant–exposure associations can attenuate slopes unless the no-measurement-error approximation is adequate or corrected.

Worked Example: A Two-Exposure IVW Fit

Suppose four independent variants have the following hypothetical associations, with equal outcome-association standard errors. Equal weights can be set to 1 without changing the fitted slopes.

Variant	$\hat{\gamma}_{j1}$	$\hat{\gamma}_{j2}$	$\hat{\Gamma}_{Yj}$
1	1	0	0.50
2	0	1	0.10
3	1	1	0.80
4	-1	1	-0.40

The exposure columns have full rank. With unit weights, the cross-product and outcome cross-product are

\hat{\boldsymbol{\Gamma}}^\mathsf{T}\hat{\boldsymbol{\Gamma}} = \begin{pmatrix}3&0\\0&3\end{pmatrix}, \qquad \hat{\boldsymbol{\Gamma}}^\mathsf{T}\hat{\boldsymbol{\Gamma}}_Y = \begin{pmatrix}1.70\\0.50\end{pmatrix}.

Thus

\hat{\boldsymbol{\theta}}_{\mathrm{IVW}} = \begin{pmatrix}3&0\\0&3\end{pmatrix}^{-1} \begin{pmatrix}1.70\\0.50\end{pmatrix} = \begin{pmatrix}1.70/3\\0.50/3\end{pmatrix} = \begin{pmatrix}0.5667\\0.1667\end{pmatrix}.

The fitted outcome associations are $0.5667, 0.1667, 0.7333,$ and $-0.4000$. Subtracting these from the observed values gives residuals $-0.0667,-0.0667,0.0667,$ and $0$, respectively. The slopes are the IVW estimates under the stated model; they are not automatically causal effects if the pleiotropy or instrument assumptions fail.

Worked Example: An Identifiable MVMR-Egger Fit

Consider five hypothetical independent variants with equal outcome-association standard errors. Their exposure associations are $(1,0),(0,1),(1,1),(-1,1),(1,-1)$, and their outcome associations are $0.50,0.30,0.70,-0.10,0.30$, respectively. These associations satisfy

\hat{\Gamma}_{Yj}=0.10+0.40\hat{\gamma}_{j1}+0.20\hat{\gamma}_{j2}.

For example, the third variant has fitted association $0.10+0.40(1)+0.20(1)=0.70$. The intercept-plus-exposure design has full rank: the five rows include $(1,0,0),(1,1,0),(1,0,1)$, which are linearly independent. Therefore its three coefficients are identifiable. The fitted intercept is $0.10$, and the two slopes are $0.40$ and $0.20$; all five residuals are zero in this deliberately exact illustration.

The intercept is an estimate of average direct effect in the fitted model, in outcome-association units; it is not a test that every variant is valid. The slopes have a causal interpretation only if the MVMR-Egger assumptions, including multivariable InSIDE and adequate conditional strength, are credible. With five variants and three fitted coefficients, the residual degrees of freedom are $5-3=2$. This illustrates why an apparently neat fit with few variants may provide little diagnostic information.

Robust Multivariable Estimators

“Robust MVMR” is not one estimator with one universal assumption. Methods protect against different problems, and their identifying assumptions must be stated method by method. Multivariable MR-RAPS extends the robust adjusted profile-score approach: in its usual formulation, it models direct effects as balanced, mean-zero random variation independent of instrument strengths and uses a robust loss to limit the influence of unusually large residuals. Its profile-score framework can also address weak-instrument bias under its model. MR-Lasso-type approaches instead penalize variant-specific direct-effect terms; their conclusions depend on the penalty and on assumptions that enough variants are valid or that invalid direct effects have a suitable sparse structure. Other robust regression choices may downweight large residuals, but downweighting alone does not imply that the remaining variants satisfy a valid-instrument assumption.

These methods are useful when an IVW fit is sensitive to outlying variants or when weak instruments or structured direct effects are concerns. They are not interchangeable. A procedure designed for a sparse set of invalid variants should not be interpreted as if it assumed random balanced pleiotropy, and a robust loss does not protect against every invalid pattern. Report the method, its assumptions, the conditional strength of the exposures, and sensitivity of the estimates to reasonable alternatives.

Conditional Strength and Heterogeneity Diagnostics

For multiple exposures, marginal F statistics do not answer whether each exposure is well predicted after accounting for the others. Sanderson–Windmeijer conditional F statistics assess the strength of the genetic instruments for each exposure conditional on the remaining exposures. They should be considered exposure by exposure. There is no universal threshold that guarantees unbiased MVMR estimates: interpretation depends on the design, number of variants, estimator, and desired precision. A low conditional F statistic warns that the corresponding exposure effect may be weakly identified, even if its marginal associations look strong.

A heterogeneity statistic assesses whether the observed outcome associations vary more around the fitted model than expected from their sampling errors. For independent variants and the no-intercept IVW fit, Cochran’s Q is

Q=\sum_{j=1}^{J} \frac{(\hat{\Gamma}_{Yj}-\hat{\boldsymbol{\gamma}}_j^\mathsf{T} \hat{\boldsymbol{\theta}}_{\mathrm{IVW}})^2} {\operatorname{SE}(\hat{\Gamma}_{Yj})^2}.

When the fitted model is correct and the stated sampling assumptions hold, Q is approximately chi-squared with $J-K$ degrees of freedom. For an intercept-and-slopes model, the corresponding residual degrees of freedom are $J-(K+1)$ when all $K+1$ parameters are estimable. A large Q signals residual heterogeneity, not its cause. It may reflect horizontal pleiotropy, model misspecification, outlying variants, or other problems; a small Q does not prove that the exclusion restriction holds.

Worked Example: Reading Conditional Strength and Q Together

Suppose a hypothetical study includes $J=30$ independent variants for two exposures. The first-stage analysis reports Sanderson–Windmeijer conditional F statistics of 18.0 for exposure 1 and 7.20 for exposure 2. The IVW analysis reports $Q=46.2$. With two fitted slopes, the reference degrees of freedom are $30-2=28$. The 5% upper critical value for a chi-squared distribution with 28 degrees of freedom is approximately 41.34, so $46.2>41.34$; Q indicates more heterogeneity than expected under the fitted IVW model at that reference level.

The conditional F statistic of 7.20 flags exposure 2 as a particular strength concern; it should not be concealed by the stronger value of 18.0 for exposure 1. Meanwhile, the elevated Q indicates model-residual inconsistency, but does not identify which variants or which mechanism produced it. These diagnostics therefore raise two distinct concerns: conditional weakness and heterogeneity. They do not show that a robust estimator will recover the causal effects, nor do they determine which robust method is appropriate.

Comparing Estimates Without Overclaiming

Suppose the same hypothetical data produce IVW estimates of $0.52$ and $0.10$, while a multivariable MR-RAPS analysis reports $0.38$ (SE $0.16$) and $0.24$ (SE $0.20$). For the robust estimates, normal-approximation 95% intervals are calculated as estimate $\pm 1.96\,SE$:

0.38\pm1.96(0.16)=(0.0664,\ 0.6936), \qquad 0.24\pm1.96(0.20)=(-0.152,\ 0.632).

If the first exposure is measured in one unit and the outcome in another, its coefficient and interval are in outcome units per one-unit increase in exposure 1, conditional on exposure 2. The second interval includes zero, so these data do not provide clear evidence of a nonzero effect of exposure 2 at the conventional two-sided 5% level under this approximation. The difference between IVW and MR-RAPS is a sensitivity signal, not proof that the latter is correct. The estimates rely on different assumptions, and uncertainty, conditional weakness, and heterogeneity all matter to their interpretation.

Common Mistakes and Practical Interpretation

Using marginal strength as a substitute for conditional strength. A variant set can predict each exposure separately yet fail to distinguish their effects jointly. Examine conditional F statistics for every exposure.
Calling any robust fit valid. Robust losses or penalties reduce sensitivity to some patterns of outliers or pleiotropy; they do not make invalid instruments valid. State the method-specific assumption.
Treating a large Q as a pleiotropy diagnosis. Q detects residual variation relative to a fitted model. It does not identify the cause, and a small Q is not proof of no pleiotropy.
Fitting an Egger intercept without checking rank. The intercept and all exposure slopes must be separately identifiable. If the intercept column and exposure columns are linearly dependent, those coefficients cannot be uniquely estimated.
Interpreting agreement as validation. Similar estimates across methods can occur even when they share a violated assumption; disagreement is a reason to examine assumptions and data, not to select the preferred answer by convenience.

Key takeaway: Multivariable IVW, MVMR-Egger, and robust multivariable methods answer related questions under different pleiotropy models. Interpret each estimate alongside its assumptions, the conditional strength of every exposure, and heterogeneity diagnostics; none of these checks alone establishes causal validity.

Check Your Understanding

Use the assumptions and calculations in this tutorial to answer the following questions.

Why can marginal instrument-strength statistics be reassuring even when one MVMR exposure is weakly identified?
What does multivariable InSIDE require, and how does the MVMR-Egger intercept differ from allowing every variant its own direct-effect parameter?
For a full-rank MVMR-Egger model with two exposures and 12 variants, how many residual degrees of freedom are available?
What does a large IVW Q statistic indicate, and what does it not establish?
Why is it incorrect to conclude that a robust multivariable estimator is valid solely because its estimate differs from IVW?