Causal mediation analysis with multiple causally non-ordered and ordered mediators based on summarized genetic data

Abstract

Causal mediation analysis investigates the mechanism linking exposure and outcome. Dealing with the impact of unobserved confounders among exposure, mediator and outcome is an issue of great concern. Moreover, when multiple mediators exist, this causal pathway intertwines with other causal pathways, rendering it difficult to estimate the path-specific effects. In this study, we propose a method (PSE-MR) to identify and estimate path-specific effects of an exposure (e.g. education) on an outcome (e.g. osteoarthritis risk) through multiple causally ordered and non-ordered mediators (e.g. body mass index and pack-years of smoking) using summarized genetic data, when the sequential ignorability assumption is violated. Specifically, PSE-MR requires a specific rank condition in which the number of instrumental variables is larger than the number of mediators. Furthermore, we illustrate the utility of PSE-MR by providing guidance for practitioners and exploring the mediation effects of body mass index and pack-years of smoking in the causal pathways from education to osteoarthritis risk. Additionally, the results of simulation reveal that the causal estimates of path-specific effects are almost unbiased with good coverage and Type I error properties. Also, we summarize the least number of instrumental variables for the specific number of mediators to achieve 80% power.

Keywords

Mediation analysis multiple mediators causally ordered mediators causally non-ordered mediators Mendelian randomization summarized genetic data

1 Introduction

Mediation analyses help uncover the mechanism of an underlying causal relationship between exposure and outcome by using mediator variables.¹ In mediation analyses, the total effect of exposure on outcome is partitioned into indirect and direct effects. Indirect effects act through mediators of interest, whereas direct effect is the effect of exposure on the outcome in the absence of the mediators. Path-specific effects (PSEs) are a broad class of indirect effects of an exposure to an outcome via one or more causal pathways with respect to a subset of mediators.² Estimating PSEs using existing methods typically requires a stringent sequential ignorability assumption³ that no unmeasured confounders exist among the exposure, mediators and outcome.⁴ However, this assumption may not hold in practice and omitting important confounders will inevitably cause biases in the results.⁵ When multiple intermediate variables (e.g. M₁ and M₂) are involved in a study, three scenarios may arise with respect to M₁ and M₂, as shown in Figure 1. In Figure 1(a), M₁ is conditionally independent of M₂ given the treatment (X).⁶ In Figure 1(b), M₁ and M₂ are causally non-ordered because they are independent of each other, conditional on the exposure (X) and unmeasured confounder (U).⁷ As shown in Figure 1(c), mediators are causally ordered. M₁ is treated as a mediator-outcome confounder affected by treatment when we are interested in mediator M₂. For causally non-ordered and ordered mediators, Imai and Yamamoto⁸ proposed an approach using a linear structural equation model. Daniel et al.⁹ considered the finest possible decomposition of the total effect when there were two causally ordered mediators and evaluated the PSEs under the counterfactual framework. VanderWeele and Vansteelandt¹⁰ regarded the multiple mediators simultaneously as joint mediators and defined the ‘joint’ natural direct and indirect effects as extensions of the usual two-way decomposition of the total effect using regression-based approach and weighting approach. However, these methods require the assumption of sequential ignorability. Recently, several methods^4,11–15 have been developed to relax this assumption. We compared the assumptions of these methods in the discussion section. However, none of them allowed for the simultaneous existence of unmeasured confounders among the exposure, mediators and the outcome.

