Sensitivity analysis for causal mediation analysis with Mendelian randomization

1.   Introduction

In causal inference, Mendelian randomization (MR) is a special instrumental variable (IV) method used to control estimate biases due to unmeasured confounders, associative selection bias (i.e., Berkson’s bias), attenuation by errors (i.e., regression dilution bias), and reverse causation^[1]. One or more genetic variants, such as single nucleotide polymorphism (SNP), serve as IV in MR. MR method was first proposed by Katan^[2], which has been extended substantially, such as two-sample MR^[3] and two-step MR^[4]. The causal directed acyclic graph (DAG) displayed in Fig. 1a represents a particular one-sample MR model, demonstrating the relationships between IV, exposure, outcome, and confounders. Ignoring unobserved confounders can substantially distort the causal effect of the exposure on the outcome using the conventional association analysis methods, and a valid IV helps to correctly identify the causal effect^[5].

Figure  1.  (a) The DAG of a simple MR model with a possibly invalid IV. C1 is the relevance assumption; C2 is the exchangeability assumption; C3 is the exclusion restriction. (b) The DAG of a mediation model with two IVs $Z_M$ and $Z_T$ and four sensitivity parameters $\beta_M^{Z_T}$, $\beta_Y^{Z_T}$, $\beta_T^{Z_M}$, and $\beta_Y^{Z_M}$. Solid lines represent causal effects, and dotted lines represent effects due to deviation of the IV assumptions.

DownLoad: Full-Size Img PowerPoint

In mediation analyses, MR has been widely used to identify which variables serve as mediators in the relationship between an exposure and an outcome, allowing the intervention of these mediators to alter the effects of the exposure. An early MR method to identify mediators was two-step epigenetic Mendelian randomization, which estimated the causal effect of an exposure on an outcome through DNA methylation^[4]. More recently, network Mendelian randomization^[6] and multivariable Mendelian randomization^[7] have been successfully applied to estimate causal effects in MR-based mediation analyses. MR-based mediation analyses have been successfully applied in a variety of topics, such as the extent to which blood pressure, cholesterol, and glucose mediate the excess risk of body mass index (BMI) on stroke^[8], the role of BMI, systolic blood pressure, and smoking behavior in mediating the effect of education on the risk of cardiovascular outcomes^[9], and whether elevated blood pressure estimates the effect of urate on cardiovascular disease risk^[10].

Certain strong assumptions are required for MR to consistently estimate causal effects, especially in mediation analyses. For example, to correctly identify the average causal effect, a valid IV should satisfy the following three assumptions (Fig. 1a): (C1) the IV and the exposure are correlated (relevance assumption); (C2) the IV is independent of all confounders of the exposure and the outcome (exchangeability or independence assumption); and (C3) the IV affects the outcome only through the exposure (exclusion restriction criteria). Refer to Angrist et al.^[5] for more discussions on IVs. The use of IVs in mediation analyses relax the assumption on the unmeasured confounders provided that some stronger assumptions are satisfied for the IVs in addition to the assumptions of linearity and no interaction on the linear model^[6]. However, the IV assumptions are typically hard to satisfy in real studies using the MR method. For example, condition (C3) can be violated due to the horizontal pleiotropy of SNPs^[11], while (C2) and (C3) cannot be fully validated^[12]. Nevertheless, we can still study the impact of the violation of IV assumptions on the MR-based causal inference results. This is the so called sensitivity analysis^[5].

To our knowledge, although sensitivity analyses for IV assumptions in MR-based average causal effect analyses have been widely investigated^[13–15], no sensitivity analysis for IV assumptions has been conducted in the literature for MR-based mediation analyses. This paper aims to fill this gap. Our contributions are three-fold. First, we propose to use only two sensitivity parameters to quantify the bias due to the violation of IV assumptions, making the sensitivity analysis very feasible. Second, we derive consistent causal effect estimators for both individual data and summary data, together with their asymptotic distributions given the two sensitivity parameters. Third, in terms of estimation direction, we demonstrate that the causal effect estimators are robust to some extent in numerical studies.

The rest of this paper is organized as follows. In Section 2, new sensitivity parameters are proposed for mediation analysis, consistent estimators are derived for both direct and indirect effects, and their asymptotic properties are established. In Section 3, numerical studies, including simulation studies, sensitivity analysis, and a real data application, are conducted to demonstrate the desired properties of the proposed estimators. Some concluding remarks are presented in Section 4.

2.   Sensitivity analyses for mediation model

First, we extend the MR model to the MR-based mediation model with possibly invalid IVs (Section 2.1). Then, we propose two sensitivity parameters to quantify the impact of deviation from the IV assumptions on indirect effect and derive consistent indirect causal effect estimators (Section 2.2). Counterparts for the direct effect are provided in Section 2.3.

2.1   Notations, assumptions, and model

$Z_T$ denotes an IV for the exposure $T$, and $Z_M$ denotes an IV for the mediator $M$. Let $Y$ denote the outcome of interest. We consider the causal mediation model with possibly invalid IVs $Z_T$ and $Z_M$, where the exposure $T$, mediator $M$, and outcome $Y$ follow a linear model without interactions. Carter et al.^[16] graphically stated some assumptions required for the validity of $Z_T$ and $Z_M$. Here, we restate these assumptions as follows. First, the IV $Z_T$ satisfies the following three assumptions:

(C1$^{\prime}$) $Z_T$ has a causal effect on the exposure $T$ (relevance assumption).

(C2$^{\prime}$) $Z_T$ is independent of all confounders of the exposure $T$, the mediator $M$, and the outcome $Y$ (independence assumption).

(C3$^{\prime}$) $Z_T$ affects the outcome $Y$ and the mediator $M$ only through the exposure $T$ (exclusion restriction criteria).

Second, the IV $Z_M$ satisfies the following three assumptions:

(C1$^{\prime\prime}$) $Z_M$ has a causal effect on the mediator $M$ (relevance assumption).

(C2$^{\prime\prime}$) $Z_M$ is independent of all confounders of the exposure $T$, the mediator $M$, and the outcome $Y$ (independence assumption).

(C3$^{\prime\prime}$) $Z_M$ affects the outcome $Y$ only through the mediator $M$ (exclusion restriction criteria).

Assumptions (C1$^{\prime}$) and (C3$^{\prime}$) for $Z_T$ may be violated by the association between $Z_T$ and $M$ or $Z_T$ and $Y$ conditional on $T$. Wang et al.^[15] proposed a single sensitivity parameter $\delta$ to measure the total violation of the IV assumption (C2) and (C3). We extend $\delta$ to two sensitivity parameters $\beta_M^{Z_T}$ and $\beta_Y^{Z_T}$ for IV $Z_T$. The parameters $\beta_M^{Z_T}$ and $\beta_Y^{Z_T}$ quantify the total violation of assumptions (C2$^{\prime}$) and (C3$^{\prime}$), one corresponding to the association between $Z_T$ and $M$ and the other corresponding to the association between $Z_T$ and $Y$. Similarly, the assumptions (C2$^{\prime\prime}$) and (C3$^{\prime\prime}$) for $Z_M$ may be violated due to the association between $Z_M$ and $T$ or $Z_M$ and $Y$ conditional on $M$. To quantify the violation of the $Z_M$ assumptions, we utilize two sensitivity parameters $\beta_T^{Z_M},\beta_Y^{Z_M}$, one corresponding to the exposure $T$ and the other to the outcome $Y$. Fig. 1b shows the mediation model with two IVs and four sensitivity parameters $\beta_M^{Z_T}$,$\beta_Y^{Z_T}$,$\beta_T^{Z_M}$, and $\beta_Y^{Z_M}$.

For individual $i$, let $\epsilon_{Yi}$, $\epsilon_{Mi},$ and $\epsilon_{Ti}$ represent the error term of the outcome, the mediator, and the exposure. Suppose $\epsilon_{Yi}$, $\epsilon_{Mi}$, and $\epsilon_{Ti}$, $i = 1,...,n$, are i.i.d. with expectations 0 and are independent of IVs $Z_{Ti}$ and $Z_{Mi}$. Denote ${\boldsymbol{T}} = (T_1,...,T_n)^{\prime}$, ${\boldsymbol{M}} = (M_1,...,M_n)^{\prime}$, ${\boldsymbol{Y}} = (Y_1,...,Y_n)^{\prime}$, ${\boldsymbol{Z}}_T = (Z_{T1},...,Z_{Tn})^{\prime}$, ${\boldsymbol{Z}}_M = (Z_{M1},...,Z_{Mn})^{\prime}$, ${\boldsymbol{\epsilon}}_Y = (\epsilon_{Y1} ,...,\epsilon_{Yn})^{\prime}$, ${\boldsymbol{\epsilon}}_M = (\epsilon_{M1} ,...,\epsilon_{Mn})^{\prime}$, ${\boldsymbol{\epsilon}}_T = (\epsilon_{T1} ,...,\epsilon_{Tn})^{\prime}$. Then, the MR-based mediation model can be written as

where the effect parameters $\beta_Y^{M}$, $\beta_Y^{T}$, $\beta_M^{T}$, $\beta_M^{Z_M}$, and $\beta_T^{Z_T}$ represent the causal effect of the mediator on the outcome, the direct effect of the exposure on the outcome, the causal effect of the exposure on the mediator, the strength of IV $Z_M$, and the strength of IV $Z_T$. We allow unmeasured confounder effects in the exposure-mediator, exposure-outcome, and mediator-outcome relationships, which are reflected by nonzero covariances of the error terms.

2.2   Indirect effect

In mediation analyses, the indirect effect $\beta_Y^M\beta_M^T$ is of great interest, which can be estimated by the product-of-coefficients method^[17,18]. In MR-based mediation analysis, the indirect effect can be consistently estimated under some IV assumptions^[6,7]. In the presence of three or more sensitivity parameters, the sensitivity analysis would be rather complicated. To effectively simplify the sensitivity analysis, we propose two new sensitivity parameters for the indirect effect.

First, we define the following relative parameters:

Then, we define the following two new sensitivity parameters for the indirect effect:

Parameters $\gamma_M$ and $\gamma_Y$ quantify the weighted bias due to the violation of IV assumptions, corresponding to effect parameters $\beta_M^T$ and $\beta_Y^M$, respectively. We derive estimators $\hat{\beta}_{M;\gamma_M}^{T}$ and $\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ for $\beta_M^T$ and $\beta_Y^M$, respectively, as follows.

Proposition 1. With given $\gamma_M$ and $\gamma_Y$, we extend the two-stage least squares (2SLS) estimator^[19] to derive new estimators of $\beta^T_M$ and $\beta^M_Y$:

and

where

Here for vectors ${\boldsymbol{A}}$, ${\boldsymbol{B}}$, and ${\boldsymbol{Z}}_C$, $\hat{\alpha}_{ {\boldsymbol{A}}\vert {\boldsymbol{Z}}_C}^{ {\boldsymbol{B}}} = \dfrac{ {\boldsymbol{B}}^{\prime} {\boldsymbol{M}}_{ {\boldsymbol{Z}}_C} {\boldsymbol{A}}}{ {\boldsymbol{B}}^{\prime} {\boldsymbol{M}}_{ {\boldsymbol{Z}}_C} {\boldsymbol{B}}}$ denotes the regression coefficient of ${\boldsymbol{B}}$ on ${\boldsymbol{A}}$ conditional on ${\boldsymbol{Z}}_C$,

and ${\boldsymbol{A}}\in \{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}, {\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\},C\in\{T,M\}.$

Refer to Appendix A (see Supporting information) for the derivation details of the two estimators given in Proposition 1. The following theorem presents asymptotic properties of the indirect effect estimator $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$. (Refer to Appendix B in Supporting information for proof.)

Theorem 1. Assume that the values of $\gamma_M$ and $\gamma_Y$ are correctly specified. Under some regularity conditions given in Appendix B (see Supporting information), the estimators $\hat{\beta}_{M;\gamma_M}^{T}$, $\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$, and $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ are consistent for $\beta_M^T$, $\beta_Y^M$, and $\beta_M^T\beta_Y^M$, respectively. Furthermore, if the values of $\gamma_M$ and $\gamma_Y$ are specified as $\gamma_M+a$ and $\gamma_Y+b$, then as $n \to \infty$,

where $\sigma_1^2 = E(Z_{T1}^2)$, $\sigma_2^2 = E(Z_{M1}^2), \sigma_{12} = E(Z_{T1}Z_{M1})$, and

The expressions of $\omega_{AB}$ and $\varOmega_{AB}$, $A,B\in\{T,M,Y\}$, are given in Appendix B (see Supporting information).

Since the asymptotic distributions given in Theorem 1 depend on some unknown parameters (i.e., $\beta^{Z_T}_Y, \beta_Y^{Z_M}, \sigma_{Y Y}, \sigma_{M M}, \sigma_{T T}$) not estimable using observed data, we cannot directly apply Theorem 1 to make inferences on the indirect effect. Instead, we can make inferences on the indirect effect given sensitivity parameters by applying the following bootstrap algorithm (Algorithm 1), provided individual level data are available.