Mendelian randomization (MR) analyses¹⁶ using summarized data have recently become popular owing to the increase in the public availability of suitable data in large sample sizes from recently published genome-wide association studies (GWAS).¹⁷ For instance, Tikkanen et al.¹⁸ performed a two-sample MR to evaluate the independent causal roles of body components (fat-free mass and fat mass) on atrial fibrillation (AF). Firstly, univariate MR was used to estimate the causal effect of fat-free mass on AF by leveraging genetic variants. Certain genetic variants may be associated with both fat-free mass and fat mass, which is problematic because fat mass is associated with AF. These genetic variants are invalid because they violate the assumption of exclusion restriction, as they unlock the pathway from genetic variants to AF, not via fat-free mass. This phenomenon is called horizontal pleiotropy, and fat mass is considered a pleiotropic trait.¹⁹ To eliminate the effect of pleiotropy on causal estimation, multivariable MR²⁰ was performed to evaluate the causal role of fat-free mass on AF independent of fat mass. Burgess et al.²¹ showed that total and direct effects in a single mediator setting could be estimated using univariable and multivariable MR analyses, respectively.²¹ This method is reviewed in Section 2.1. In Section 2.2, we extend the analysis from a single mediator to multiple mediators (PSE-MR) for both causally non-ordered and ordered mediators, even in the presence of pleiotropy. In Section 3, we apply our method to estimate PSEs from education to osteoarthritis (OA) risk through body mass index (BMI) and pack-years of smoking. In Section 4, we conduct a series of simulations to evaluate the performance of the PSE-MR in different scenarios. Finally, we provide a short review to compare PSE-MR and existing mediation analyses, and then discuss the potential application of PSE-MR. The R package PSEMR for implementing PSE-MR is provided in GitHub (https://github.com/hhoulei/PSEMR).

2 Methods

2.1 Motivating example and notations

Globally the most common form of arthritis is OA. It accounts for 2.4% of all years lived with disability (YLD) and ranks as a leading contributor to global YLDs.^22,23 Some studies^23–26 have highlighted education, obesity and smoking as the common mechanisms for underlying OA. Education is a common socioeconomic factor and may influence obesity²⁷ as well as smoking.²⁸ A MR study revealed that a higher BMI increases the likelihood of becoming a smoker and increases smoking heaviness in current smokers using data from multiple cohorts.²⁹ Hence, we aim to examine whether education affects OA through its influence on BMI and pack-years of smoking in Europe (Figure 2).

Hereafter, we let X, Y, $G_{j}$ , $M_{i}$ and U denote the exposure (education), outcome (OA), jth (j = 1,…, J) genetic variant, ith (i = 1,…, n) mediator (here, we have two mediators: $M_{1}$ — BMI, $M_{2}$ — pack-years of smoking), and potential confounders among the mediators, exposure and outcome, respectively. The parameters $α_{i}$ and $δ_{i}$ denote the causal effect from exposure X to the ith mediator $M_{i}$ and the ith mediator $M_{i}$ to outcome Y, respectively. $r_{p q}$ denotes the direct effect of $M_{p}$ on $M_{q}, p, q \in (1, 2, \dots, n), p \neq q$ , that is, $r_{12}$ represents the causal effect of BMI ( $M_{1}$ ) on pack-years of smoking ( $M_{2}$ ). $θ_{0}$ denotes the direct causal effect of X (education) on Y (OA).

We define the PSEs using the counterfactual model,³⁰ and the PSE of X on Y along the path $π$ with respect to the comparison exposure value x and baseline value $x^{'}$ on the mean difference scale can be defined in terms of the difference in expectations of two nested counterfactuals

E {Y}_{p_{π} (Y)} - E {Y [x^{'}]}

where

p_{π} (Y)

is the distribution of the nested potential outcome defining the PSE along the causal path

π

. This expression can also be rewritten as the sum and product of the effects of edges along the path. Taking Figure 2 as an example, under the linearity condition, the PSE of X on Y along the path

X \to M_{1} \to M_{2} \to Y

, which corresponds to the path education→BMI→pack-years of smoking→OA in the motivation example, can be written as

E {Y [x^{'}, M_{1} (x), M_{2} (x^{'}, M_{1} (x))]} - E {Y [x^{'}]} = α_{1} r_{12} δ_{2}

The second term

E {Y [x^{'}]}

is equal to

E {Y [x^{'}, M_{1} (x^{'}), M_{2} (x^{'}, M_{1} (x^{'}))]}

. The total effect of education on OA was defined as

E {Y [x]} - E {Y [x^{'}]} = θ_{0} + α_{1} δ_{1} + α_{2} δ_{2} + α_{1} r_{12} δ_{2}

. Thus we aim to estimate the edges’ effects along the path

π

using genetic variants as instrumental variables, then PSEs and total effect can be identified and estimated.

In this study, we assume that all variables are continuous, and that the relationships between variables ( $G_{j} (j = 1, \dots, J), X, M_{i} (i = 1, \dots, n)$ ) are linear and homogeneous across the population. Then we declare that the consistency and composition assumptions in causal mediation analyses hold³¹ (refer S1 Appendix, section 1). Note that we relax the cross-world independence assumption and the assumption of no unmeasured confounders among the X (education), mediators M (BMI and pack-years of smoking), and outcome Y (OA), which are required in most mediation studies. PSE-MR analyses are implemented using GWAS summarized data (including the beta-coefficients ${\hat{β}}_{X_{j}}$ , ${\hat{β}}_{M_{i j}}$ , ${\hat{β}}_{Y_{j}}$ and their standard errors se( ${\hat{β}}_{X_{j}}$ ), se( ${\hat{β}}_{M_{i j}}$ ), se( ${\hat{β}}_{Y_{j}}$ )) on the associations of each genetic variant $G_{j} (j = 1, \dots, J)$ with the exposure X, mediators $M_{i} (i = 1, \dots, n)$ and outcome Y, obtained from linear regressions on non-overlapping data consortia. If exposure X or outcome Y is binary, such as OA in the motivating example, then these summarized association estimates can be replaced with association estimates (log(OR)) obtained from logistic regression. We assume that all genetic variants are uncorrelated (not in linkage disequilibrium) and this assumption can be relaxed using generalized weighted linear regression.³²

2.2 PSE-MR in a single mediator setting

Initially, we review the method proposed by Burgess et al.²¹ In this section, we only consider one mediator, taking $M_{1}$ (BMI) as an example, and estimate total, indirect (through $M_{1}$ (BMI), $θ_{I} = α_{1} \times δ_{1}$ ) and direct effect (not through $M_{1}$ (BMI), $θ_{D} = θ_{0}$ ) of X(education) on Y(OA) using genetic variants $G_{j} (j = 1, \dots, J)$ (Figure 1(d)).

2.2.1 PSE-MR based on IVW (PSE-IVW)

The genetic variants used to estimate the total effect of X(education) on Y(OA) must satisfy the standard assumptions of MR: they are associated with X(education), not associated with U, and there is no pathway from any $G_{j} (j = 1, \dots, J)$ to Y(OA) except via X(education) (i.e. $γ_{1 j} = γ_{0 j} = 0$ ).¹⁹ Based on these assumptions, the inverse-variance weighting method (IVW) can provide an estimate of the total effect $θ_{T}$ of X(education) on Y(OA) using the following weighted regression with the intercept set to zero

β_{Y} = θ_{T} β_{X} + ε_{T}, ε_{T} \sim N (0, Σ_{Y}), Σ_{Y} = d i a g (s e (β_{Y_{1}})^{2}, s e (β_{Y_{2}})^{2}, \dots, s e (β_{Y_{J}})^{2})

(1)

where

β_{Y} = (β_{Y_{1}}, β_{Y_{2}}, \dots, β_{Y_{J}})

and

β_{X} = (β_{X_{1}}, β_{X_{2}}, \dots, β_{X_{J}})

. The total effect

θ_{T}

of X(education) on Y(OA) can be decomposed into a direct effect

θ_{D} = θ_{0}

and an indirect effect via

M_{1}

(BMI)

θ_{I} = α_{1} \times δ_{1}

Under the framework of multivariable MR, all variants used to estimate the direct effect of X(education) on Y(OA) must satisfy the assumptions of multivariable MR: they are associated with the X(education) and/or $M_{1}$ (BMI), not associated with U, and there is no pathway from any $G_{j} (j = 1, \dots, J)$ to Y(OA) except via X(education) and/or $M_{1}$ (BMI) (that is $γ_{0 j} = 0$ ).¹⁹ The weighted regression model can be expanded by including genetic associations with the mediator

β_{Y} = θ_{0} β_{X} + δ_{1} β_{M_{1}} + ε_{D}, ε_{D} \sim N (0, Σ_{Y})

(2)

where

β_{M_{1}} = (β_{M_{11}}, β_{M_{12}}, \dots, β_{M_{1 J}})

. The indirect effect

θ_{I}

of X(education) on Y(OA) can be calculated as

{\hat{θ}}_{I D} = {\hat{θ}}_{T} - {\hat{θ}}_{D}

(difference indirect effect). It is equivalent to

{\hat{θ}}_{I P} = {\hat{α}}_{1} \times {\hat{δ}}_{1}

(product indirect effect), where

δ_{1}

can be estimated by equation (2) and

α_{1}

can be estimated by the following weighted regression with the intercept set to zero

β_{M_{1}} = α_{1} β_{X} + ε_{M_{1}}, ε_{M_{1}} \sim N (0, Σ_{M_{1}}), Σ_{M_{1}} = d i a g (s e (β_{M_{11}})^{2}, s e (β_{M_{12}})^{2}, \dots, s e (β_{M_{1 J}})^{2})

(3)

All variants used to estimate

α_{1}

must satisfy the three assumptions of MR: they are associated with X(education), not associated with U, and there is no pathway from any

G_{j} (j = 1, \dots, J)

M_{1}

(BMI) except via X(education) (i.e.

γ_{1 j} = 0

). The standard error of the difference and product indirect effects are presented in S1 Appendix, section 5. The total effect can also be estimated from individual-level data using the two-stage least squares (2SLS) method.²¹

2.2.2 PSE-MR of a single mediator based on MR-Egger (PSE-Egger)

The method proposed by Burgess et al. (2017) has some limitations. None of the effects can be unbiasedly estimated if the direct effects of $G_{j}$ on $M_{1}$ (BMI) or $G_{j}$ on Y (OA) exist (pleiotropy). Thus, we allow for direct effects of $G_{j}$ on $M_{1}$ (BMI) ( $γ_{1 j} \neq 0$ ) as well as $G_{j}$ on Y(OA) ( $γ_{0 j} \neq 0$ ). Without the limitation of the intercept set to zero, the causal effect of X(education) on Y(OA) can be obtained by the MR-Egger regression. To satisfy the InSIDE assumption³³ for MR-Egger, we require

β_{X} ⊥ γ_{0} ⊥ γ_{1},

(4)

where

γ_{0} = (γ_{01}, γ_{02}, \dots, γ_{0 J})

and

γ_{1} = (γ_{11}, γ_{12}, \dots, γ_{1 J})

. Under the assumptions mentioned in Section 2.1.1, the total, direct and indirect effects can be estimated using the following weighted linear regressions:

β_{Y} = γ_{T} + θ_{T} β_{X} + ε_{T}, ε_{T} \sim N (0, Σ_{Y}) .

(5)

β_{Y} = γ_{D} + θ_{0} β_{X} + δ_{1} β_{M_{1}} + ε_{D}, ε_{D} \sim N (0, Σ_{Y})

(6)

β_{M_{1}} = γ_{M_{1}} + α_{1} β_{X} + ε_{M_{1}}, ε_{M_{1}} \sim N (0, Σ_{M_{1}})

(7)

The intercept terms

γ_{D}

and

γ_{M_{1}}

, which differs from zero (directional horizontal pleiotropy) are indicators of the direct effects of genetic variants on the OA and genetic variants on the mediator BMI, respectively. The estimation of standard error for the difference and product indirect effects is presented in S1 Appendix.

2.3 Extending PSE-MR to multiple mediators setting

In this section, we extend the PSE-MR method to multiple mediators. If there are n mediators $M_{1}, M_{2}, \dots, M_{n}$ in the causal pathway from X to Y, PSEs can also be identified. In the multiple mediators setting, we consider two relationships among mediators: causally non-ordered and causally ordered. In both cases, total effect $θ_{T}$ of X on Y can also be estimated by equation (1). The genetic variants used to estimate $θ_{T}$ must satisfy the standard assumptions of MR, as mentioned in section 2.2.1. Notably, we can estimate PSEs from a series of regressions. First, we extend the three assumptions of MR in section 2.2.1 to the multiple mediator’s cases:

Assumption I:
For each $G_{j} (j = 1, \dots, J)$ in the kth (k = 1,…, n + 1) regression, $G_{j} ⧸ ⊥ (X, M_{1 : k - 1})$ .
Assumption II:
For each $G_{j} (j = 1, \dots, J)$ in the kth (k = 1,…, n + 1) regression, $G_{j} ⊥ U$ .
Assumption III:
For each $G_{j} (j = 1, \dots, J)$ in the kth (k = 1,…, n) regression, $G_{j} ⊥ M_{k} | (X, M_{1 : k - 1}, U)$ , for each $G_{j}$ in the (n + 1)th regression, $G_{j} ⊥ Y | (X, M_{1 : n}, U)$ ; where $M_{1 : k - 1} = {\begin{matrix} \emptyset, k = 1 \\ {M_{1}, \dots, M_{k - 1}}, k > 1 \end{matrix}$ .

Assumption I requires that $G_{j}$ should be strongly associated with X and/or $M_{i}$ , otherwise, weak instrumental variable bias will occur.³⁴ The ‘rule of thumb’ advocates that the F statistic and conditional F statistic of each instrumental variable should be at least 10 to avoid this bias^35–37 (refer S1 Appendix, section 2.3). In general, we choose the single-nucleotide polymorphisms (SNPs) associated with exposure (e.g. education) or mediators (e.g. BMI and pack-years of smoking) with a p-value less than $5 \times 10^{- 8}$ . Assumption II requires that $G_{j}$ must be independent of U. If not, the InSIDE assumption may be violated. Swanson and VanderWeele³⁸ suggested that the E-value can be used to examine the independence between $G_{j}$ and U, that is, to evaluate the sensitivity of estimates to confounders between $G_{j}$ and Y (refer S1 Appendix, section 5). A small E-value implies little unmeasured confounding would be needed to explain away an effect estimate from $G_{j}$ (SNPs) on Y(OA). Assumption III requires that there is no pleiotropy. In other words, $G_{j}$ must affect Y through X and/or $M_{i}$ , and the pathways $G_{j} \to Y$ (not via X and $M_{i}$ ) and/or $G_{j} \to M_{i}$ are not allowed. In the motivating example, it means SNPs affect the outcome (OA) only through the exposure (education) or mediators (BMI and pack-years of smoking). Egger's test examines this assumption by adding an intercept in IVW regression. If the hypothesis testing of the intercept term is significant, then pleiotropy exists. We examine and relax this assumption later. Besides, the heterogeneity test as a common sensitivity analysis in MR study tests whether there are outliers, which may violate Assumption II or III.
2.3.1 PSE-MR for causally non-ordered mediators

Firstly, we consider causally non-ordered mediators (Figure 3(a)), where n mediators are independent of each other condition on X. The causal effects of edges in the causal pathway from X to Y can be estimated using the following weighted regressions with the intercept set to zero

B = Ψ A + ε, ε \sim N (0, Σ_{B}), Σ_{B} = d i a g (s e (β_{M_{1}})^{2}, \dots, s e (β_{M_{n}})^{2}, s e (β_{Y})^{2})

(8)

where

B = [\begin{matrix} β_{M_{1}} \\ β_{M_{2}} \\ \dots \\ β_{M_{n}} \\ β_{Y} \end{matrix}], Ψ = [\begin{matrix} α_{1} & 0 & 0 & \dots & 0 \\ α_{2} & 0 & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots \\ α_{n} & \dots & \dots & \dots & \dots \\ θ_{0} & δ_{1} & δ_{2} & \dots & δ_{n} \end{matrix}], A = [\begin{matrix} β_{X} \\ β_{M_{1}} \\ β_{M_{2}} \\ \begin{matrix} \dots \\ β_{M_{n}} \end{matrix} \end{matrix}], ε = [\begin{matrix} ε_{1} \\ ε_{2} \\ ε_{3} \\ \dots \\ ε_{n + 1} \end{matrix}]

These estimations can also be obtained from individual-level data using the 2SLS method.²¹

Secondly, we relax Assumption III by allowing for the direct effect from $G_{j}$ to mediators $M_{i}$ ( $γ_{1 j}, γ_{2 j}, \dots, γ_{n j} \neq 0$ ), as well as $G_{j}$ to Y ( $γ_{0 j} \neq 0$ ) (Figure 3(a)). Under the InSIDE assumption $β_{X} ⊥ γ_{1} ⊥ γ_{2} ⊥ \dots ⊥ γ_{n} ⊥ γ_{0}$ , total effect $θ_{T}$ can also be estimated using equation (5). The causal effects of the edges in the causal pathway from X to Y can be estimated using the following linear regression equations:

B = γ + Ψ A + ε, ε \sim N (0, Σ_{B}),

(9)

where

γ = {[\begin{matrix} γ_{1} & γ_{2} & \dots & γ_{n} & γ_{0} \end{matrix}]}^{T}

. Intercept terms

γ_{0}

and

γ_{i} (i = 1, \dots, n)

that differ from zero are indicators of direct effects between genetic variants and outcome, as well as genetic variants and mediators, respectively. In addition, to obtain a consistent estimation, we require that the number of SNPs to be larger than the number of mediators. Detailed theoretical derivations are presented in S1 Appendix, section 3.

2.3.2 PSE-MR for causally ordered mediators

When all mediators are causally ordered (Figure 3(b)), the causal effects of edges in the causal pathway from X to Y can be estimated by the weighted regressions in equations (8) and (9) by substituting $Ψ *$ for $Ψ$ , where

Ψ * = [\begin{matrix} α_{1} & 0 & 0 & \dots & 0 \\ α_{2} & r_{12} & 0 & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots \\ α_{n} & r_{1 n} & r_{2 n} & \dots & 0 \\ θ_{0} & δ_{1} & δ_{2} & \dots & δ_{n} \end{matrix}] .

The causal effect

r_{p q}

from

M_{p}

M_{q}, p, q \in (1, 2, \dots, n), p \neq q

can also be identified. Details of the theoretical derivation are presented in S1 Appendix, section 3.

In this case, the cross-world independence assumption is violated if pleiotropy exists. This is because, taking two ordered mediators as an example (Figure 1(f)), $M_{1}$ (BMI) acts as the X-induced confounder between $M_{2}$ (pack-years of smoking) and Y(OA) when we are interested in the mediator $M_{2}$ . If $M_{1}$ (BMI) is missed, the InSIDE assumption would be $β_{X} ⊥ (r_{12} γ_{1} + γ_{2}) ⊥ (γ_{0} + δ_{1} γ_{1})$ , which may not be satisfied because of the common factor $γ_{1}$ except when $δ_{1}$ or $r_{12}$ was zero. If we were missing $M_{2}$ (pack-years of smoking), the estimation of new direct and indirect effects were unbiased regardless of the existence of pleiotropy. The InSIDE assumption $β_{X} ⊥ γ_{1} ⊥ (γ_{0} + δ_{2} γ_{2})$ is satisfied. Therefore, PSE-MR allows the mediator to be missed if it doesn’t have direct causal effect on Y or downstream mediators.

If interactions exist, the linearity assumption is violated. When there are interactions between exposure and mediators on the outcome, our method is infeasible for summary data but is still available for individual data. We can divide the entire population into several subgroups according to mediators from small to large and estimate the PSE in different subgroups from equations (8) and (9), respectively. Because the causal effect from exposure to outcome is not constant at different values of mediators (refer S1 Appendix, section 4.1). Similarly, when there are interactions between mediators (refer S1 Appendix, section 4.4), our method is also available for individual data and PSEs can be estimated by dividing the entire population into multiple subgroups according to the different values of mediators. Specifically, when there are interactions between mediators (BMI and pack-years of smoking) and confounders (refer S1 Appendix, section 4.3), or interactions between exposure and unmeasured confounders (refer S1 Appendix, section 4.2), not all the PSEs of exposure (education) on the outcome (OA) can be identified because of the uncertainty of unmeasured confounders. For the former, the direct effects of exposure on mediators and outcomes can be identified. For the latter, direct effects from mediators to outcome can be identified. We illustrate this in detail in S1 Appendix, section 4.

3 Application

As an illustrative example, we attempt to reveal the causal mechanism from education to OA. To clearly illustrate this example, a guide for practitioners was prepared to apply PSE-MR (Figure 4). The use of this pipeline to perform PSE-MR is demonstrated in the practical example given below.

For step A, our analysis focused on the diagram in Figure 2, which was obtained by the prior knowledge mentioned at the beginning of Section 2.1. Further, we prepared the summary data for education, BMI, pack-years of smoking and OA (step B): (1) genetic associations with education (years of schooling) in 1.1 million participants from European were obtained from the Social Science Genetic Association Consortium (SSGAC)³⁹; (2) genetic associations with BMI in 315,347 participants were obtained from the single large multiethnic Genetic Epidemiology Research on Adult Health and Aging (GERA) cohort⁴⁰; (3) genetic associations with pack-years of smoking and OA in 142,387 and 417,596 participants (39,427 cases and 378,169 controls) were obtained from the UK Biobank cohort, respectively.⁴¹ Summary data of binary OA included SNP information, log(OR) and their standard errors from logistic regression (step D). Summary data of continuous phenotypes (education, BMI and pack-years of smoking) included SNP information, beta-coefficients and their standard errors from linear regression (step E). To ensure the homogeneity of the population (step F), all populations from the above database were mainly European ancestry.

Then we chose SNPs as valid instrumental variables for PSE-MR. These SNPs were filtered by the P-value ( $P \leq 5 \times 10^{- 8}$ ) and linkage disequilibrium (LD) ( $r^{2} < 0.1$ , clumping window = 1000 kb) (step G). Finally, we identified 470, 382 and 527 SNPs associated with education, for single mediator (BMI, pack-years of smoking) and two mediators analyses, respectively. The number of instrumental variables was far greater than the number of mediators, which indicated that the rank condition of PSE-MR was satisfied (step H).

Next, we performed several sensitivity analyses to examine whether the three core assumptions of PSE-MR were satisfied. For Assumption I, we calculated F statistics for genetic instruments and F statistics were greater than 10 (steps J and K). SNPs with F statistics smaller than 10 cannot be instrumental variables because they can induce weak instrumental bias. E-values between SNPs and OA were small enough to satisfy Assumption II (steps L and M, refer S1 Appendix, section 6, eTables 1 to 3). SNPs with large E-values must be removed because the InSIDE assumption may be violated. Heterogeneity tests revealed that there was strong heterogeneity for the Wald ratios⁴² of education − BMI, education − pack-years of smoking and education-OA (step N, refer S1 Appendix, section 6, eTables 4, 6, 8). Following step O, we should remove outliers by MR-PRESSO,⁴³ 297, 315 and 306 SNPs were left for single mediator (BMI, pack-years of smoking) and two mediators analyses, respectively (refer S1 Appendix, section 6, eTables 5, 7, 9). In addition, we evaluated whether horizontal pleiotropy was present using the Egger's test (step P). The results revealed no significant effects of BMI, pack-years of smoking, or OA (refer S1 Appendix, section 6, eTable 10), indicating the absence of directional pleiotropy. Thus, we have been well prepared to perform PSE-MR, and PSE-IVW as the primary method. If Egger's test revealed significant horizontal pleiotropy, we chose PSE-Egger as the primary method.

We performed single and multiple mediators analyses to explore the causal pathway from education to OA via PSE-MR. Figure 5 shows the path-specific effects of education on OA (step Q). Single mediator analysis suggested that BMI and pack-years of smoking were mediators in the causal pathway from education to OA. Then we performed PSE-MR analysis with multiple mediators to test whether education had indirect effects on OA risk through BMI and pack-years of smoking. We found that education was a protective factor against the risk of OA and a significant direct effect was obtained after adjusting for genetic associations with BMI and pack-years of smoking. Indirect effects through BMI and pack-years of smoking explained a considerable proportion of the causal effect of education on OA, and their total mediation proportion (MP) was 53.97%. In conclusion, three mediated pathways exist from education to OA: education→BMI→OA (MP: 21.52%), education→pack-years of smoking→OA (MP: 38.66%), and education→BMI→pack-years of smoking→OA (MP: 1.68%). These results (Figure 5) are consistent with the results of other observational and MR studies^23,44 and previously described biological mechanisms.⁴⁵

4 Simulation

4.1 Settings

To validate the utility of the PSE-MR method for estimating PSEs, we considered the multiple mediators with causally non-ordered ( $r_{12} = 0$ ) and ordered mediators ( $r_{12} = 0.8$ ). For ith individual (i = 1, 2, ..., N), The data-generating models were

\begin{aligned} U_{i} = \sum_{j = 1}^{J} κ_{j} G_{i j} + ε_{U_{i}} \\ X_{i} = \sum_{j = 1}^{J} ϖ_{j} G_{i j} + U_{i} + ε_{X_{i}} \\ M_{2 i} = α_{2} X_{i} + U_{i} + \sum_{j = 1}^{J} γ_{2 j} G_{i j} + ε_{M_{2 i}} \\ M_{1 i} = α_{1} X_{i} + r_{12} M_{2 i} + U_{i} + \sum_{j = 1}^{J} γ_{1 j} G_{i j} + ε_{M_{1 i}} \\ Y_{i} = θ_{0} X_{i} + δ_{1} M_{1 i} + δ_{2} M_{2 i} + U_{i} + \sum_{j = 1}^{J} γ_{0 j} G_{i j} + ε_{Y_{i}} \end{aligned}

where the error variables

ε_{X_{i}}, ε_{M_{1 i}}, ε_{M_{2 i}}, ε_{Y_{i}}

were independently generated from an

N (0, 1)

distribution; and the instrument strength parameters

ϖ_{j}

were generated from an

U (0.05, 2)

distribution.

γ_{0 j}, γ_{1 j}

and

γ_{2 j}

represent the direct effects of genetic variants on the outcome,

M_{1}

and

M_{2}

, respectively.

We classified all the genetic variants as $G_{j}$ to five categories: (1) $G_{X}$ : only associated with X; (2) $G_{M_{1}}$ : only associated with $M_{1}$ ; (3) $G_{M_{2}}$ : only associated with $M_{2}$ ; (4) $G_{X M_{1}}$ : associated with X and $M_{1}$ ; (5) $G_{X M_{12}}$ : associated with X, $M_{1}$ and $M_{2}$ . We generated data on 30 genetic variants for category (1), 10 genetic variants for categories (2) to (5), an exposure (X), two mediators ( $M_{1}$ and $M_{2}$ ), and an outcome (Y) for 500,000 individuals (UK Biobank scale sample size) in 100 repeated datasets. Firstly, we observed the performance of PSE-MR with a different number of instrumental variables and sample sizes. Secondly, according to the specification of direct and indirect effects, all simulations included four scenarios: no direct effect, no indirect effect, a direct effect and a directionally concordant indirect effect, as well as a direct effect and a directionally discordant indirect effect. Thus we have 10 and 11 parameter settings for the analysis of two non-ordered and ordered mediators, respectively (refer S2 Appendix, sections 3 to 6). Thirdly, for pleiotropy, the data were simulated to consider the following three cases (detailed settings for $γ_{0 j}, γ_{1 j}$ and $γ_{2 j}$ are shown in S2 Appendix, eTables 8 and 16):

Case (a): Balanced pleiotropy, InSIDE assumption satisfied.

Case (b): Directional pleiotropy, InSIDE assumption satisfied.

Case (c): Directional pleiotropy, InSIDE assumption not satisfied.

We believe that the PSE-IVW performs well in Case (a). As a solution to consider pleiotropy, the PSE-Egger can also perform well in Case (b). We consider how much of an impact it has on the estimate when the InSIDE assumption is violated (Case(c)). Finally, we also performed additional simulations for sensitivity analyses, where the population homogeneity assumption is violated, the causal order is misspecified, one of the mediators is missing, bi-directional causal effects between exposure and mediators exist, exposure and mediators are time-varying. PSE-MR may not be robust in these scenarios because it may induce necessary assumptions of PSE-MR violation. Additionally, we determined the optimal number of genetic variants when we considering multiple mediators. The details of the simulation are presented in S2 Appendix.

We used the following metrics to evaluate the performance of our methods: relative bias, mean square error (MSE), coverage ratio, Type I error rate for a null causal effect and empirical power to detect a non-null effect (i.e. the proportion of confidence intervals excluding zero). Standard errors of the PSEs were calculated by bootstrap.

4.2 Results

We found that the causal estimates of PSEs were unbiased with good Type I error properties. As the sample size increased, bias and standard errors decreased, while power improved. Higher power and lower bias were observed as the number of instrumental variables increased (refer S2 Appendix, sections 1, 3 and 5, eTables 1, 3, 6, 7, 13 and 14).

For the two non-ordered mediators, PSE-IVW showed good performance of in standard MR when estimating the total, direct and indirect effects as well as the three PSEs (Table 1). As the sample size and the number of genetic variants increased, the bias was reduced, and the Type I error was more stable at approximately 0.05 (refer S2 Appendix, section 3). The performance of PSE-MR based on IVW and MR-Egger with two non-ordered mediators in Cases (a) and (b), are listed in eTables 9 to12 (refer S2 Appendix, section 4). In Case (a), we observed that the bias was close to zero and the Type I error rates were approximately 0.05 in PSE-MR. PSE-Egger had less bias and more stable Type I error rates than IVW when directional pleiotropy existed in at least one pathway from G to Y (Case (b)). MR-Egger performed better than IVW in terms of bias, even when the InSIDE assumption was not satisfied (Case (c)). When the pleiotropic effects through confounders (violating the InSIDE assumption) were 2.5 times larger than the direct pleiotropic effects (satisfying the InSIDE assumption), estimates from PSE-Egger were much less biased and rejection rates of the causal null hypothesis were much closer to the nominal 5% rate than those from PSE-IVW were. In all cases, PSE-Egger had a smaller MSE and a more stable Type I error rate (0.05) than PSE-IVW when the PSE was zero. Estimators of indirect effects based on the product method had more stable Type I error rates (0.05) than those based on the difference method. The results for the two ordered multiple mediators were similar to those of the two non-ordered mediators (Table 2 and eTables 17 to 24 in S2 Appendix, section 6). In addition, the magnitude of $r_{q p}$ did not influence the performance of PSE-MR. The details are presented in eTable 15 (refer S2 Appendix, section 5).

Table 1.
Simulation results of PSE-IVW with two non-ordered mediators in standard MR.

True Relative bias MSE Power/Type I error rate Coverage ratio

TE DE IE IE1 IE2 TE DE IE_d IE_p IE1 IE2 TE DE IE_d IE_p IE1 IE2 TE DE IE_d IE_p IE1 IE2 TE DE IE_d IE_p IE1 IE2

1.6 0.8 0.8 0.5 0.3 0.012 −0.004 0.027 0.027 0.016 0.034 0.003 0.001 0.003 0.002 0.001 0.001 1.000 1.000 1.000 1.000 1.000 1.000 0.953 0.973 0.933 0.960 0.947 0.947

0.6 0.8 −0.2 −0.5 0.3 0.045 0.014 −0.080 −0.086 −0.024 0.007 0.003 0.001 0.002 0.002 0.001 0.001 1.000 1.000 0.973 1.000 1.000 1.000 0.967 0.927 0.953 0.953 0.967 0.973

0.7 0.5 0.2 0.5 −0.3 0.013 0.001 0.044 0.059 0.002 −0.047 0.003 0.001 0.003 0.002 0.001 0.001 1.000 1.000 0.993 1.000 1.000 1.000 0.973 0.953 0.973 0.953 0.940 0.867

−0.3 0.5 −0.8 −0.5 −0.3 −0.082 0.028 −0.013 −0.016 −0.013 −0.029 0.004 0.000 0.003 0.002 0.001 0.001 1.000 1.000 1.000 1.000 1.000 1.000 0.920 0.900 0.913 0.960 0.967 0.920

0 −0.8 0.8 0.5 0.3 0.027 0.004 0.038 0.034 0.027 0.032 0.004 0.001 0.003 0.002 0.001 0.001 0.133 1.000 1.000 1.000 1.000 1.000 0.927 0.927 0.933 0.933 0.960 0.927

−1.6 −0.8 −0.8 −0.5 −0.3 −0.018 −0.019 −0.018 −0.016 −0.016 −0.027 0.003 0.000 0.002 0.002 0.001 0.001 1.000 1.000 1.000 1.000 1.000 1.000 0.973 0.880 0.973 0.960 0.940 0.933

0.8 0 0.8 0.5 0.3 0.027 0.000 0.026 0.027 0.019 0.029 0.004 0.001 0.003 0.002 0.001 0.001 1.000 0.053 1.000 1.000 1.000 1.000 0.920 0.940 0.900 0.853 0.947 0.893

1.3 0.8 0.5 0.5 0 0.007 −0.005 0.026 0.031 0.009 0.011 0.004 0.000 0.003 0.002 0.001 0.001 1.000 1.000 1.000 1.000 0.100 1.000 0.920 0.973 0.907 0.920 0.887 0.900

0.8 0.8 0 0 0 0.017 0.004 0.010 0.011 0.007 0.004 0.002 0.000 0.001 0.001 0.001 0.000 1.000 1.000 0.120 0.153 0.093 0.113 0.947 0.987 0.920 0.920 0.933 0.933

0.8 0.8 0 0 0 0.031 0.009 0.018 0.017 0.007 0.010 0.003 0.000 0.002 0.002 0.001 0.001 1.000 1.000 0.080 0.160 0.107 0.173 0.953 0.960 0.953 0.967 0.927 0.947

True	Relative bias	MSE	Power/Type I error rate	Coverage ratio
1.6	0.8	0.8	0.5	0.3	0.012	−0.004	0.027	0.027	0.016	0.034	0.003	0.001	0.003	0.002	0.001	0.001	1.000	1.000	1.000	1.000	1.000	1.000	0.953	0.973	0.933	0.960	0.947	0.947
0.6	0.8	−0.2	−0.5	0.3	0.045	0.014	−0.080	−0.086	−0.024	0.007	0.003	0.001	0.002	0.002	0.001	0.001	1.000	1.000	0.973	1.000	1.000	1.000	0.967	0.927	0.953	0.953	0.967	0.973
0.7	0.5	0.2	0.5	−0.3	0.013	0.001	0.044	0.059	0.002	−0.047	0.003	0.001	0.003	0.002	0.001	0.001	1.000	1.000	0.993	1.000	1.000	1.000	0.973	0.953	0.973	0.953	0.940	0.867
−0.3	0.5	−0.8	−0.5	−0.3	−0.082	0.028	−0.013	−0.016	−0.013	−0.029	0.004	0.000	0.003	0.002	0.001	0.001	1.000	1.000	1.000	1.000	1.000	1.000	0.920	0.900	0.913	0.960	0.967	0.920
0	−0.8	0.8	0.5	0.3	0.027	0.004	0.038	0.034	0.027	0.032	0.004	0.001	0.003	0.002	0.001	0.001	0.133	1.000	1.000	1.000	1.000	1.000	0.927	0.927	0.933	0.933	0.960	0.927
−1.6	−0.8	−0.8	−0.5	−0.3	−0.018	−0.019	−0.018	−0.016	−0.016	−0.027	0.003	0.000	0.002	0.002	0.001	0.001	1.000	1.000	1.000	1.000	1.000	1.000	0.973	0.880	0.973	0.960	0.940	0.933
0.8	0	0.8	0.5	0.3	0.027	0.000	0.026	0.027	0.019	0.029	0.004	0.001	0.003	0.002	0.001	0.001	1.000	0.053	1.000	1.000	1.000	1.000	0.920	0.940	0.900	0.853	0.947	0.893
1.3	0.8	0.5	0.5	0	0.007	−0.005	0.026	0.031	0.009	0.011	0.004	0.000	0.003	0.002	0.001	0.001	1.000	1.000	1.000	1.000	0.100	1.000	0.920	0.973	0.907	0.920	0.887	0.900
0.8	0.8	0	0	0	0.017	0.004	0.010	0.011	0.007	0.004	0.002	0.000	0.001	0.001	0.001	0.000	1.000	1.000	0.120	0.153	0.093	0.113	0.947	0.987	0.920	0.920	0.933	0.933
0.8	0.8	0	0	0	0.031	0.009	0.018	0.017	0.007	0.010	0.003	0.000	0.002	0.002	0.001	0.001	1.000	1.000	0.080	0.160	0.107	0.173	0.953	0.960	0.953	0.967	0.927	0.947

TE: total effect; DE: direct effect; IE: indirect effect; IE1: X→M₁→Y; IE2: X→M₂→Y; IE_d: indirect effect calculated by difference method; IE_p: indirect effect calculated by product method.

Table 2.

Simulation results of PSE-IVW with two ordered mediators in standard MR.

True						Relative bias							MSE							Power/Type I error rate							Coverage ratio
TE	DE	IE	IE1	IE2	IE3	TE	DE	IE_d	IE_p	IE1	IE2	IE3	TE	DE	IE_d	IE_p	IE1	IE2	IE3	TE	DE	IE_d	IE_p	IE1	IE2	IE3	TE	DE	IE_d	IE_p	IE1	IE2	IE3
1.84	0.8	1.04	0.5	0.3	0.24	0.027	0.002	0.046	0.046	0.026	0.053	0.073	0.006	0.001	0.006	0.005	0.001	0.001	0.001	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.980	0.960	0.973	0.953	0.980	0.920	0.967
0.64	0.6	0.04	−0.5	0.3	0.24	0.021	0.015	0.109	0.009	−0.006	−0.014	0.042	0.005	0.000	0.004	0.003	0.001	0.001	0.000	1.000	1.000	0.160	0.167	1.000	1.000	1.000	0.947	0.953	0.973	0.980	0.973	0.953	0.940
0.76	0.8	−0.04	0.5	−0.3	−0.24	0.038	0.001	−0.070	0.030	0.019	−0.068	0.001	0.005	0.000	0.004	0.004	0.001	0.001	0.001	1.000	1.000	0.053	0.033	1.000	1.000	1.000	0.927	0.980	0.927	0.927	0.960	0.873	0.933
−0.24	0.8	−1.04	−0.5	−0.3	−0.24	−0.131	0.017	−0.017	0.015	−0.014	−0.029	0.002	0.006	0.000	0.004	0.004	0.001	0.001	0.000	0.820	1.000	1.000	1.000	1.000	1.000	1.000	0.920	0.920	0.947	0.920	0.953	0.940	0.953
0.24	−0.8	1.04	0.5	0.3	0.24	0.119	0.001	0.028	0.030	0.017	0.031	0.050	0.005	0.001	0.004	0.003	0.001	0.001	0.001	0.987	1.000	1.000	1.000	1.000	1.000	1.000	0.927	0.933	0.920	0.940	0.967	0.940	0.940
−1.84	−0.8	−1.04	−0.5	−0.3	−0.24	−0.015	−0.019	−0.013	0.014	−0.015	−0.018	−0.003	0.005	0.000	0.004	0.003	0.001	0.001	0.000	1.000	1.000	1.000	1.000	1.000	1.000	1.000	0.953	0.887	0.967	0.960	0.960	0.940	0.947
1.04	0	1.04	0.5	0.3	0.24	0.025	−0.004	0.028	0.032	0.014	0.042	0.053	0.006	0.001	0.006	0.004	0.001	0.001	0.001	1.000	0.027	1.000	1.000	1.000	1.000	1.000	0.933	0.973	0.893	0.893	0.907	0.940	0.927
1.34	0.8	0.54	0	0.3	0.24	0.017	0.003	0.038	0.023	0.006	0.017	0.051	0.006	0.000	0.004	0.003	0.001	0.001	0.001	1.000	1.000	1.000	1.000	1.000	0.060	1.000	0.960	0.947	0.953	0.953	0.973	0.947	0.947
0.8	0.8	0	0	0	0	0.019	0.001	0.014	0.012	0.009	0.003	0.000	0.001	0.000	0.001	0.001	0.001	0.000	0.000	1.000	1.000	0.093	0.127	0.060	0.100	0.007	0.967	0.973	0.947	0.973	0.967	0.960	1.000
0.8	0.8	0	0	0	0	0.042	0.006	0.029	0.028	0.009	0.011	0.007	0.006	0.000	0.005	0.004	0.001	0.001	0.000	1.000	1.000	0.180	0.173	0.140	0.127	0.140	0.973	0.960	0.933	0.927	0.920	0.907	0.907
0.24	−0.8	1.04	0.5	0.3	0.24	0.101	0.005	0.027	0.028	0.013	0.029	0.053	0.006	0.000	0.005	0.003	0.001	0.001	0.001	0.973	1.000	1.000	1.000	1.000	1.000	1.000	0.967	0.980	0.967	0.960	0.973	0.927	0.933

TE: total effect; DE: direct effect; IE: indirect effect; IE1: X→M₁→Y; IE2: X→M₂→Y; IE3: X→M₁→M₂→Y; IE_d: indirect effect calculated by difference method; IE_p: indirect effect calculated by product method.

The estimation of the direct effect is unbiased regardless of whether bidirectional causal effects between exposure and mediators exist, or the causal order is misspecified, although the estimation of PSEs is biased. Heterogeneous populations sometimes introduce bias in causal estimation for non-ordered and ordered mediators. If upstream mediators (e.g. M₁) are missing, M₁ is the confounder of M₂ and Y and is affected by X (i.e. X – induced unmeasured confounder of M₂ and Y). Thus the assumption of cross-world independence is violated. In addition, if we can obtain information at each time point, PSE-MR can be applied to time-varying exposure and mediators and it can also deal with the bi-directional relationship between exposure and mediators (refer S1 Appendix, sections 7 to 13). We also summarized the least number of instrumental variables for the specific number of mediators to achieve 80% power (Figure 6). The details are listed in eTable 41 to 42.

5 Discussion

In this study, we develop a method PSE-MR to identify and estimate PSEs from an exposure on an outcome through mediator(s) using MR when there are unmeasured confounders among the exposure, mediators and the outcome. We extend PSE-MR from a single mediator setting to multiple mediators setting for both causally ordered and non-ordered mediators, and outline the assumptions required to obtain a causal effect. PSE-IVW can be used to explore the role of multiple mediators in causal pathways between the exposure and outcome. The PSE-Egger can be viewed as a sensitivity analysis to provide robustness against both measured and unmeasured pleiotropy and to strengthen the evidence from the PSE-IVW analysis.

We made a short review to compare the assumptions of PSE-MR with traditional mediation analysis methods in Table 3. Almost all mediation methods require consistency and composition assumptions. If the cross-world independence assumption is violated, semi-parametric method¹¹ is a viable choice for performing mediation analysis. If we know the confounders that are affected by exposure, effect decomposition proposed by VanderWeele et al.⁴⁶ and Mittinty et al.,⁴⁷ and PSE-MR can be regarded as appropriate solutions. VanderWeele and Vansteelandt⁴⁸ and Valeri and VanderWeele⁴⁹ provided illustrations and expressions of direct and indirect effects in the presence of additive interaction of the exposure and mediators on the outcome. Recently, VanderWeele⁵⁰ proposes a four-way decomposition that, also allowed the existence of interaction. For PSE-MR, factorial MR⁴² can be used to estimate causal effects when there are interactions of multiple exposures on the mediators if there are individual data. Notably, multicollinearity among the genetic associations of exposure and mediators may affect the estimation. Therefore, semi-parametric¹¹ and non-parametric⁵¹ methods were proposed to avoid strict linear conditions.

Table 3.

Comparison of the assumption in PSE-MR and typical causal mediation analysis.

	PSE-MR	Existing causal mediation analysis	Solution for violating assumptions
Similarities	(1) Consistency assumption: M(x) = M and Y(x) = Y if X = x; Y(x,m) = Y if X = x, M = m.
	(2) Composition assumption: Y(x) = Y(x,M(x)) if X = x.
	(3) A cross-world independence assumption: $Y (x, m) ⊥ M (x ) \| C$ for all (x, x, m).		^11,46,47,52
	(4) There is no additive interaction of the exposure and the mediators on the outcome (Y).		^11,48–50
	(5) Linearity.		^11,51
Differences	(1) The IV $G_{j}$ is associated with the exposure X but independent of all the unmeasured confounders.	(1) No-unmeasured confounders of the X-Y relation, that is, $Y (x, m) ⊥ X \| C$ for all (x, m).	^4,11
	(2) The Exclusion restriction assumption or the InSIDE assumption.	(2) No-unmeasured confounders of the M-Y relation, that is, $Y (x, m) ⊥ M (x) \| X = x, C$ for all (x,m).	^2,4,11–15
	(3) There is no additive interaction of the exposure or the mediators and confounders on the mediator (M) and the outcome (Y).	(3) No-unmeasured confounders of the X-M relation, that is, $M (x) ⊥ X \| C$ for all x.	^4,11
	(4) Data for exposure, mediators and outcome can from different datasets with homogeneous population.	(4) Data for exposure, mediators and outcome must from the same dataset.
		(5) Exposure, mediators and confounders not vary with time.	^52–54

Several methods have been used for time-varying mediation analyses. Vansteelandt et al.⁵³ provided a new method to perform a mediation analysis of time-to-event endpoints accounting for repeatedly measured mediators subject to time-varying confounding. Zheng and van der Laan⁵⁴ proposed a longitudinal mediation analysis with time-varying mediators and exposure with application to survival outcomes. PSE-MR can violate the stringent sequential ignorability assumption. Although several methods have attempted to relax this assumption, only some of these assumptions can be relaxed. For example, many methods only allowed for the existence of confounders between mediators and outcome, and assumed that exposure was randomized,^2,12–14 which was difficult to satisfy in observational data; even semi-parametric methods^4,11 did not allow for the existence of three types of confounders. All these methods require individual data, except for PSE-MR, which can be performed if only summary data are available. This is an advantage of the PSE-MR.

PSE-MR can estimate the direct effect between exposure and outcome and the indirect effects through mediators when the sequential ignorability assumption¹¹ in mediation analyses is relaxed. The proposed method requires other independence assumptions. While Assumptions I and III are testable, there is no accepted method to test Assumption II. Several sensitivity analyses can be performed to examine this assumption, such as E-value^38,55 and heterogeneity tests. The validity of the multiple mediator PSE-Egger and its ability to estimate consistent causal effects rely on the InSIDE assumption^20,56 being satisfied. When the direct genetic associations with exposure are independent of the direct genetic associations with mediators and outcome, the InSIDE assumption is satisfied. Whereas the InSIDE assumption is plausible in some cases, it will not always be valid. For example, heterogeneous populations and misspecification of multiple mediators could bias the mediation effect estimation. When $γ_{k j}$ is not independent of each other or $γ_{0 j}$ is not independent of $γ_{k j}$ for k = 1,…,n (e.g. we are missing one of the multiple mediators), the direct effect is downward-biased and the indirect effect is upward-biased. According to our simulation, we found that PSE-IVW is more robust in estimating causal effects than PSE-Egger for heterogeneous populations and misspecified multiple mediators. However, PSE-Egger can be applied to test directional pleiotropy, and can provide less biased estimates when the InSIDE assumption is violated.

For multiple causally ordered mediator settings, PSE-MR can be widely used for time-varying exposures and mediators. Labrecque and Swanson⁵⁷ suggested that if the genetic associations of the exposure and mediators were time-varying, the lifetime effect estimate could be biased if information on the exposure and mediators was obtained only at one time point. However, if we can obtain the information on the exposure and mediators at different time points, PSE-MR can provide unbiased estimates of the lifetime effects of exposure and mediators on the outcome and other PSEs (refer S2 Appendix, section 9). Thus PSE-MR can estimate each PSE, including causal relationships (which may potentially be bi-directional) in a non-experimental setting.

In conclusion, we propose a method of causal mediation analysis with causally ordered and non-ordered mediators based on summarized genetic data and provide a new perspective for mediation analysis.

Supplemental Material

sj-docx-1-smm-10.1177_09622802221084599 - Supplemental material for Causal mediation analysis with multiple causally non-ordered and ordered mediators based on summarized genetic data

Supplemental material, sj-docx-1-smm-10.1177_09622802221084599 for Causal mediation analysis with multiple causally non-ordered and ordered mediators based on summarized genetic data by Lei Hou, Yuanyuan Yu, Xiaoru Sun, Xinhui Liu, Yifan Yu, Hongkai Li and Fuzhong Xue in Statistical Methods in Medical Research

Supplemental Material

sj-docx-2-smm-10.1177_09622802221084599 - Supplemental material for Causal mediation analysis with multiple causally non-ordered and ordered mediators based on summarized genetic data

Supplemental material, sj-docx-2-smm-10.1177_09622802221084599 for Causal mediation analysis with multiple causally non-ordered and ordered mediators based on summarized genetic data by Lei Hou, Yuanyuan Yu, Xiaoru Sun, Xinhui Liu, Yifan Yu, Hongkai Li and Fuzhong Xue in Statistical Methods in Medical Research

Footnotes

Acknowledgements

We would like to thank Editage () for English language editing.

Authors’ contributions

HL and FX conceived the study. LH, HL contributed to theoretical derivation with assistance from YY, XL and XS. LH and YY contributed to the data simulation. LH and HL contributed to the application. LH and HL wrote the manuscript with input from all other authors. All authors reviewed and approved the final manuscript.

Ethics approval and consent to participate

Ethical approval was not sought, because this study involved analysis of publicly available summary-level data from GWASs, and no individual-level data were used.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

HL was supported by the National Natural Science Foundation of China (Grant 82003557). FX was supported by the National Natural Science Foundation of China (Grant 82173625) and the Shandong Provincial Key Research and Development project (2018CXGC1210).

Availability of data and materials

GWAS summary data for education, body mass index, pack-years of smoking, and osteoarthritis are publicly available at the MR-BASE platform (). Code to implement the method and reproduce all simulations and analyses is available on Github (https://github.com/hhoulei/PSEMR).

Supplemental Material

Supplemental material for this paper is available online.

ORCID iDs

Hongkai Li

Fuzhong Xue

References

Lee

Herbert

McAuley

. Mediation analysis. JAMA 2019; 321: 697–698.

Miles

BCH

Shpitser

Kanki

, et al. On semiparametric estimation of a path-specific effect in the presence of mediator-outcome confounding. Biometrika 2020; 107: 159–172.

Imai

Keele

Yamamoto

. Identification, inference and sensitivity analysis for causal mediation effects. Stat Sci 2010; 25: 51–71.

Fulcher

Shi

Tchetgen Tchetgen

. Estimation of natural indirect effects robust to unmeasured confounding and mediator measurement error. Epidemiology 2019; 30: 825–834.

Mccandless

Somers

. Bayesian sensitivity analysis for unmeasured confounding in causal mediation analysis. Stat Methods Med Res 2019; 28: 515–531.

MacKinnon

. Contrasts in multiple mediator models. Multivariate applications in substance use research: New methods for new questions, 141–160. 2000.

Avin

Shpitser

Pearl

. Identifiability of path-specific effects. 2005.

Imai

Yamamoto

. Identification and sensitivity analysis for multiple causal mechanisms: Revisiting evidence from framing experiments. Polit Anal 2013; 21: 141–171.

Daniel

De Stavola

Cousens

, et al. Causal mediation analysis with multiple mediators. Biometrics 2015; 71: 1–14.

10.

VanderWeele

Vansteelandt

. Mediation analysis with multiple mediators. Epidemiol Method 2014; 2: 95–115.

11.

Tchetgen

EJT

Shpitser

. Semiparametric theory for causal mediation analysis: Efficiency bounds, multiple robustness, and sensitivity analysis. Ann Stat 2012; 40: 1816–1845.

12.

Luo

Geng

. Causal mediation analysis for survival outcome with unobserved mediator–outcome confounders. Comput Stat Data Anal 2016; 93: 336–347.

13.

VanderWeele

Chiba

. Sensitivity analysis for direct and indirect effects in the presence of exposure-induced mediator-outcome confounders. Epidemiol Biostat Public Health 2014; 11: e9027.

14.

Smith

VanderWeele

. Mediational E-values: Approximate sensitivity analysis for unmeasured mediator–outcome confounding. Epidemiology 2019; 30: 835–837.

15.

Ding

Vanderweele

. Sharp sensitivity bounds for mediation under unmeasured mediator-outcome confounding. Biometrika. 2016;103:483–490.

16.

Burgess

Small

Thompson

. A review of instrumental variable estimators for Mendelian randomization. Stat Methods Med Res 2017; 26: 2333–2355.

17.

Burgess

Butterworth

Thompson

. Mendelian randomization analysis with multiple genetic variants using summarized data. Genet Epidemiol 2013; 37: 658–665.

18.

Tikkanen

Gustafsson

Knowles

, et al. Body composition and atrial fibrillation: A Mendelian randomization study. Eur Heart J 2019; 40: 1277–1282.

19.

Burgess

Thompson

. Multivariable Mendelian randomization: The use of pleiotropic genetic variants to estimate causal effects. Am J Epidemiol 2015; 181: 251–260.

20.

Rees

Wood

Burgess

. Extending the MR-Egger method for multivariable Mendelian randomization to correct for both measured and unmeasured pleiotropy. Stat Med 2017; 36: 4705–4718.

21.

Burgess

Thompson

Rees

, et al. Dissecting causal pathways using Mendelian randomization with summarized genetic data: application to age at menarche and risk of breast cancer. Genetics 2017; 207: 481–487.

22.

Hawker

. Osteoarthritis is a serious disease. Clin Exp Rheumatol 2019; 37: 3–6.

23.

Gill

Karhunen

Malik

, et al. Cardiometabolic traits mediating the effect of education on osteoarthritis risk: A Mendelian randomization study. Osteoarthritis Cartilage 2021; 29: 365–371.

24.

Wang

Bai

, et al. Osteoarthritis and the risk of cardiovascular disease: A meta-analysis of observational studies. Sci Rep 2016; 6: 1–7.

25.

Carter

Gill

Davies

, et al. Understanding the consequences of education inequality on cardiovascular disease: Mendelian randomisation study. Br Med J 2019; 365: l1855.

26.

Tillmann

Vaucher

Okbay

, et al. Education and coronary heart disease: Mendelian randomisation study. Br Med J 2017: 358: j3542.

27.

Kim

Roesler

von dem Knesebeck

. Causation or selection – examining the relation between education and overweight/obesity in prospective observational studies: A meta-analysis. Obesity Rev: Off J Int Assoc Study Obesity 2017; 18: 660–672.

28.

de Walque

. Does education affect smoking behaviors? Evidence using the Vietnam draft as an instrument for college education. J Health Econ 2007; 26: 877–895.

29.

Taylor

Richmond

Palviainen

, et al. The effect of body mass index on smoking behaviour and nicotine metabolism: A Mendelian randomization study. Hum Mol Genet 2019; 28: 1322–1330.

30.

Shpitser

. Counterfactual graphical models for longitudinal mediation analysis with unobserved confounding. Cogn Sci 2013; 37: 1011–1035.

31.

Taguri

Featherstone

Cheng

. Causal mediation analysis with multiple causally non-ordered mediators. Stat Methods Med Res 2018; 27: 3–19.

32.

Burgess

Dudbridge

Thompson

. Combining information on multiple instrumental variables in Mendelian randomization: Comparison of allele score and summarized data methods. Stat Med 2016; 35: 1880–1906.

33.

Bowden

Davey Smith

Burgess

. Mendelian randomization with invalid instruments: Effect estimation and bias detection through Egger regression. Int J Epidemiol 2015; 44: 512–525.

34.

Angrist

Krueger

. Split-sample instrumental variables estimates of the return to schooling. J Bus Econ Stat 1995; 13: 225–235.

35.

Martens

Pestman

de Boer

, et al. Instrumental variables: Application and limitations. Epidemiology 2006; 17: 260–267.

36.

Liu

Yin

Ziegler

, et al. Quantum walks: The mean first detected transition time. Phys Rev Res 2020; 2: 033113.

37.

Staiger

Stock

. Instrumental variables regression with weak instruments (No. t0151). National Bureau of Economic Research. 1994.

38.

Swanson

VanderWeele

. E-values for Mendelian randomization. Epidemiology 2020; 31: e23–e24.

39.

Lee

Wedow

Okbay

, et al. Gene discovery and polygenic prediction from a genome-wide association study of educational attainment in 1.1 million individuals. Nat Genet 2018; 50: 1112–1121.

40.

Hoffmann

Choquet

Yin

, et al. A large multiethnic genome-wide association study of adult body mass index identifies novel loci. Genetics 2018; 210: 499–515.

41.

Tachmazidou

Hatzikotoulas

Southam

, et al. Identification of new therapeutic targets for osteoarthritis through genome-wide analyses of UK Biobank data. Nat Genet 2019; 51: 230–236.

42.

Rees

Foley

Burgess

. Factorial Mendelian randomization: Using genetic variants to assess interactions. Int J Epidemiol 2020; 49: 1147–1158.

43.

Verbanck

Chen

Neale

, et al. Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases. Nat Genet 2018; 50: 693–698.

44.

Dubé

Liu

Driban

, et al. The relationship between smoking and knee osteoarthritis in the osteoarthritis initiative. Osteoarthritis Cartilage 2016; 24: 465–472.

45.

Felson

Zhang

. Smoking and osteoarthritis: A review of the evidence and its implications. Osteoarthritis Cartilage 2015; 23: 331–333.

46.

VanderWeele

Vansteelandt

Robins

. Effect decomposition in the presence of an exposure-induced mediator-outcome confounder. Epidemiology 2014; 25: 300–306.

47.

Mittinty

Lynch

Forbes

, et al. Effect decomposition through multiple causally nonordered mediators in the presence of exposure–induced mediator–outcome confounding. Stat Med 2019; 38: 5085–5102.

48.

VanderWeele

Vansteelandt

. Conceptual issues concerning mediation, interventions and composition. Stat Interface 2009; 2: 457–468.

49.

Valeri

VanderWeele

. Mediation analysis allowing for exposure–mediator interactions and causal interpretation: Theoretical assumptions and implementation with SAS and SPSS macros. Psychol Methods 2013; 18: 137–150.

50.

VanderWeele

. A unification of mediation and interaction: A four-way decomposition. Epidemiology 2014; 25: 749–761.

51.

Loh

Moerkerke

Loeys

, et al. Non-linear mediation analysis with high-dimensional mediators whose causal structure is unknown. arXiv preprint arXiv:2001.07147. 2020.

52.

Tchetgen

EJT

VanderWeele,

. On identification of natural direct effects when a confounder of the mediator is directly affected by exposure. Epidemiology 2014; 25: 282–291.

53.

Vansteelandt

Linder

Vandenberghe

, et al. Mediation analysis of time-to-event endpoints accounting for repeatedly measured mediators subject to time-varying confounding. Stat Med 2019; 38: 4828–4840.

54.

Zheng

van der Laan

. Longitudinal mediation analysis with time-varying mediators and exposures, with application to survival outcomes. J Causal Inference 2017; 5: 20160006.

55.

Hou

, et al. Exploring the causal pathway from bilirubin to cvd and diabetes in the uk biobank cohort study: Observational findings and Mendelian randomization studies. Atherosclerosis 2020; 320: 112-121.

56.

Hou

, et al. Exploring the causal pathway from ischemic stroke to atrial fibrillation: A network Mendelian randomization study. Mol Med 2020; 26: 7.

57.

Labrecque

Swanson

. Commentary: Mendelian randomization with multiple exposures: The importance of thinking about time. Int J Epidemiol 2020; 49: 1158–1162.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

1.32 MB

1.38 MB