Note that the indirect effect estimator $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ can only be calculated using individual data, and it cannot be calculated using only summary data reported in the literature. Fortunately, the estimators given in Proposition 1 can be simplified so that they depend only on regression coefficient estimates obtained from summary data, as detailed below (refer to Appendix C in the Supporting information for a derivation).

Algorithm 1 A bootstrap algorithm for inferring indirect effects
Require: original sample ${\boldsymbol{X}}_1, {\boldsymbol{X}}_2,..., {\boldsymbol{X}}_n$
Ensure: t-test statistic and confidence interval of the indirect effect $\beta_M^T\beta_Y^M$
1: compute $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ using ${\boldsymbol{X}}_1, {\boldsymbol{X}}_2,..., {\boldsymbol{X}}_n$
2: for $i=1$ to $B$ do
3: 　generate sample ${\boldsymbol{Y}}_1^{i}, {\boldsymbol{Y}}_2^{i},..., {\boldsymbol{Y}}_n^{i}$ by sampling with replacement from ${\boldsymbol{X}}_1, {\boldsymbol{X}}_2,..., {\boldsymbol{X}}_n$
4: 　compute indirect effect estimate $\hat{\beta}_{M;\gamma_M}^{T,(i)}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,(i)}$ using the sample ${\boldsymbol{Y}}_1^{i}, {\boldsymbol{Y}}_2^{i},..., {\boldsymbol{Y}}_n^{i}$
5: end for
6: compute the standard error $\hat{\sigma}_{\beta}$ of the estimates $\hat{\beta}_{M;\gamma_M}^{T,(i)}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,(i)},i=1,2,...,B$;
7: return Z-test statistic $B^{1/2}\hat\beta/\hat\sigma_\beta$ and confidence interval
$\left(\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}-B^{-\frac{1}{2}}Z_{\alpha/2}/\hat{\sigma}_{\beta},\ \hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}+B^{-\frac{1}{2}}Z_{\alpha/2}/\hat{\sigma}_{\beta}\right).$

 | Show Table

DownLoad: CSV

Corollary 1. If we additionally assume that IVs $Z_T$ and $Z_M$ are independent, then $\hat{\beta}_{M;\gamma_M}^{T}$ and $\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ given in Proposition 1 can be simplified as

and

respectively, where $\hat{\beta}_{M;\gamma_M}^{Z_M,\perp} = \hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_M}-\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_M}\cfrac{\hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_T}}{\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_T}}+\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_M}\gamma_M$, and $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ is the regression coefficient of ${\boldsymbol{B}}$ on ${\boldsymbol{A}}$ for ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$ and ${\boldsymbol{B}}\in \{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$. Furthermore, if the true values of $\gamma_M$ and $\gamma_Y$ are specified to be $\gamma_M+a$ and $\gamma_Y+b$, respectively, then

in probability as $n \to \infty$.

In practice, we can use $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}$ instead of $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$ if only summary data are available. We establish some asymptotic properties of $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}$ in the following corollary:

Corollary 2. Assume that 1) the values of $\gamma_M$ and $\gamma_Y$ are correctly specified, 2) IVs $Z_T$ and $Z_M$ are independent, and 3) the regression coefficient estimators $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$, ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$, ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$, are independent. Then, $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}$ is consistent for $\beta_M^T\beta_Y^M$, and its variance can be approximated by

where $\sigma_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ is the standard error of $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ for ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$ and ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$. The detailed expressions of $\frac{\partial}{\partial \hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}}(\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp})$ are given in Appendix D (see Supporting information).

Using Corollary 2, we can easily construct a $(1-\alpha)$-confidence interval of the indirect effect as follows:

where $\hat\sigma_{\gamma_M,\gamma_Y}^{2,\perp}$ is a consistent estimator of $\sigma_{\gamma_M,\gamma_Y}^{2,\perp}$ using a plug-in rule. Furthermore, the Z-test statistic for the hypothesis test problem $H_0:\beta_M^T\beta_Y^M = 0\leftrightarrow H_a:\beta_M^T\beta_Y^M\not = 0$ takes the form $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}/\sqrt{\hat\sigma_{\gamma_M,\gamma_Y}^{2,\perp}}$, whose null limit distribution is standard normal.

2.3   Direct effect

In Section 2.2, we derive consistent estimators of the indirect effect given two sensitivity parameters and its asymptotic normality. We consider estimating the direct effect in this subsection. The direct effect can be estimated in two distinct ways. One is to employ an estimator by the fact that the direct effect is defined as the difference between the total effect (denoted by $\beta$) of exposure on the outcome and the indirect effect $\beta_M^T\beta_Y^M$ (difference method). The other is to extend the multivariable Mendelian randomization (MVMR)^[7], which utilizes a 2SLS estimator to estimate the direct effect (MVMR method).

The difference method first estimates $\beta_M^T\beta_Y^M$, refer to Proposition 1. The 2SLS estimator of the total effect $\beta$ can be consistent given four sensitivity parameters $\beta_M^{Z_T}$, $\beta_Y^{Z_T}$, $\beta_T^{Z_M}$, and $\beta_Y^{Z_M}$ under model (1). Here we propose to use sensitivity parameters $\gamma_M$ and $\rho_Y^{Z_T}$ to obtain the total effect. Specifically, define ${\boldsymbol{M}}_{ {\boldsymbol{Z}}_M} = {\boldsymbol{I}}_n- {\boldsymbol{Z}}_M( {\boldsymbol{Z}}_M^{\prime} {\boldsymbol{Z}}_M^{\prime})^{-1} {\boldsymbol{Z}}_M^{\prime}$, then the 2SLS estimator of the total effect can be represented as (refer to Appendix E in the Supporting information for a derivation)

where ${\boldsymbol{M}}_{ {\boldsymbol{Z}}_M}$ is defined in (3). Note that $\widehat{\beta}_{\mathrm{2SLS}}\to \beta+\beta_Y^M\gamma_M+\rho_Y^{Z_T}$ in probability as $n\to\infty$, we propose to use the direct effect estimator

which is consistent for $\beta_Y^T$ given three sensitivity parameters $\gamma_Y$, $\gamma_M$, and $\rho_Y^{Z_T}$.

No sensitivity analysis for IV assumptions has been conducted in the literature for the MVMR method, and we extend this method to model (1) with possibly invalid IVs. In what follows, we show that the MVMR method requires only two sensitivity parameters $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$ to estimate the direct effect, compared to the difference method. Denote

where ${\boldsymbol{\alpha}}$ is a sensitivity parameter vector. Model (1) can be rewritten as

Under the above model, the MVMR estimator can be written as

where

which requires specifying two sensitivity parameters $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$ in estimating the direct effect $\beta_Y^{T}$. Then, we have the following asymptotic results:

Theorem 2. With correctly specified $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$, under some regularity conditions presented in Appendix F (see Supporting information), the MVMR estimator $\hat{\beta}_{X;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{\mathrm{MVMR}}$ is consistent for ${\boldsymbol{\beta}}_{X}$. Furthermore, if the true values of $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$ are specified to be $\beta_Y^{Z_T}+c$ and $\beta_Y^{Z_M}+d$, respectively, then as $n \to \infty$,

where $\alpha_M^{Z_M} = \beta_M^{Z_M}+\beta_M^T\beta_T^{Z_M}$, $\alpha_M^{Z_T} = \beta_M^{Z_T}+\beta_M^T\beta_T^{Z_T}$, $\sigma_1^2 = E(Z_{T1}^2)$, $\sigma_2^2 = E(Z_{M1}^2)$, and $\sigma_{12} = E(Z_{T1}Z_{M1})$. The expressions of $\omega_{TT}$, $\omega_{TM}$, and $\omega_{MM}$ are given in Appendix F (see Supporting information).

Since the asymptotic distributions given in Theorem 2 depend on some unknown parameters (i.e., $\beta_M^{Z_T}, \sigma_{Y Y}, \sigma_{M M}, \sigma_{T T}$) that are not estimable using observed data, we cannot directly apply Theorem 2 to make inferences on the direct effect given sensitivity parameters. Instead, we can apply a bootstrap algorithm similar to the one given in Section 2.2. Furthermore, the MVMR estimators of the direct effect $\beta_Y^{T}$ can be simplified so that it depends only on regression coefficient estimates obtained from summary data, as detailed below (refer to Appendix G in the Supporting information for a derivation).

Corollary 3. If we additionally assume that IVs $Z_T$ and $Z_M$ are independent, then $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR}}$ can be simplified as

where $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ is the regression coefficient of ${\boldsymbol{B}}$ on ${\boldsymbol{A}}$ for ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$ and ${\boldsymbol{B}}\in \{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$. Furthermore, if the values of $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$ are specified as $\beta_Y^{Z_T}+c$ and $\beta_Y^{Z_M}+d$, respectively, then,

in probability as $n \to \infty$.

In practice, we can use $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp}$ instead of $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR}}$ if only summary data are available. We establish asymptotic properties of $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp}$ as follows.

Corollary 4. Assume that 1) the values of $\beta_Y^{Z_T}$ and $\beta_Y^{Z_M}$ are correctly specified, 2) IVs $Z_T$ and $Z_M$ are independent, and 3) the regression coefficient estimators $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$, ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$, ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$, are independent. Then, $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp}$ is consistent for $\beta_Y^T$, and its variance can be approximated by

where $\sigma_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ is the standard error of $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$ for ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$ and ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$, and the detailed expressions of $\frac{\partial}{\partial \hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}}(\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp})$ are given in Appendix H (see Supporting information).

Using Corollary 4, we can easily construct a $(1-\alpha)$-confidence interval of the direct effect as follows:

where $\hat\sigma_{\gamma_M,\gamma_Y}^{2,\mathrm{MVMR},\perp}$ is a consistent estimator of $\sigma_{\gamma_M,\gamma_Y}^{2,\mathrm{MVMR},\perp}$ by applying the plug-in rule. Furthermore, the Z-test statistic for hypothesis problem $H_0$ : $\beta_Y^T = 0\leftrightarrow H_a:\beta_Y^T\not = 0$ takes the form $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp}/\sqrt{\hat\sigma_{\gamma_M,\gamma_Y}^{2,\mathrm{MVMR},\perp}}$, whose null limit distribution is standard normal.

3.   Numerical studies

3.1   Simulation studies with correctly specified sensitivity parameters

We conduct four simulation studies to demonstrate our method’s effectiveness in estimating the indirect effect and direct effect. The two-step MR method^[4] and the MVMR method are also included for the purpose of comparison. The first simulation study focuses on the performance of our method and the two-step MR method using individual data to estimatie the indirect effect. The second simulation study compares the MVMR method with the two-step MR method for estimating the direct effect. The third simulation study focuses on the performance of our method using summary data for estimating the indirect effect. The fourth simulation study focuses on the performance of our method using summary data for estimating the direct effect.

Data are generated according to the following model:

where $i = 1,\ldots,n$, $Z_T$ and $Z_M$ are independent biallelic genetic variants with a common minor allele frequency 0.3, which are both coded according to an addiditive mode of inheritance. In this linear model, the correlation coefficients of the residuals quantify the effects due to the confounders between the exposure, mediator, and outcome, as specified by

With a sample size of $n =$1000, 2000, or 4000, we set $\beta_T^{Z_T} = 0.4$ and $\beta_M^{Z_M} = 0.5$, which results in an average F statistic greater than 10 (F statistic reflects the ‘strength’ of the IV), indicating that the potential bias from weak IV is acceptable ^[20].

In the first simulation study, we consider three sets of sensitivity parameters $(\gamma_M,\gamma_Y)$ = (0, 0), (0.25, 0.25), and (0.75, 0.75), representing no violation, mild violation, and strong violation, respectively, of IV assumptions. No violation indicates that all the IV assumptions are satisfied.

Table 1 summarizes the estimation results based on 1000 replicates of simulations. Our method provides a smaller bias and standard error (SE) than the two-step MR method. In the mild violation and strong violation scenarios, the two-step estimates are biased, while the estimates of our method are still virtually unbiased. As the two-step MR method does not consider the case where the IV assumptions are violated, it is not surprisingly biased when the sensitivity parameters are not 0. In the no violation scenario, the bias and SE of our method are also lower than those of the two-step MR method. This is because our method incorporates more information than the two-step MR estimator. In all scenarios, the bias of our method decreases with the sample size, as expected. In the scenario in which both $\beta_M^T$ and $\beta_Y^M$ equal 0, the coverage probability (CP) of confidence intervals approaches 1 because the standard error of the product-of-coefficients estimator is generally conservative^[21]. In summary, our method outperforms the two-step MR method.

Table  1.  Summary of simulation results for the indirect effect using individual data.

Sample size	$(\beta_{M}^T,\beta_{Y}^M)$	Our method w/o independence2					Our method w/ independence3					Two-step MR
		Bias	SE	SEE	CP		Bias	SE	SEE	CP		Bias	SE	SEE	CP
No violation ($(\gamma_M,\gamma_Y)=(0,0)$)
1000	(0.5,0.5)	−0.006	0.061	0.070	0.968		−0.004	0.078	0.091	0.967		−0.007	0.083	0.093	0.957
2000		−0.004	0.042	0.046	0.956		−0.002	0.057	0.061	0.945		−0.005	0.055	0.062	0.953
4000		0.000	0.030	0.032	0.954		0.000	0.043	0.043	0.942		−0.002	0.041	0.042	0.948
1000	(0,0)	−0.002	0.009	0.019	1.000		0.000	0.014	0.028	1.000		0.000	0.013	0.026	1.000
2000		−0.001	0.005	0.009	1.000		0.000	0.007	0.013	1.000		0.000	0.006	0.012	1.000
4000		−0.001	0.002	0.004	1.000		0.000	0.003	0.006	1.000		0.000	0.003	0.006	1.000
Mild violation ($(\gamma_M,\gamma_Y)=(0.25,0.25)$)
1000	(0.5,0.5)	−0.007	0.080	0.086	0.964		−0.010	0.091	0.102	0.970		0.157	0.102	0.106	0.732
2000		−0.003	0.054	0.056	0.957		0.000	0.063	0.067	0.964		0.157	0.070	0.073	0.448
4000		−0.003	0.037	0.038	0.950		0.000	0.046	0.046	0.963		0.158	0.053	0.050	0.137
1000	(0,0)	−0.003	0.015	0.025	1.000		−0.001	0.017	0.033	1.000		0.013	0.029	0.034	0.997
2000		−0.001	0.007	0.012	1.000		−0.001	0.009	0.015	1.000		0.013	0.020	0.021	0.983
4000		−0.001	0.003	0.006	1.000		0.000	0.004	0.007	1.000		0.012	0.014	0.014	0.951
Strong violation ($(\gamma_M,\gamma_Y)=(0.75,0.75)$)
1000	(0.5,0.5)	−0.123	0.450	5.030	0.977		−0.070	0.567	2.933	0.978		0.659	0.134	0.134	0.000
2000		−0.045	0.168	0.300	0.927		−0.046	0.173	0.334	0.957		0.660	0.089	0.093	0.000
4000		−0.017	0.102	0.110	0.938		−0.020	0.098	0.116	0.957		0.664	0.066	0.065	0.000
1000	(0,0)	−0.031	0.078	3.006	0.999		−0.049	0.200	7.663	0.999		0.225	0.077	0.085	0.283
2000		−0.014	0.033	0.108	0.999		−0.017	0.064	0.102	0.998		0.223	0.057	0.058	0.035
4000		−0.007	0.016	0.029	1.000		−0.009	0.030	0.051	1.000		0.224	0.032	0.032	0.000
1 Note: $\beta_{M}^T\beta_{Y}^M$, true indirect effect; Bias, mean of estimated indirect effects minus true value; SE, standard error of indirect effect estimates; SEE, mean estimated standard errors; CP, coverage probability of confidence interval at level $95\%$. 2 Independence assumption between IVs $Z_T$ and $Z_M$ IS NOT imposed on the estimator $\hat{\beta}_{M;\gamma_M}^{T}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M}$. 3 Independence assumption between IVs $Z_T$ and $Z_M$ IS imposed on the estimator $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}$.

 | Show Table

DownLoad: CSV

In the second simulation study, we compare the performance of two indirect effect estimators, i.e., the MVMR method and the two-step MR method. We consider three sets of sensitivity parameters: $(\beta_T^{Z_M},\beta_M^{Z_T},\beta_Y^{Z_M},\beta_Y^{Z_T})$ = (0, 0, 0, 0), (0.1, 0.125, 0.2, 0.2), (0.3, 0.375, 0.6, 0.6), representing no violation, mild violation, and strong violation, respectively, of IV assumptions. Simulation results based on 1000 replicates are summarized in Table 2. The bias, SE, and SEE of the MVMR method are lower than those of the two-step MR method in the no violation scenario. In other scenarios, the MVMR method is still nearly unbiased, while the two-step MR method is biased. Furthermore, the SE of the MVMR method increases with sensitivity parameters, as expected. In the strong violation scenario, the SE and SEE of the MVMR method are relatively large. Nevertheless, the SE, SEE, and CP become acceptable when the sample size is 4000.

Table  2.  Summary of simulation results for the direct effect using individual data.

Sample size	$\beta_{Y}^T$	MVMR method					Two-step MR method
		Bias	SE	SEE	CP		Bias	SE	SEE	CP
No violation ($(\beta_T^{Z_M},\beta_M^{Z_T},\beta_Y^{Z_M},\beta_Y^{Z_T})=(0,0,0,0)$)
1000	0.5	0.003	0.099	0.105	0.960		−0.009	0.126	0.135	0.970
2000		0.001	0.073	0.073	0.955		−0.004	0.085	0.090	0.957
4000		−0.005	0.051	0.051	0.949		−0.002	0.062	0.062	0.946
1000	0	0.001	0.100	0.105	0.972		−0.009	0.131	0.133	0.957
2000		0.000	0.073	0.073	0.955		−0.005	0.091	0.090	0.952
4000		−0.001	0.051	0.051	0.955		−0.002	0.063	0.062	0.943
Mild violation ($(\beta_T^{Z_M},\beta_M^{Z_T},\beta_Y^{Z_M},\beta_Y^{Z_T})=(0.1,0.125,0.2,0.2)$)
1000	0.5	−0.005	0.127	0.130	0.979		0.344	0.100	0.105	0.098
2000		0.000	0.085	0.088	0.973		0.343	0.072	0.071	0.005
4000		0.000	0.059	0.062	0.966		0.346	0.048	0.050	0.000
1000	0	0.002	0.119	0.130	0.978		0.380	0.097	0.101	0.032
2000		0.004	0.086	0.089	0.960		0.376	0.069	0.070	0.001
4000		−0.003	0.063	0.062	0.950		0.376	0.047	0.048	0.000
Strong violation ($(\beta_T^{Z_M},\beta_M^{Z_T},\beta_Y^{Z_M},\beta_Y^{Z_T})=(0.3,0.375,0.6,0.6)$)
1000	0.5	0.018	0.519	5.196	1.000		0.080	0.196	0.205	0.947
2000		−0.003	0.277	1.283	0.997		0.095	0.139	0.139	0.885
4000		0.006	0.180	0.203	0.989		0.092	0.097	0.097	0.837
1000	0	0.006	0.477	4.538	0.999		0.381	0.165	0.172	0.401
2000		−0.008	0.265	0.993	0.996		0.380	0.117	0.116	0.110
4000		−0.003	0.188	0.204	0.985		0.372	0.081	0.081	0.004
$\beta_{Y}^T$, true direct effect; Bias, mean estimated direct effects minus true value; SE, standard error of indirect effect estimates; SEE, mean estimated standard errors; CP, coverage probability of confidence intervals at level $95\%$.

 | Show Table

DownLoad: CSV

In the third simulation study, we investigate the performance of our method in various situations for estimating the indirect effect. A strong assumption is required to construct a confidence interval of the indirect effect based on Corollary 2. That is, the six regression coefficient estimators: $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$, ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$, ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$, are independent. We call this a fully-independent assumption. The fully-independent assumption is hard to satisfy in real studies, so we also relax this strong assumption by generating data under the so-called partially-independent assumption: $(\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_T},\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_M})$, $(\hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_T},\hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_M})$, and $(\hat{\alpha}_{ {\boldsymbol{Y}}}^{ {\boldsymbol{Z}}_T},\hat{\alpha}_{ {\boldsymbol{Y}}}^{ {\boldsymbol{Z}}_M})$ are independent. In what follows, we investigate the performance of our method with data generated under these two assumptions separately.

Data are generated according to model (12), with sensitivity parameters being $\beta_T^{Z_M} = 0.1$, $\beta_M^{Z_T} = 0.125$, and $\beta_Y^{Z_M} = \beta_Y^{Z_T} = 0.2$. To obtain summary data satisfying the above assumptions, we split the generated data into six subsets of equal size, which are used to calculate six regression coefficient estimators $\hat{\alpha}_{ {\boldsymbol{A}}}^{ {\boldsymbol{B}}}$, ${\boldsymbol{A}}\in\{ {\boldsymbol{T}}, {\boldsymbol{M}}, {\boldsymbol{Y}}\}$ and ${\boldsymbol{B}}\in\{ {\boldsymbol{Z}}_T, {\boldsymbol{Z}}_M\}$, respectively. These six regression coefficient estimators satisfy the fully-independent assumption. Similarly, three subsets are randomly selected from the above six subsets to obtain regression coefficient estimators that satisfy the partially-independent assumption. The first one is used to calculate $\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_T}$ and $\hat{\alpha}_{ {\boldsymbol{T}}}^{ {\boldsymbol{Z}}_M}$, the second one is used to calculate $\hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_T}$ and $\hat{\alpha}_{ {\boldsymbol{M}}}^{ {\boldsymbol{Z}}_M}$, and the third one is used to calculate $\hat{\alpha}_{ {\boldsymbol{Y}}}^{ {\boldsymbol{Z}}_T}$ and $\hat{\alpha}_{ {\boldsymbol{Y}}}^{ {\boldsymbol{Z}}_M}$.

Simulation results based on 1000 replicates are reported in Table 3. It appears that the estimated indirect effects in the partially-independent situation are comparable to those in the fully-independent situation. When both $\beta_Y^M$ and $\beta_M^T$ are nonzero, the estimates are nearly unbiased, and CPs are close to the nominal level 0.95. When both $\beta_Y^M$ and $\beta_M^T$ equal 0, the estimates are still unbiased, but SEEs are larger than SEs, and the CPs are close to 1. This finding is similar to that of the simulation study with individual level data. In conclusion, the violation of the fully-independent assumption in the inference of the indirect effect using summary data has little impact on the analysis results, so we can apply Corollary 2 to infer the indirect effect by relaxing assumption 2). That is, we can use summary data from three instead of independent studies, as done in our real data application described in Section 3.3.

Table  3.  Summary of simulation results for the indirect effect using summary data.

Sample size	Fully-independent assumption2					Partially-independent assumption3
	Bias	SE	SEE	CP		Bias	SE	SEE	CP
$\beta_Y^M=0.5,\beta_M^T=0.5$
1000	0.015	0.157	0.155	0.925		0.008	0.149	0.152	0.917
2000	0.005	0.106	0.102	0.928		0.009	0.102	0.103	0.934
4000	0.004	0.070	0.071	0.947		0.002	0.068	0.070	0.940
$\beta_Y^M=0.5,\beta_M^T=0$
1000	0.006	0.074	0.073	0.984		0.003	0.074	0.073	0.989
2000	−0.002	0.048	0.048	0.973		0.002	0.048	0.049	0.975
4000	0.002	0.035	0.034	0.959		0.001	0.033	0.034	0.961
$\beta_Y^M=0,\beta_M^T=0.5$
1000	0.000	0.064	0.066	0.999		0.004	0.067	0.066	1.000
2000	0.000	0.043	0.043	0.989		−0.002	0.042	0.042	0.989
4000	0.000	0.028	0.029	0.976		0.000	0.029	0.029	0.978
$\beta_Y^M=0,\beta_M^T=0$
1000	0.001	0.015	0.020	1.000		0.000	0.016	0.020	1.000
2000	0.000	0.007	0.009	1.000		0.000	0.008	0.009	1.000
4000	0.000	0.004	0.005	1.000		0.000	0.004	0.005	1.000
1 Sample size, size of each subset. $\beta_{M}^T\beta_{Y}^M$, true indirect effect. Bias, mean estimated indirect effects minus true value; SE, standard error of indirect effect estimates; SEE, mean of estimated standard errors; CP, coverage probability of confidence intervals at level $95\%$. 2 The six estimators $\hat{\alpha}_{\boldsymbol{A}}^{\boldsymbol{B}}$, ${\boldsymbol{A}}\in\{{\boldsymbol{T}},{\boldsymbol{M}},{\boldsymbol{Y}}\}$, ${\boldsymbol{B}}\in\{{\boldsymbol{Z}}_T,{\boldsymbol{Z}}_M\}$, are independent. 3 The three estimators $(\hat{\alpha}_{{\boldsymbol{T}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{T}}}^{{\boldsymbol{Z}}_M})$, $(\hat{\alpha}_{{\boldsymbol{M}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{M}}}^{{\boldsymbol{Z}}_M})$, and $(\hat{\alpha}_{{\boldsymbol{Y}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{Y}}}^{{\boldsymbol{Z}}_M})$ are independent.

 | Show Table

DownLoad: CSV

In the fourth simulation study, summary data are generated under the fully-independent assumption or the partially-independent assumption. The data generation process is similar to that of the third simulation study, with sensitivity parameters where $\beta_T^{Z_M} = 0.1$, $\beta_M^{Z_T} = 0.125$, and $\beta_Y^{Z_M} = \beta_Y^{Z_T} = 0.2$. Simulation results based on 1000 replicates are reported in Table 4.

Table  4.  Summary of simulation results for the direct effect using summary data.

Sample size	Fully-independent assumption2					Partially-independent assumption3
	Bias	SE	SEE	CP		Bias	SE	SEE	CP
$\beta_Y^T=0.5$
1000	0.002	0.196	0.199	0.969		−0.001	0.190	0.193	0.958
2000	0.004	0.130	0.132	0.960		0.000	0.131	0.134	0.954
4000	−0.001	0.093	0.094	0.959		−0.001	0.092	0.094	0.954
$\beta_Y^T=0$
1000	−0.005	0.133	0.138	0.956		−0.004	0.132	0.137	0.960
2000	0.002	0.096	0.095	0.960		0.001	0.097	0.095	0.956
4000	−0.003	0.067	0.067	0.960		−0.001	0.065	0.066	0.955
$\beta_Y^T=-0.5$
1000	−0.002	0.141	0.137	0.958		−0.006	0.138	0.137	0.954
2000	−0.007	0.098	0.096	0.951		−0.005	0.097	0.095	0.956
4000	−0.002	0.067	0.067	0.947		−0.002	0.065	0.066	0.953
1 Sample size, size of each subset; $\beta_{Y}^T$, true direct effect. Bias, mean estimated indirect effects minus true value; SE, standard error of indirect effect estimates; SEE, mean of estimated standard errors; CP, coverage probability of confidence intervals at level $95\%$. 2 The six estimators $\hat{\alpha}_{\boldsymbol{A}}^{\boldsymbol{B}}$, ${\boldsymbol{A}}\in\{{\boldsymbol{T}},{\boldsymbol{M}},{\boldsymbol{Y}}\}$, ${\boldsymbol{B}}\in\{{\boldsymbol{Z}}_T,{\boldsymbol{Z}}_M\}$, are independent. 3 The three estimators $(\hat{\alpha}_{{\boldsymbol{T}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{T}}}^{{\boldsymbol{Z}}_M})$, $(\hat{\alpha}_{{\boldsymbol{M}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{M}}}^{{\boldsymbol{Z}}_M})$, and $(\hat{\alpha}_{{\boldsymbol{Y}}}^{{\boldsymbol{Z}}_T},\hat{\alpha}_{{\boldsymbol{Y}}}^{{\boldsymbol{Z}}_M})$ are independent.

 | Show Table

DownLoad: CSV

Again, the estimated direct effects in the partially-independent situation are comparable to those in the fully-independent situation. In all scenarios, the estimates are nearly unbiased, SEEs are close to SEs, and CPs are close to the nominal level 0.95. This finding is similar to that of the third simulation study, that is, the violation of the fully-independent assumption using summary data has little impact on the inference results of the direct effect. Therefore, we can apply Corollary 4 to infer the direct effect by relaxing assumption 2).

3.2   Sensitivity analysis with misspecified sensitivity parameters

We conduct a sensitivity analysis to illustrate the robustness of the estimator $\hat{\beta}_{M;\gamma_M+a}^{T,\perp}\hat{\beta}_{Y;\gamma_M+a,\gamma_Y+b}^{M,\perp}$ for the indirect effect and the MVMR estimator $\hat{\beta}_{Y;\beta_Y^{Z_T}+c,\beta_Y^{Z_M}+d}^{T,\mathrm{MVMR},\perp}$ for the direct effect. According to (6) and (11), the asymptotic means of the two estimators are

and

respectively.

For any given $\beta^T_M$ and $\beta^M_Y$, we are interested in determining the range of $a$ and $b$ such that the asymptotic mean given in (13) has the same sign as that of the true indirect effect $\beta^T_M\beta^M_Y$. This issue can be addressed by solving an inequality. In this way, we can assess the robustness of this estimator for the indirect effect. Similarly, we can assess the robustness of the MVMR estimator given in (11) for the direct effect. In what follows, we illustrate the robustness of the two estimators through two examples.

In the first example, we set $\gamma_M = \gamma_Y = \beta_T^{Z_M} = \beta_M^{Z_T} = \beta_Y^{Z_M} = \beta_Y^{Z_T} = 0$, $\beta^{Z_T}_T = 0.4$, $\beta^{Z_M}_M = 0.5$, and $\beta_Y^M = \beta_M^T = 0.5$ so that the corresponding indirect effect is positive (i.e., $\beta_Y^M\beta_M^T = 0.25$). In this example, the asymptotic mean of $\hat{\beta}_{M;\gamma_M+a}^{T,\perp}\hat{\beta}_{Y;\gamma_M+a,\gamma_Y+b}^{M,\perp}$ reduces to $(\beta^T_M-a)(\beta^M_Y+b)$. A sufficient condition for a positive asymptotic mean is $a<\beta^T_M = 0.5$ and $b>-\beta^M_Y = -0.5$.

In the second example, we set $\beta_Y^T = 0.5$ (the true direct effect), $\beta_M^T = \beta_Y^M = 0$, $\beta^{Z_T}_T = 0.4$, $\beta^{Z_M}_M = 0.5$, and $(\beta_T^{Z_M},\beta_M^{Z_T},\beta_Y^{Z_M},\beta_Y^{Z_T})$ = (0.1, 0.1, 0, 0). In this example, the asymptotic mean of $\hat{\beta}_{Y;\beta_Y^{Z_T}+c,\beta_Y^{Z_M}+d}^{T,\mathrm{MVMR},\perp}$ reduces to $\beta_Y^T-\dfrac{50}{19}c+\dfrac{10}{19}d$. As a result, the sufficient and necessary condition for a positive direct effect is $50c-10d>-9.5$.

3.3   A real data application

We apply the proposed estimators to study the causal relationship between BMI and stroke mediated by systolic blood pressure (SBP). In the mediation model we considered, BMI is the exposure, SBP is the mediator, and stroke is the outcome, and we want to estimate the causal indirect and direct effects. SNPs rs2867125 on chromosome 2 and rs2681492 on chromosome 12 are chosen as IVs for BMI and SBP, respectively, since they have been shown to be strongly associated with BMI^[22] and SBP^[23], respectively. Summary level data (i.e., point estimates, standard errors, and p-values for the associated regression parameters) are extracted from three independent GWAS^[24–26] in the IEU OpenGWAS project^[27] (Dataset IDs: ieu-b-4816, ukb-b-20175, and ebi-a-GCST005838), as presented in Table 5. We conduct sensitivity analyses to demonstrate the impact of the violation of IV assumptions on causal inference results.

Table  5.  Summary of the genetic effects of two genetic IVs on BMI, SBP, and stroke.

Variable	rs2867125				rs2681492
	ES	SE	P-val		ES	SE	P-val
BMI	0.297	0.024	2.6E−34		−0.034	0.026	2.0E−1
SBP	0.014	0.003	1.6E−7		−0.038	0.003	7.9E−43
Stroke	0.001	0.011	9.0E−1		−0.036	0.010	2.7E−4
The summary data for the relationship between BMI, SBP, and stroke are taken from three independent GWAS in the IEU OpenGWAS project[27]. ES, genetic effect estimate; SE, standard error of effect estimate; P-val, p-value for significance testing of genetic effect.

 | Show Table

DownLoad: CSV

Since SNPs rs2867125 and rs2681492 locate in different chromosomes, the independence assumption is satisfied for the two SNPs. Furthermore, summary level data are obtained from three independent datasets. Thus, the summary level statistics $\hat{\beta}_{M;\gamma_M}^{T,\perp}\hat{\beta}_{Y;\gamma_M,\gamma_Y}^{M,\perp}$ and $\hat{\beta}_{Y;\beta_Y^{Z_T},\beta_Y^{Z_M}}^{T,\mathrm{MVMR},\perp}$ are applicable to this dataset according to the simulation results under the partially dependent assumption in Section 3.1.

First we report the analysis results under the IV assumptions. With $\gamma_Y = 0$ and $\gamma_M = 0$, the estimates of $\beta_M^T$ and $\beta_Y^M$ are 0.048 and 0.996 (SEs = 0.010 and 0.292), respectively, and the corresponding estimated indirect effect is $0.048 \times 0.996 = 0.048$ (SE = 0.017). As a result, the causal indirect effect of BMI on stroke mediated by SBP is significant at level 5% (p-value = 0.006), coinciding with previous finding^[8]. Similarly, with $\beta_Y^{Z_T} = \beta_Y^{Z_M} = 0$, the estimate of $\beta_Y^T$ is $-0.043$ (SE = 0.043), so that the direct effect is not significant at level 5% (p-value = 0.309).

Next we report sensitivity analysis results. By varying the specified values of sensitivity parameters, we can assess the robustness of the indirect and direct effect estimates. We set the range of specified sensitivity parameters as $-0.125\le \gamma_M,\gamma_Y,\beta_Y^{Z_T},\beta_Y^{Z_M} \le 0.125$ in assessing the robustness of indirect effect and direct effect. The contours of the estimated effects and p-values are depicted in Fig. 2.

Figure  2.  Contour plots of the corresponding significance test p-values (column 1) and mean estimated causal effects (column 2) with various specified sensitivity parameters. Row A, indirect effect; row B, direct effect.

DownLoad: Full-Size Img PowerPoint

As shown in Fig. 2a1 and Fig. 2a2, the indirect effect is significantly positive at level 0.05 (estimated indirect effect $>0.02$) when $\gamma_M < 0.015$. It appears that the result is largely determined by $\gamma_M$ instead of $\gamma_Y$. This is because the Z-test statistic of the indirect effect is approximately

where $a$ and $b$ are specified values of $\gamma_M$ and $\gamma_Y$, respectively, according to Corollaries 1–2 and Table 5. This statistic is very sensitive to $\gamma_M$ but not $\gamma_Y$ when $\gamma_M$ and $\gamma_Y$ are small.

As shown in Fig. 2b1 and Fig. 2b2, the direct effect analysis result is very sensitive to a linear combination of $\beta^{Z_T}_Y$ and $\beta^{Z_M}_Y$, and the result is more sensitive to $\beta^{Z_T}_Y$ than $\beta^{Z_M}_Y$. In fact, the estimated direct effect appears to be significant at level 0.05 when $-3.515\beta^{Z_T}_Y-1.295\beta^{Z_M}_Y>0.093$ or $-3.515\beta^{Z_T}_Y-1.295\beta^{Z_M}_Y < -0.093$, where the coefficient of $\beta^{Z_T}_Y$ is larger than that of $\beta^{Z_M}_Y$. This is because the Z-test statistics of the direct effect is approximately

according to Corollaries 3–4 and Table 5, and this statistic is very sensitive to $3.515\beta^{Z_T}_Y+1.295\beta^{Z_M}_Y$, a linear combination of $\beta^{Z_T}_Y$ and $\beta^{Z_M}_Y$.

4.   Conclusions

In this paper, we focus on sensitivity analysis for IV assumptions in MR-based mediation analysis. We use two sensitivity parameters to quantify causal indirect effect estimation bias due to violation of the IV assumptions. Similarly, we use two sensitivity parameters for the direct effect. With these sensitivity parameters, we derive consistent estimators of indirect effect and direct effect using both individual data and summary data. Numerical studies demonstrate that our method is valid with correctly specified sensitivity parameters, and it is robust to misspecifying sensitivity parameters to some extent.

Multiple genetic IVs are widely used in MR studies since they are valuable for detecting or avoiding bias if the validity assumptions are violated for several or all IVs^[28]. It is too complicated to conduct sensitivity analysis with too many IVs. Instead, we use allele score, a single variable that summarizes multiple genetic variants associated with the exposure, to reduce a SNP vector to a one-dimensional variable. Specifically, the average allele is defined as an unweighted sum of all alleles (unweighted score) or the sum of weights for each allele corresponding to estimated genetic effect sizes (weighted score).

The allele score was frequently applied in MR studies with multiple IVs. For instance, Burgess and Thompson^[29] showed that allele scores enable valid causal effect estimates with a large number of genetic variants; Pierce et al.^[30] demonstrated that employing an allele score alleviates the ‘weak-IV’ bias in the presence of confounders; and Burgess et al.^[6] utilized the weighted allele score to evaluate the indirect effect of BMI on uric acid mediated by C-reactive protein. To conduct MR-based mediation analyses for individual data, we can use the allele score to replace multiple genetic instruments and derive sensitivity parameters similar to $\gamma_M$ and $\gamma_Y$.

Sample structure is a major confounding factor that is largely ignored by existing summary-level MR methods. When sample structure is present, the estimates derived from summary data become correlated, resulting in biased estimates^[31]. To detect average causal effects with well-calibrated statistical inference, Hu et al.^[31] proposed MR-APSS to account for pleiotropy and sample structure simultaneously by leveraging genome-wide information. Our current method cannot be trivially extended to address the problem due to sample structure. Further investigation is warranted to extend MR-APSS to mediation analysis.

Conflict of interest

The authors declare that they have no conflict of interest.