 Research
 Open Access
 Published:
Estimating the natural indirect effect and the mediation proportion via the product method
BMC Medical Research Methodology volume 21, Article number: 253 (2021)
Abstract
Background
The natural indirect effect (NIE) and mediation proportion (MP) are two measures of primary interest in mediation analysis. The standard approach for mediation analysis is through the product method, which involves a model for the outcome conditional on the mediator and exposure and another model describing the exposure–mediator relationship. The purpose of this article is to comprehensively develop and investigate the finitesample performance of NIE and MP estimators via the product method.
Methods
With four common data types with a continuous/binary outcome and a continuous/binary mediator, we propose closedform interval estimators for NIE and MP via the theory of multivariate delta method, and evaluate its empirical performance relative to the bootstrap approach. In addition, we have observed that the rare outcome assumption is frequently invoked to approximate the NIE and MP with a binary outcome, although this approximation may lead to nonnegligible bias when the outcome is common. We therefore introduce the exact expressions for NIE and MP with a binary outcome without the rare outcome assumption and compare its performance with the approximate estimators.
Results
Simulation studies suggest that the proposed interval estimator provides satisfactory coverage when the sample size ≥500 for the scenarios with a continuous outcome and sample size ≥20,000 and number of cases ≥500 for the scenarios with a binary outcome. In the binary outcome scenarios, the approximate estimators based on the rare outcome assumption worked well when outcome prevalence less than 5% but could lead to substantial bias when the outcome is common; in contrast, the exact estimators always perform well under all outcome prevalences considered.
Conclusions
Under samples sizes commonly encountered in epidemiology and public health research, the proposed interval estimator is valid for constructing confidence interval. For a binary outcome, the exact estimator without the rare outcome assumption is more robust and stable to estimate NIE and MP. An R package mediateP is developed to implement the methods for point and variance estimation discussed in this paper.
Background
Biomedical and epidemiological studies evaluating the impact of exposure on health outcomes have received considerable attention in recent years. In addition to estimating the total effect of the exposure on the disease outcome, it is also of interest to explore potential pathways and mechanisms underlying these exposuredisease relationships. An important tool for exploring these pathways is mediation analysis [1, 2], which examining the extent to which the third variable (usually termed as mediator) mediates the observed relationship between the exposure and the outcome. In practice, the identification of a mediator can facilitate the translation of results from analytic epidemiology to improve prevention and treatment. For example, if inflammation were to be confirmed as an important mediator of the obesity–lethal prostate cancer relationship, new approaches to the prevention, prediction, and the treatment of prostate cancer would likely emerge [3]. As another example, the antiretroviral intervention referred to as “Early Access to ART for All” (EAAA) was found effective to improve the retention in care compared to the standard of care among HIVpositive patients in the MaxART study [4]. If we can identify an important mediator in the pathway between the intervention and HIV care retention, new strategies targeting the mediator may be developed to further improve outcomes.
The main goal of mediation analysis is decomposition of the total effect (TE) of the exposure into two components, a natural indirect effect (NIE) through the mediator and a natural direct effect (NDE) whose impact derives solely from the exposure. Formal definitions of NIE and NDE were first given in [5, 6] under a causal framework. Meanwhile, researchers sometimes calculate the ratio of the NIE and the TE to capture the relative importance of the mediator in explaining the pathway through which the exposure has an effect on the outcome [7]. This ratio is called mediation proportion (MP) or equivalently, the proportional mediated. There has been an increasing number of epidemiological studies that report the MP to explain the exposuredisease mechanism when conducting mediation analysis [8–13], and some of these studies include only a small to moderate sample size not exceeding 1000. Given the increasing use of MP, a comprehensive evaluation of the finitesample operating characteristics of the point and interval estimators for MP is needed to better inform practice.
Many statistical methods have been developed to identify and estimate the natural indirect effect from observational data, including but not limited to the nonparametric approach [14, 15], the weightingbased semiparametric approach [16–18], as well as the parametric outcome regression approach [2, 19–21]. This paper will focus on the parametric outcome regression approach. The parametric outcome regression approach includes two major variants, the difference method [21, 22] and the product method [2]. The difference method is a more traditional approach for mediation analysis in public health studies. It considers two regression models for the outcome with and without adjusting for the mediator. However, because those two regression models include the same response variable, they are not necessarily compatible. This model compatibility issue may restrict the application of the difference method, especially when the outcome is noncontinuous and common [21]. The product method, on the other hand, considers one regression model for the outcome and another regression model for the mediator, circumventing the model compatibility issue. Because of this potential advantage, the product method has been increasingly adopted for mediation analysis in epidemiology and public health [2, 7, 20, 23]. However, the empirical performance of the product method has not been extensively evaluated under realistic epidemiologic settings with different types of mediators and outcomes. In this article, we conduct an extensive simulation study to address this knowledge gap and to better inform practical application of the product method.
While the estimation and inference for the mediation measures via the difference method have been extensively studied [1, 19, 21, 22], several important issues remain less clear when using the product method for mediation analysis. First, the bootstrap approach has been suggested by a number of authors to construct valid confidence intervals of the mediation measures [24–27]. The explicit expressions of the closedform asymptotic variance and interval estimators and their empirical performance have not previously been sufficiently detailed. Second, when the outcome is binary, traditional mediation analyses have make the rare disease or outcome assumption and estimate the approximate NIE and MP under that assumption [2, 28]. Several recent publications have proposed exact mediation estimators for a common binary outcome [29–31], but the empirical performance of these exact estimators has not been extensively evaluated. Specifically, the relative performance of these exact estimators compared to the approximate estimators under the rare and common outcome scenarios has not been studied in detail, and not at all with a binary outcome and a continuous mediator. Third, while the empirical performance of the NIE estimator has been investigated in prior studies [19, 26, 32], there has been relatively little empirical evaluations for the MP estimates via the product method. Given that the MP is of primary interest in epidemiology and medicine, a comprehensive empirical evaluation of the MP estimator based on the product method is valuable to inform practice.
To address the aforementioned issues, we conducted an extensive Monte Carlo study to evaluate the point and interval estimators for NIE and MP under the four common data types: Case #1, a continuous outcome and a continuous mediator; Case #2, a continuous outcome and a binary mediator; Case #3, a binary outcome and a continuous mediator; as well as Case #4, a binary outcome and a binary mediator. We reviewed the counterfactual outcome framework for estimating the NIE and MP via the product method, and importantly, derived the closedform variance and interval estimators through the method of generalized estimating equations and the multivariate delta method. We provided a comparison between the closedform interval estimator and the bootstrap interval estimators via simulations with different sample sizes. In the scenario with a binary outcome and a continuous mediator, we also developed the exact NIE and MP expressions, obviating the rare outcome assumption, and evaluate the performance of the exact versus approximate expressions under rare and common outcome assumptions. Thus, an important goal of our study is to unify and supplement existing evidence on the empirical performance of the product method, especially when closedform variance is considered and when the rare outcome assumption fails to hold. To better elucidate our contribution to the literature, Table 1 summarizes the scenarios considered in our study as compared to several prior studies [19, 26, 29, 30, 32–36], across the four common data types.
The remainder of this paper is organized as follows. First, we describe the product method to obtain the point and interval estimators of NIE and MP in the aforementioned four cases with different data types. Then, we describe our Monte Carlo simulation study that investigates the performance of the point and interval estimators for NIE and MP, and report our findings. To illustrate the product method, we study how much the effect of an EAAA intervention on 12month retention is mediated by the 6month visit adherence in the MaxART study [4]. Finally, we offer a brief discussion.
Methods
Mediation measures and the product method
Assume that we have an outcome of interest, Y, an exposure, X, and a mediator, M, where each variable can be continuous or binary. We also observe W, a vector of covariates, associated with outcome and measured before the exposure, some of which may be confounders of the estimated exposureoutcome association and/or the mediatoroutcome association. A directed acyclic graph illustrating the causal relationship between those variables is shown in Fig. 1. We are interested in identifying two causal effects, the NIE and NDE. To identify these causal effects, we will follow the counterfactual framework used in the classic causal inference literature [5, 6]. Specifically, we will follow the notation in [21] and define Y_{x} and M_{x} as the outcome Y and mediator M, respectively, that would have been observed when setting X=x. Similarly, let Y_{x,m} be the value of outcome Y that would have been observed when setting X=x and M=m.
Based on the classical counterfactual framework, we made several standard assumptions. The first assumption is consistency, which assumes that (i) the potential outcome Y_{x} and potential mediator M_{x} equal their respective observed variables Y and M when we set X=x, and (ii) the potential outcome Y_{x,m} equals to the observed Y if X=x and M=m are observed [37]. Intuitively, item (i) states such that an individual’s potential outcome (and potential mediator) under his observed exposure assignment is precisely equal to their observed outcome (and observed mediator); and item (ii) states that an individual’s potential outcome given his observed history of exposure and mediator is precisely equal to his observed outcome. The consistency assumption requires that the exposure and mediator are defined and measured unambiguously. For example, the consistency assumption is plausible in randomized experiments studying welldefined medical interventions with precisely measured mediator observations. The consistency assumption may be more questionable in some observational studies where it difficult to conceptualize how the exposure of interest could be manipulated [38]. The second assumption is the composition assumption, which requires that \(Y_{x,M_{x}}\), i.e., the potential outcome Y when X=x is equal to the potential outcome Y when X=x is observed and M is set to its value corresponding to when X=x [37]. In order to identify NIE and NDE, we also require that (A.1) Y_{x,m}⊥XW, (A.2) Y_{x,m}⊥MW,X, (A.3) M_{x}⊥XW, and (A.4) \(\phantom {\dot {i}\!}Y_{x,m} \perp M_{x^{*}}  \boldsymbol W\) for all x, x^{∗}, and m. That is, (A.1)(A.3) assume that the exposureoutcome, mediatoroutcome, and exposuremediator relationships are not confounded conditional on covariates W. (A.4) is sometimes termed crossworld independence, which stresses that none of the confounders in the mediatoroutcome relationship can be affected by exposure X [31]. Intuitively, (A.1)(A.4) requires the investigators to measure all preexposure covariates that may confound the exposuremediatoroutcome mechanisms.
Under this framework, when changing the exposure levels from x^{∗} to x conditional on W=w, Nevo et al. [21] defined the NIE and NDE on a gfunction scale:
where g(.) is a prespecified monotone function. In this article, g(.) is set to the link function of the model for conditional mean Y, which will be discussed later. Given (1), the TE is defined as summation of NIE and NDE, that is,
Finally, the MP is given by the ratio of NIE and TE.
The above mediation measures are derived while fixing W to w while changing the exposure level from x^{∗} to x, and therefore are conditional causal parameters. Alternatively, definitions of mediation measures have been proposed by averaging over the distribution of W. As discussed in [14, 15], under assumptions (A.1)–(A.4), the NIE, NDE, and TE are given by \(\text {NIE}(x^{*},x) = E[Y_{x,M_{x}}]  E[Y_{x,M_{x^{*}}}]\), \(\text {NDE}(x^{*},x) = E[Y_{x,M_{x^{*}}}]  E[Y_{x^{*},M_{x^{*}}}]\), and \(\text {TE}(x^{*},x) = E[Y_{x,M_{x}}]  E[Y_{x^{*},M_{x^{*}}}]\), respectively. Estimation and inference for the marginal mediation measures has been discussed in [15, 16, 39], based on nonparametric or semiparametric models for the outcome and mediator, as commonly appear in practice, and will not be pursued here. The focus of this paper is estimation and inference mediation measures conditional on W=w, as shown in (1).
We now introduce the product method for estimating those mediation effects. Specifically, we assume the following conditional mean model for the outcome (Y),
where g(.) is a link function, β_{1} is the exposure effect on the outcome conditional on the effects of the mediator and confounders, β_{2} represents the relationship between the mediator and outcome conditional on the effect of the exposure and confounders. Common link functions include the identity function when the outcome is continuous and logistic function when the outcome is binary. Because previous empirical evidence suggests that there are few interaction effects between an intervention/exposure and covariates that replicated across studies in public health and epidemiology [40], we assume there are no mediator–exposure interactions in (2). Mediation analyses in the presence of mediator–exposure interaction effects are studied elsewhere, for example, [23, 28, 29]. Researchers can empirically verify this assumption in their data before applying these methods.
In addition to the outcome model (2), the product method additionally requires fitting the following model for the mediator:
where γ_{1} represents the association between the exposure to the mediator conditional on the effects of the covariates, h is a link function, which can be a identity function and a logistic function when the mediator is continuous and binary, respectively. For simplicity of notation but without loss of generality, we assume that the mediator model and the outcome model share the same set of covariates. We can set some elements in β_{3} or γ_{2} to zero when the covariate sets in the outcome and mediator models are not exactly the same.
Next, we provide expressions of the mediation measures under the scenarios of continuous and binary mediator, separately.
Continuous mediator
When the mediator is continuous and h is an identity link function, model (3) becomes
First consider Case #1, where both the outcome and mediator are continuous and g is an identity link function. Case #1 is the most basic scenario in epidemiologic and public health studies. For example, Pandey and Shrivastava (2017) [41] studied the mediation effect of social support on the association between hardiness and immune response, where the immune response was treated as a continuous outcome and measured by CD4+ Tlymphocyte counts; the level of social support was considered as a continuous mediator measured by a social support score developed by Pandey and Shrivastava (2016) [42]. As indicated by [28], if the identification assumptions hold and the outcome as well as the mediator models are correctly specified, the NIE, NDE and TE can be expressed as β_{2}γ_{1}(x−x^{∗}), β_{1}(x−x^{∗}) and (β_{2}γ_{1}+β_{1})(x−x^{∗}), respectively. The mediation proportion is given as \(\frac {\beta _{2}\gamma _{1}}{\beta _{2}\gamma _{1}+\beta _{1}}\). Of note, Case #1 only requires that the conditional mean in (4) is correctly specified and does not require further distributional assumptions for the mediator; e.g., even if MX,W is not normally distributed or Var(MX,W) is heteroskedastic, the above expressions are valid.
Scenarios with a binary outcome and a continuous mediator (Case #3) is also a common scenario in public health research. For instance, VanderWeele et al. (2012) [43] examined the extent to which the effect of two genetic variants on lung cancer is mediated by cigarette smoking, where the presence of lung cancer is considered as a binary outcome and the amount of smoking (cigarettes/day) is considered as a continuous mediator. In this data type, g(.) is usually set to logistic link function, and the NIE, NDE, and TE can be defined on the log odds ratio scale. In Web Appendix A (see additional file 1), we derive exact expressions for the mediation measures under the conditions that MX,W in model (4) follows a normal distribution with a constant variance σ^{2}. Specifically, define an integrand
we have that
both of which involve onedimensional logisticnormal integrals that do not have closedform solutions. Gaynor et al. (2019) [29] uses a probit function to approximate the logistic function in the integral and obtained closedform expressions for the mediation measures (Table 2). However, the probit approximation tends to be inaccurate as the outcome prevalence deviates from 50%, as discussed in [29] and Web Appendix A (see additional file 1). Web Figure 1 visualizes the percent bias of the probit approximation approach. It is show that the probit approximation method exhibits high bias for calculating NIE if the baseline prevalence is less than 10% or greater than 80%. Instead of using probit function to approximate the logisticnormal integrals, we consider in this paper the GaussHermite Quadrature (GHQ) approach [44] to numerically calculate the integrals.
When analyzing data with a dichotomous health outcome, a frequent scenario is that the health outcome is rare in that the probability of the outcome is small, e.g. less than 10%. In epidemiologic cohort studies of cancer and heart disease, the outcome is extremely rare. Although there is no clear cutoff to consider the outcome as rare, some authors have considered the rare outcome assumption to be plausible when the outcome rate is less than 2% (e.g. [45]). In mediation analysis, the mediation formulas can be simplified when the outcome is rare. Specifically, we can approximate the NIE, NDE, and TE by β_{2}γ_{1}(x−x^{∗}), β_{1}(x−x^{∗}), and (β_{2}γ_{1}+β_{1})(x−x^{∗}), respectively, under the rare outcome assumption[23, 28]. It follows that \(\text {MP}\approx \frac {\beta _{2}\gamma _{1}}{\beta _{2}\gamma _{1}+\beta _{1}}\). Here we can obtain simple and closedform expressions for the mediation measures because the logistic link function in model (2) approximates a log link function, and the log binomial model is collapsible. Therefore, the logisticnormal integrals in the mediation expressions can be approximated by the lognormal integrals that have closedform solutions [23]. In Web Appendix A (see additional file 1), we provide more details on the validity of the approximate expressions when the outcome is rare. In order to distinguish the approximate expressions under a rare outcome assumption from the exact expressions, we write those approximate mediation measures as NIE^{(a)}, NDE^{(a)}, TE^{(a)}, and MP^{(a)}, respectively.
Binary mediator
In this section, we describe the product method estimators when M is a binary variable, based on the following logistic regression model for M:
First, consider the case when outcome is continuous (Case #2) and g(.) is the identity link function. As an example of this data type, Li et al. (2007) [46] evaluated the mediation effect of chronic cerebral infarction in the causal relationship between the apolipoprotein E ε4 allele and the cognitive function, in which the cognitive function is a continuous measure evaluated by 19 neuropsychological performance tests and the chronic cerebral infarction serves as a binary mediator. Throughout this Section, we define the exponentiated linear component \(\kappa (x,\boldsymbol w)=\exp (\gamma _{0}+\gamma _{1} x + \boldsymbol \gamma _{2}^{T} \boldsymbol w)\). If the identification assumptions (A.1)–(A.4) hold and the outcome as well as mediator models are correctly specified, we have
and NDE=β_{1}(x−x^{∗}), as shown in [32]. As a result, the mediation proportion is given by \(\text {MP} = \frac {\text {NIE}}{\text {NDE}+\text {NIE}} \).
Finally, consider a binary mediator and a binary outcome (Case #4), where we fit the mediator model (7) and outcome model (2) with logistic link functions. In the illustrative example shown later, we examine how much the effect of the EAAA intervention on 12month retention is mediated by 6month visit adherence in the MaxART study. Here, the 12month retention and 6month visit adherence are considered as a binary outcome and mediator, respectively. In Case #4, as long as the identification assumptions (A.1)–(A.4) hold and (2) and (7) are correctly specified, the NIE and NDE on a log odds ratio scale is given by [30, 31]:
where \(\eta (x,\boldsymbol w)=\exp (\beta _{0}+\beta _{1} x+\beta _{3}^{T} \boldsymbol w)\). Given NIE and NDE, the MP is given by \(\frac {\text {NIE}}{\text {NIE+NDE} }\). If the outcome is rare, the following approximate NIE and NDE have been widely used [2]:
and NDE^{(a)}=β_{1}(x−x^{∗}). As a result, \(\text {MP}^{(a)}=\frac {\text {NIE}^{(a)}}{\text {NIE}^{(a)}+\text {NDE}^{(a)} }\). When the outcome is rare, we provide an explanation that the approximate expressions are valid by exploiting the similarity between the logistic and log link functions (see Web Appendix B in additional file 1).
Point and interval estimates for NIE and MP
In the previous two sections, we provided expressions for NIE and MP in Cases #1–4, which are functions of the unknown parameters in the outcome model \(\boldsymbol \beta = [\beta _{0},\beta _{1},\beta _{2},\boldsymbol \beta _{3}^{T}]^{T}\), and unknown parameters in the mediator model, \(\boldsymbol \gamma =[\gamma _{0},\gamma _{1},\boldsymbol \gamma _{2}^{T}]^{T}\). Exact mediation expressions for mediation measures in Case #3 also involve the variance of the error term in the mediator model (7), σ^{2}. Let θ denote all unknown parameters used in the expressions of NIE and MP, which is [β^{T},γ^{T},σ^{2}]^{T} for the exact expressions in Case #3 and [β^{T},γ^{T}]^{T} for other cases. Hereafter, we will rewrite the expressions of NIE and MP as \(\mathscr {NIE}(\boldsymbol \theta)\) and \(\mathscr {MP}(\boldsymbol \theta)\), respectively, to emphasize that those expressions are functions of θ.
In practice, the coefficients in the outcome model, \(\hat {\boldsymbol \beta }\), can be obtained by solving the following generalized estimating equation (GEE) [47]:
where \(\{Y_{i},X_{i},M_{i},\boldsymbol W_{i}\}_{i=1}^{n}\) are n observations of {Y,X,M,W} and V_{i} is a working variance term for the outcome Y_{i}. The optimal asymptotic efficiency of \(\hat {\boldsymbol \beta }\) will be obtained when V_{i}=Var(Y_{i}). When V_{i} is misspecified, the resultant \(\hat {\boldsymbol \beta }\) is still consistent under mild regularity conditions but could be less efficient [48]. Similarly, the coefficients in the mediator model, \(\hat {\boldsymbol \gamma }\), can be obtained by solving the GEE below
where \(V_{i}^{*}\) is the working variance for M_{i}. Also, misspecification of \(V_{i}^{*}\) only impacts the efficiency of \(\hat {\boldsymbol \gamma }\). When evaluating the exact expressions in Case #3, \(\hat \sigma ^{2}\) is needed and can be estimated by solving the GEE
After obtaining \(\hat {\boldsymbol \theta }\) through the above estimating equations, we can calculate \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) by plugging those parameter estimates in their expressions introduced in the previous sections; i.e., \(\widehat {\text {NIE}}=\mathscr {NIE}(\hat {\boldsymbol \theta })\) and \(\widehat {\text {MP}}=\mathscr {MP}(\hat {\boldsymbol \theta })\). In what follows, we develop the closedform asymptotic variance expressions of \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) based on the multivariate delta method, and review the nonparametric bootstrap approach to obtain the confidence interval estimators.
The multivariate delta method first estimates the variancecovariance matrix of \(\hat {\boldsymbol \theta }\), abbreviated by \(\hat {\boldsymbol \Sigma }_{\boldsymbol \theta }\), based on the theory of estimating equations [47]. Specifically, \(\text {Var}(\hat {\boldsymbol \beta })\), \(\text {Var}(\hat {\boldsymbol \gamma })\), and \(\text {Var}(\hat {\sigma }^{2})\) can be estimated by the robust sandwich variance estimators [47] based on their respective estimating equations. If the working variance terms, V_{i} and \(V_{i}^{*}\), are correctly specified when estimating \(\hat {\boldsymbol \beta }\) and \(\hat {\boldsymbol \gamma }\), one can also approximate \(\text {Var}(\hat {\boldsymbol {\beta }})\) and \(\text {Var}(\hat {\boldsymbol {\gamma }})\) by the negative inverse information matrix, \(\left \{\frac {\partial U(\boldsymbol {\hat {\boldsymbol {\beta }}})}{\partial {\boldsymbol \beta }^{T}}\right \}^{1}\) and \(\left \{\frac {\partial U(\hat {\boldsymbol \gamma })}{\partial \boldsymbol {\gamma }^{T}}\right \}^{1}\), respectively. In general, to obtain \(\hat {\boldsymbol \Sigma }_{\boldsymbol \theta }\), we also need to estimate the covariances between \(\hat {\boldsymbol {\beta }}\), \(\hat {\boldsymbol {\gamma }}\), and \(\hat {\sigma }^{2}\). However, as we state below, these asymptotic covariances are 0 since the estimating equations are asymptotically uncorrelated. We formalize this result below.
Result 1.
Estimators \(\hat {\boldsymbol \beta }\), \(\hat {\boldsymbol \gamma }\), and \(\hat \sigma ^{2}\) obtained by solving U(β)=0, U(γ)=0, and U(σ^{2})=0 are asymptotically uncorrelated and have zero asymptotic covariances.
Detailed proof of Result 1 is provided in Web Appendix C (see additional file 1). Following this result, we have that \(\hat {\boldsymbol \Sigma }_{\boldsymbol \theta } = \left [\begin {array}{ccc} \widehat {\text {Var}}(\hat {\boldsymbol \beta }) & \boldsymbol 0 & \boldsymbol 0\\ \boldsymbol 0 & \widehat {\text {Var}}(\hat {\boldsymbol \gamma }) & \boldsymbol 0\\ \boldsymbol 0 & \boldsymbol 0 & \widehat {\text {Var}}(\hat {\sigma }^{2}) \\ \end {array}\right ]\)for the exact expressions in Case #3 and \(\hat {\boldsymbol \Sigma }_{\boldsymbol \theta } = \left [\begin {array}{cc} \widehat {\text {Var}}(\hat {\boldsymbol \beta }) & \boldsymbol 0\\ \boldsymbol 0 & \widehat {\text {Var}}(\hat {\boldsymbol \gamma }) \\ \end {array}\right ]\) for the remaining Cases. Finally, the variances of \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) is obtained through the multivariate delta method [49]:
Given these variance estimators, the 95% confidence intervals of NIE and MP can be computed by normal approximation as \(\Big \{\widehat {\text {NIE}}  1.96 \times \sqrt {\widehat {\text {Var}}(\widehat {\text {NIE}})}, \widehat {\text {NIE}} + 1.96 \times \sqrt {\widehat {\text {Var}}(\widehat {\text {NIE}})}\Big \}\) and \(\Big \{\widehat {\text {MP}}  1.96 \times \sqrt {\widehat {\text {Var}}(\widehat {\text {MP}})}, \widehat {\text {MP}} + 1.96 \times \sqrt {\widehat {\text {Var}}(\widehat {\text {MP}})}\Big \}\), respectively.
Alternatively, the nonparametric percentile bootstrap approach [50, 51] can be used to approximate the empirical distributions of \(\widehat {\text {NIE}}\) by resampling the dataset with replacement and reestimating all model parameters. The values for NIE are then calculated for each bootstrap dataset. This step of resampling and calculating NIE is repeated for a large of times (e.g., 1,000 times), and then the bootstrap distribution of NIE is obtained from the collection of estimates based on resampled datasets. Finally, the percentile bootstrap approach employs the 2.5% and 97.5% percentiles of the bootstrap sample distribution to obtain a 95% confidence interval of NIE. Similarly, the percentile bootstrap approach can be used to obtain a 95% confidence interval of MP. When implementing the product method, several studies [24–27] suggested using the bootstrap approach to calculate confidence intervals for the mediation measures, because the finitesample distribution of the mediation measure estimators, especially the MP estimator, can be skewed and may not be well approximated by the normal distribution. In what follows, we will compare the empirical performance between the closedform sandwich confidence interval estimator and the bootstrap interval estimator, and assess when each approach provides accurate empirical coverage for all 4 cases.
Results
Simulation design
We conducted simulation studies under a range of scenarios likely to be encountered in practice to assess the performance of the point and interval estimators of NIE and MP. To represent a wide range of sample sizes commonly seen in practice, we considered 4 levels of sample sizes: 150, 500, 1,000, and 5,000 for the continuous outcome cases (Cases #1 and #2). With a baseline outcome prevalence of 3%, in order to obtain 15, 30, 150, 600 disease cases in expectation, we considered sample sizes of 500, 1,000, 5,000, 20,000 for the binary outcome scenarios (Cases #3 and #4). Here, we did not consider a sample size at 150 for the binary outcome scenarios since the expect number of cases (about 5 cases in each simulated dataset) would be too small to obtain stable NIE and MP estimates. For each data type and sample size, we considered the exposure (X) as a binary variable and for simplicity assumed that there were no confounders in both the mediator and outcome models. When the outcome was continuous, we set TE ∈(0.25,0.5,1), indicating small, medium and large total effects on the difference scale. When the outcome was binary, we set the TE measured on the log odds ratio scale to be log(1.2), log(1.5), and log(2). For each value of TE, we considered MP ∈(0.05,0.2,0.5). All the mediation measures are defined for a change in X from 0 to 1. With a smaller TE, it is expected that a larger sample size is needed to obtain a reasonably precise mediation estimator because the information available or estimating even the overall exposure effect is rather limited. Furthermore, for a fixed TE, a smaller MP or equivalently a smaller NIE would require a larger sample size for reasonably precise mediation analysis since the mediation effect size becomes limited.
Case #1, continuous outcome and continuous mediator. We first generated the exposure X∼Bernoulli(0.5). Then, we simulated the mediator M∼N(γ_{0}+γ_{1}X,1), where γ_{0}=0. γ_{1} was chosen to 0.408, corresponding to the exposuremediator correlation, Corr(X,M)=0.2. Finally, given X and M, we simulated Y∼N(β_{0}+β_{1}X+β_{2}M,1), where β_{0} was fixed as 0. We let β_{1}=(1−MP)×TE and \(\beta _{2}=\frac {\text {MP}\times \text {TE}}{\gamma _{1}}\) based on the relationships that NDE=(1−MP)×TE=β_{1} and NIE=MP×TE=β_{2}γ_{1}.
Case #2, continuous outcome and binary mediator. First, we simulated X∼Bernoulli(0.5). Then, we generated M conditional on X by the logistic regression logit(P(M=1X))=γ_{0}+γ_{1}X. Noting that \(\gamma _{0}=\log \left (\frac {P(M=1X=0)}{1P(M=1X=0)}\right)\), we chose the values of γ_{0} such that the baseline prevalence of the mediator P(M=1X=0)=0.2. Similar to Case #1, we chose γ_{1}=0.903 to give an exposuremediator correlation 0.2. Finally, we simulated Y∼N(β_{0}+β_{1}X+β_{2}M,1), where β_{0}=0, β_{1}=(1−MP)×TE from the definition NDE=(1−MP)×TE=β_{1}, and \(\beta _{2} = {\text {MP}\times \text {TE}}/\left \{\frac {e^{\gamma _{0}+\gamma _{1} }}{ 1+ e^{\gamma _{0}+\gamma _{1} }}  \frac {e^{\gamma _{0} }}{ 1+ e^{\gamma _{0}}}\right \}\) by noticing \(\text {MP}\times \text {TE} = \text {NIE} = \beta _{2} \left \{\frac {e^{\gamma _{0}+\gamma _{1} }}{ 1+ e^{\gamma _{0}+\gamma _{1} }}  \frac {e^{\gamma _{0} }}{ 1+ e^{\gamma _{0}}}\right \}\).
Case #3, binary outcome and continuous mediator. We first generated X and M using the same procedure in Case #1. Now, given X and M, we used the following logistic regression model to simulate Y,
Since \(\beta _{0}=\log \left (\frac {P(Y=1X=0,M=0)}{1P(Y=1X=0,M=0)}\right)\), we selected β_{0} so the baseline outcome prevalence was 3%. Then, we selected β_{1} and β_{2} by numerically solving the system of equations \(\text {MP}\times \text {TE} = \mathscr {NIE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) and \((1\text {MP})\times \text {TE} = \mathscr {NDE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\), where \(\mathscr {NIE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) and \(\mathscr {NDE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) refer to the exact expressions of NIE and NDE, given in (5) and (6).
Case #4, binary outcome and binary mediator. We first generated X and M using the same procedure in Case #2, then generated the outcome Y using the logistic regression model (10). The values of β_{0}, β_{1} and β_{2} were obtained as follows. We chose β_{0} such that the baseline outcome prevalence is around 3%. We then chose β_{1} and β_{2} by solving the system of equations \(\text {MP}\times \text {TE} = \mathscr {NIE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) and \((1\text {MP})\times \text {TE} = \mathscr {NDE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\), where \(\mathscr {NIE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) and \(\mathscr {NDE}(\beta _{0},\beta _{1},\beta _{2},\gamma _{0},\gamma _{1})\) refer to the exact expressions of NIE and NDE, given in (8) and (9).
For each data type, we obtained \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\), the variance estimates of \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) by the multivariate delta method, and 95% confidence intervals of NIE and MP by the multivariate delta method and percentile bootstrap approach. We did not use the bootstrap approach for the sample size of 20,000 in Cases #3 and #4, due to the computational cost of performing the number of bootstrapping samples with a large dataset across 5,000 replications. With a binary outcome, we first evaluated the performance of the estimators based on the approximate expressions (i.e., NIE^{(a)} and MP^{(a)}), and compared the performance between the approximate and the exact estimators for the mediation measures by varying the baseline outcome prevalence from 1% to 50%.
The percent bias (Bias(%)) was used to evaluate the accuracy of the point estimates of NIE and MP. It was calculated as the median bias relative to the true value over 5,000 replications. We employed the “median” rather than the “average” value over all the replications to avoid the undue influence of outliers on the results when the sample sizes were not large. The variance ratio was defined as the ratio between the median of the estimated variance and the empirical variance, and is used to determined the accuracy of the variance estimator obtained by the multivariate delta method. The accuracy of the interval estimator is determined by calculating the empirical coverage rate (CR) of 95% confidence interval across 5,000 replications. The simulation results for the point, variance and interval estimates of NIE and MP under the four data types are shown in Table 3 (Case #1), Web Table 1 (Case #2), Web Table 2 in additional file 1 (Case #3) and Web Table 3 in additional file 1 (Case #4), respectively. Detailed findings from the simulation study are reported below.
Estimation of NIE
When the outcome is continuous, the point estimates of NIE generally had minimal bias for sample sizes equal or greater than 500, for all values of TE and MP considered. When the sample was as small as 150, the NIE point estimates had relatively small bias (percent bias within ±10%) as long as MP≥0.2 and TE≥0.5. With binary outcome, however, the point estimates did not show sufficiently small percent bias until the sample size was 1,000, as shown in Web Table 2 in additional file 1 (Case #3) and Web Table 3 in additional file 1 (Case #4). When the outcome was relatively rare and sample size was not large, there was not enough data to accurately estimate the NIE by the product method. In Case #3 (Web Table 2 in additional file 1), we also observed that the percent bias of \(\widehat {\text {NIE}}^{(a)}\) was notable when TE= log(2) and MP=0.5 even when the sample size was 20,000, indicating that with a rare outcome, bias persisted when TE and MP were also large.
The variance ratio of the NIE under all data types and sample sizes were very close to 1, indicating that the variance estimator derived by the multivariate delta method was accurate. Compared to the multivariate delta method, the bootstrap provided more accurate NIE interval estimates especially when the sample size is small. When the outcome was continuous, the bootstrap coverage rates were close to or greater than 95% even when the sample size was 150. In general, both the bootstrap and the multivariate delta method had confidence interval coverage rates very close to the nominal value when the sample size ≥500, except for the scenario when both TE and MP were large in Case #3, in which case both methods exhibited lower coverage rates than nominal as sample size increased. This is because the rare outcome assumption may be inadequate and so the resulting bias of \(\widehat {\text {NIE}}^{(a)}\) are more pronounced with a larger sample size. Overall, these findings suggested that the multivariate delta method is a valid approach for estimating the variance and confidence intervals for the NIE when the sample size is at least 500 with a continuous outcome.
For the cases of binary outcome, we also evaluated the performance of the NIE estimates based on the exact expressions (i.e., \(\widehat {\text {NIE}}\)). The results are shown in Web Tables 2 and 3 in (see additional file 1) or Cases #3 and #4, respectively. As long as the sample sizes are greater than or equal to 1,000, the estimator based on the exact expressions consistently carried small percent bias and nominal coverage rate.
Estimation of MP
In contrast to the NIE estimates, the point estimates of MP often diverged from the true MP values with smaller sample sizes. When the outcome was binary, the percent bias in MP was usually larger than 30% for sample sizes ≤1,000, and sometimes even larger than 70% for sample size ≤500, as shown in Web Table 2 (Case #3) and Web Table 3 (Case #4) in additional file 1. The MP point estimates were more stable with a continuous outcome, in which case the bias is negligible when the sample size was at least 500. For the binary outcome scenarios, the MP percent bias was not close to 0 until the sample size was 5,000 or the number of cases ≥200. For all data types, the percent bias of \(\widehat {\text {MP}}\) appears larger when the TE is small. For the same TE, a smaller percent bias occurred when the MP was larger.
The variance estimators of \(\widehat {\text {MP}}\) obtained from the multivariate delta method were smaller than their corresponding empirical variances, and this phenomenon became more noticeable with small sample sizes and a small TE. In general, the multivariate delta method provided accurate variance estimators when the sample size was at least 500 for the continuous outcome scenarios. When the outcome was binary, the multivariate delta method only provided accurate variance estimators when the sample size was at least 5000 and TE ≥ log(1.5).
Similar to the NIE interval estimates, the bootstrap provided more accurate MP interval estimates than the multivariate delta method, especially with smaller sample sizes, because the distribution of \(\widehat {\text {MP}}\) with small sample sizes can deviate from normality, which was assumed in the multivariate delta method. For all data types and sample sizes, the bootstrap approach generally had coverage rates higher than 95% and began to provide close to nominal coverage rates for sample size of 500 for scenarios with continuous outcome and 5,000 for the scenarios with binary outcome. On the other hand, the multivariate delta method did not provide accurate confidence intervals for smaller sample sizes, and its coverage rates sometimes dropped below 80% in some settings with a binary outcome. With a continuous outcome, a sample size of 1,000 was needed for the multivariate delta method to provide satisfactory interval estimates, whereas the sample size requirement for the binary outcome scenarios was substantially larger. As shown in Web Table 2 (Case #3) and Web Table 3 (Case #4) in additional file 1, the multivariate delta method provided satisfactory coverage rates when TE≥ log(1.5) and sample size ≥5,000 with a binary outcome. With a small TE (TE = log(1.2)), a sample size of at least 20,000 will be needed for the multivariate delta method to provide MP interval estimates with nominal coverage. The performance of the delta method was also sensitive to the magnitude of TE. When the TE was small, the MP confidence interval tended to be much wider than it should be when the MP was also small, but tends to be narrower than it should be for larger MP.
The impact of outcome prevalence with a binary outcome
We conducted additional simulations to compare the mediation analyses based on the approximate expressions and exact expressions (\(\widehat {\text {NIE}}^{(a)}\) v.s. \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}^{(a)}\) v.s. \(\widehat {\text {MP}}\)) when the baseline outcome prevalence was varied from 1% to 50%. We considered TE ∈{log(1.2), log(2)}, MP ∈{0.1,0.5} and a large sample size of 20,000 to alleviate concerns on smallsample biases.
Figure 2 and Web Figure 2 in additional file 1 present the results for MP and NIE estimates, respectively, for the binary outcome and binary mediator scenario (i.e., Case #4). The estimates based on the exact causal mediation expressions, \(\widehat {\text {MP}}\) and \(\widehat {\text {NIE}}\), provided accurate point and interval estimates when the outcome prevalence >1%. When TE= log(1.2) and the baseline outcome prevalence was 1%, \(\widehat {\text {MP}}\) and \(\widehat {\text {NIE}}\) showed significant negative percent bias, as the number of cases (about 200) were quite small. The NIE and MP estimates based on the approximate expressions did not generally exhibit satisfactory performance when the outcome prevalence was high. Specifically, when TE > log(1.2), the percent bias of \(\widehat {\text {MP}}^{(a)}\) diverged from 0 and the coverage rate of \(\widehat {\text {MP}}^{(a)}\) by the delta method substantially decreased from 95%. Similarly, \(\widehat {\text {NIE}}^{(a)}\) was also quite different from its true value when the baseline outcome prevalence ≥10% (See Web Figure 1 in additional file 1).
With a binary outcome and continuous mediator (Case #3), \(\widehat {\text {MP}}\) and \(\widehat {\text {NIE}}\) provided accurate point and interval estimates among all levels of baseline outcome prevalence, as shown in Web Figures 3 and 4 (see additional file 1). The estimator \(\widehat {\text {MP}}^{(a)}\) also provided very robust point and interval estimates (Web Figure 2 in additional file 1) for the common outcome scenarios, where its percent bias was less than 1% among all TEs, MPs and outcome prevalences considered. However, the performance of \(\widehat {\text {NIE}}^{(a)}\) was sensitive to the baseline prevalence and the magnitude of MP. For example, when MP=0.5, the percent bias of \(\widehat {\text {NIE}}^{(a)}\) increased as baseline outcome prevalence increased, and the confidence interval coverage rates rapidly declined. When the true MP=0.1, the percent bias of \(\widehat {\text {NIE}}^{(a)}\) was negligible over the range of outcome prevalences considered.
Homoscedasticity and normality assumptions in case #3
In the above simulations, we simulated MX under a homoscedastic normal distribution, as assumed in deriving the mediation measure expressions in Case #3. As a sensitivity analysis, here we evaluate the performance of the product method in Case #3 based on estimands (5) and (6), when M is either heteroskedastic or not normally distributed.
First, we assessed the robustness of the proposed method when the homoscedasticity assumption is violated. Specifically, given X, we generated the mediator M by \(MX \sim N(\gamma _{0} + \gamma _{1} X,\sigma _{x}^{2})\), where \(\sigma _{x}^{2}=\eta _{0}+\eta _{1} X\). Here, η_{1} measures the extent to which the variance depends upon the exposure and η_{0} was chosen such that the expectation of \(\sigma _{x}^{2}\) equals to 1 in all cases considered, thereby holding the overall variance constant to study heteroscedasticity. When η_{1}=0, homoskedasticity holds, and a higher value of η_{1} represents a stronger degree of heteroskedasticity. We set η_{1}∈{0,0.25,0.5,0.75,1} such that the variance of the mediator in the exposed group is κ∈{1,1.26,1.67,2.20,3} times that in the unexposed group. We used the proposed exact and approximate estimators to estimate NIE and MP. The results are given in Web Table 4 (see additional file 1). We found that both the exact and approximate estimators are relatively insensitive to the violations of the homoscedasticity assumption. In particular, the exact estimators exhibited minimal bias and nominal confidence interval coverage under all degrees of heteroskedasticity considered and the approximate estimator also exhibited good performance when the outcome was rare (outcome prevalence at 3%). These results suggest that the proposed estimators are robust to heteroscedasticity, at least the forms of this we considered.
Second, we examine the sensitivity of proposed estimators when the error term is not normally distributed. We follow the data generating process for Case #3, except that we simulated M=γ_{0}+γ_{1}X+ε, where ε follows a gamma distribution. Specifically, we let \(\epsilon \sim \frac {bE[b]}{\sqrt {\text {Var}(b)}}\), where b follows a gamma distribution with density \(f(b) = \frac {1}{\Gamma (k)\theta ^{k}} b^{k1} e^{b/\theta }\) (b>0), k and θ are shape and scale parameters. We subtract b with its expectation E[b]=kθ and then divide it by its standard deviation \(\sqrt {k}\theta \) in order to fix the mean and variance of ε at 0 and 1, matching the first two moments of the standard normal distribution. Here, we chose k=(2/s)^{2} and θ=s/2 such that the coefficient of skewness of ε was s. We choose s to be 1, 1.5 and 2, representing different degrees of skewness, and then still use our exact NIE and MP estimator based on (5) and (6), and the corresponding approximate estimators for analysis. The simulation results are presented in Web Table 5 (see additional file 1). We also included a scenario that ε follows a standard normal distribution as a benchmark. We observed that both the exact and approximate method are robust with regard to violations of normality assumption, where all the point and variance estimators and confidence interval coverage rates are comparable across the two data generating processes. This finding suggests that our exact NIE and MP estimators based on (5) and (6) are insensitive to moderate skewness of the outcome distribution, even though the derivation assumes normality.
In the presence of a binary confounder
To assess the robustness of our simulation findings with observed confounding, we conducted additional simulation studies in the presence of a binary confounder variable that is associated with the exposure, mediator and outcome. We investigated the NIE and MP estimators in terms of percent bias, variance ratio, and 95% confidence interval coverage under each data type, by varying the confounderexposure, confoundermediator, and confounderoutcome associations. The simulation design and results are described in detail in Web Appendix D and Web Tables 6–7 in additional file 1. Overall, we observed that the performance of the product method are robust under confounding of different magnitudes and structures. The percent bias, variance ratio and 95% confidence interval coverage rate are comparable across the scenarios with and without a binary confounder.
Application to the maxART study
We performed mediation analysis in the MaxART study [4], which is a steppedwedge cluster randomized trial among HIVpositive participants in Eswatini. The primary objective of the study was to understand the impact of early access to antiretroviral therapy (EAAA) versus standard of care (SoC). From September 2014 to August 2017, the MaxART Consortium randomly assigned 14 participating clinics in pairs to shift from SoC to EAAA at randomly chosen prespecified dates. Further details of the design of this MaxART study can be found in [52]. The MaxART study previously found that EAAA improved retention in HIV care [4], but the mechanisms underlying the interventionretention relationship is unknown. In this illustrative example, we investigated the extent to which the effect of the intervention (SoC v.s. EAAA) on 12month retention in HIV care was mediated by visit adherence at 6 months. Participants were classified as retained in HIV care for 12 months if, at the end of the 12th month post enrollment, the participant was alive and had not discontinued treatment, where either the last clinic visit was less than 90 days from the end of study or next scheduled visit date was within 30 days from the end of study. In order to obtain 12month retention, we required (1) the participant’s enrollment date to be longer than 12 months from end of the study and/or (2) if initially receiving SoC treatment, the participant’s transition date to EAAA was longer than 12 months from enrollment. Participants who did not meet the above two requirements were excluded. Finally, 1,731 participants were used in our illustrative analysis, with 1,335 individuals retained in care for 12 months and 396 individuals not retained, of whom 1,014 individuals received SoC and 717 individuals received EAAA. Baseline characteristics of participants are given in Web Table 8 (see additional file 1). For purpose of illustration, we do not addressing clustering of subjects within clinics; we do note, however, minimal degree of clustering has been reported previously in the previous analysis of the MaxART study [4].
The hypothesized mediator considered here was 6month visit adherence, which measures whether a participant’s frequency of clinical visits coincides with the MaxART protocol over the first 6 months following enrollment. According to the MaxART protocol, participants are expected to have a followup visit in every 30 days. It follows that at the end of the 6th month the participants should have completed 6 or more visits. Here, the definition of 6month visit adherence completion 5 or more clinical visits by the end of the 6th month after enrollment (yes=1, no=0). Based on this definition, 831 participants adhered to the MaxART visit schedule in the first 6 months, whereas 900 participants did not.
We considered two scenarios for confounding adjustment (i.e., W) in the outcome and mediator models. In Scenario I, we only adjusted for the steptime. In Scenario II, we adjusted for all factors that may have been confounders of the interventionretention relationship, visit adherenceretention relationship, or interventionvisit adherence relationship. From clinical knowledge and prior analyses of these data, the comprehensive set of potential confounders included steptime, age at study enrollment (<20 yrs, [20,30) yrs, [30,40) yrs, [40,50) yrs, [50,60) yrs, ≥60 yrs), sex, marital status (married, divorced/widowed, single), education (illiterate/primary, secondary, high school, and tertiary), CD4 counts (<350 cells/ul, [350,500] cells/ul, >500 cells/ul), WHO stage (I, II, III and IV stages), BMI (<18.5, [18.5,25), [25,30), ≥30kg/m^{2}), screened for TB symptoms (yes, no), viral load (<5000 copies/ml, [5000,30000] copies/ml, >30000 copies/ml), treatment support (yes, no), level of clinic (hospital, clinic with maternity ward, clinic without maternity ward), time from HIV tested positive to enrollment (<1 yr, 13 yrs, >3 yrs), and clinic volume (low: < median, high ≥ median). For simplicity and consistency with the primary analysis of the MaxART study [4], the missing indicator method [53] was used to account for missing covariates, where missing data was treated as a separate group for each of the confounding variables in the models.
We first implemented the product method based on the approximate mediation measure expressions assuming a rare outcome. We coded nonretention as 1 and retention as 0. We calculated the NIE, TE and MP comparing EAAA to SoC, conditional on the mode for each model covariate. Although we estimated the mediation measures at the mode of each model covariate, these measures could also be calculated at other values of the covariates as well. Results are given in Table 4. In Scenario I, the steptimeadjusted model, we found that the intervention was protective against 12month nonretention, with odds ratios of 0.23 and 0.55 for \(\widehat {\text {TE}}^{(a)}\) and \(\widehat {\text {NIE}}^{(a)}\) respectively. Because the 95% confidence intervals, either by the delta method or bootstrap, for both parameters excluded the null, we conclude that both effects were significantly different from zero. The steptimeadjusted \(\widehat {\text {MP}}^{(a)}\) was 40.8% (95% CI by delta method: (0.27, 0.54)), implying that over 40% of the intervention effect was mediated by 6month visit adherence. In multivariateadjusted analyses (Scenario II), stronger NIE and TE effects were obtained, corresponding to odds ratios of 0.08 and 0.38, respectively and the multivariateadjusted MP estimate was also around 40%. The bootstrap and multivariate confidence intervals were very close in this analysis, although the width of the bootstrap confidence interval was slightly smaller than that using the delta method variance for NIE and TE, but slightly larger than the delta method for MP.
Because the outcome prevalence in the MaxART study is around 23%, the above analysis based on the rare outcome approximation may be biased, as suggested by our simulation study. Thus, we repeated the mediation analysis using the exact expressions given in Table 2. Generally speaking, the results of the mediation analysis accounting for common outcome prevalence were similar to those obtained using the rare outcome approximation. In both the steptimeadjusted model and the multivariateadjusted model, the adjusted \(\widehat {\text {TE}}\) were slightly weaker than previous results assuming a rare outcome (steptimeadjusted \(\widehat {\text {TE}}=0.24\) and multivariateadjusted \(\widehat {\text {TE}}= 0.10\) on the odds ratio scale, compared to 0.23 and 0.08, respectively). As a result, \(\widehat {\text {MP}}\) slightly increased in the adjusted analyses (steptimeadjust ed \(\widehat {\text {MP}}=44\%\), multivariateadjusted \(\widehat {\text {MP}}=42\%\)), compared to 41% and 40%, respectively, with the rare outcome approximation. In summary, for both the results based on and not based on the rare disease assumption, over 40% of the intervention effect on the 12month retention in care was mediated by 6month visit adherence.
Discussion
The difference and product methods are two popular approaches for estimating NIE and MP in mediation analysis [2]. While there has been comprehensive empirical evaluations of the difference method [21], there were only a few empirical evaluations of the product method as shown in Table 1. For this reason, we conducted a comprehensive simulation study to evaluate the performance of \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) obtained by the product method under various scenarios likely to be encountered in practice. We also provided the \(\widehat {\text {NIE}}\) and \(\widehat {\text {MP}}\) estimators without the rare outcome assumption for a binary outcome, and examined extent to which the current approximate mediation analysis was robust to violations of the rare outcome assumption.
The estimators investigated in our work has been implemented in an R package mediateP freely available at the Comprehensive R Archive Network (CRAN; https://cran.rproject.org); installations and instructions for the R package are given in Appendix A to facilitate their applications. Comparing to the current software for assessing conditional mediation measures (e.g., the SAS and SPSS macros given by [28] and the SAS macro and R package GEEmediate given by [21]), the mediateP package gives the exact NIE and MP estimates without rare outcome assumption, when a binary outcome is modeled by a logistic regression.
We demonstrated that \(\widehat {\text {NIE}}\) had very little bias and the variance estimate for \(\widehat {\text {NIE}}\) were quite close to the true values estimated under all scenarios considered from Monte Carlo simulations. In general, the multivariate delta method provided accurate variance estimates and valid interval estimates once the sample size was at least 500, and the bootstrap remained accurate even when the same size was even smaller. We found that larger sample sizes were needed to obtain valid MP point and interval estimates. Specifically, when the outcome was continuous, a sample size of 500 was required for valid point and interval estimates. In the binary outcome scenarios with a rare outcome, a sample size of 5000, and 200 cases or more, were required to obtain satisfactory MP point estimate and bootstrap interval estimates. We observed that the multivariate delta method provided valid MP confidence intervals when sample size ≥20,000 and number of cases ≥500 in binary outcome scenario.
We confirmed that the bootstrap method provided better interval estimates compared the multivariate delta method in smaller sample size scenarios, as may be found in some social science applications. However, the bootstrap method requires substantially more computational time to fit the mediation models in order to obtain the empirical \(\widehat {\text {NIE}}\) or \(\widehat {\text {MP}}\) distribution, which may be computationally burdensome in large epidemiological cohort studies. While we recommend bootstrap for the interval estimation with small sample size, when the sample size is larger, sample size ≥500 for the continuous outcome scenarios and sample size ≥20,000 and number of cases ≥500 for studies of binary outcome, we recommend the multivariate delta method for obtaining valid and computationally efficient confidence intervals. To facilitate application, our R package mediateP implements both variance estimators.
In addition, our simulation study also showed that the accuracy of \(\widehat {\text {MP}}\) also depends on the effect size of the TE. When the sample size is too small, a smaller TE is associated with a more biased MP point estimates and interval estimates with undercoverage, especially in binary outcome scenarios. In many epidemiological studies, when there is reason to believe that the NIE or NDE is not close to zero, a relatively smaller sample size may be adequate for obtaining valid point and interval estimates of the MP. For example, when the outcome is binary, we suggest that a sample size ≥20,000 and number of cases ≥500 is needed for the multivariate delta method to obtain satisfactory MP confidence intervals with close to nominal coverage rate. If there is reason to believe that the TE is not too small (TE ≥ log(1.5)), we found that the product method could accurately estimate MP with a sample size of at least 5000 and number of cases at least 150.
In the binary outcome scenarios (Cases #3 and #4), expressions of NIE^{(a)} and MP^{(a)} defined on a log odds ratio scale have been commonly used in biomedical and epidemiological studies [54–56]. Those expressions can be extended to include a log link function in the outcome model (2) and bypass the rare outcome assumption, and the corresponding mediation measure expressions can be defined on a log risk ratio scale. Using the logistic outcome model, our simulation study suggests that, when the outcome prevalence was less than 5%, the rare outcome assumption worked well for both \(\widehat {\text {NIE}}^{(a)}\) and \(\widehat {\text {MP}}^{(a)}\). When the outcome prevalence was greater than 5% and with a binary mediator (Case #4), the percent bias of \(\widehat {\text {NIE}}^{(a)}\) and \(\widehat {\text {MP}}^{(a)}\) was substantial (usually greater than 10%) and the exact expressions are recommended. However, with a binary outcome and continuous mediator (Case #3), we found that \(\widehat {\text {MP}}^{(a)}\) always provides satisfactory point and interval estimates even when the outcome was common.
Our simulation study has several limitations and future work is needed to supplement the conclusions in this manuscript. For simplicity, our simulation study did not consider the impact of unmeasured confounders, future research is needed to examine the robustness of product method in the presence of unmeasured confounders. In addition, when the outcome is binary, the mediation measures can be defined on the odds ratio scale [2]. Web Appendix E (see Additional file 1) show the relationship between mediation measures defined on a log odds ratio scale and odds ratio scale. In summary, we have found that asymptotic inference performs well for the product method in sample sizes typically found in epidemiology and public health settings. In addition, for common binary outcomes, exact expressions are needed to obtain unbiased estimates and strategies for point and variance estimation have been provided here.
Conclusions
We conducted simulation studies to examine the empirical performance of product method for calculating the point and interval estimates of two commonly used mediation measures: the natural indirect effect (NIE) and mediation proportion (MP). Our simulations confirm that the proposed multivariate delta method based on the asymptotic theory is valid for constructing confidence interval under sample sizes commonly encountered in public health. Given that the multivariate delta method is computationally more efficient than bootstrap, we therefore recommend to use the multivariate delta method to calculate the confidence interval when sample size is adequate (sample size ≥500 in scenarios of continuous outcome and sample size ≥20,000 number of cases ≥500 in scenarios of binary outcome). The rare outcome assumption in the scenarios of binary outcome generally performs well when baseline outcome prevalence less than 5% but can result in substantial bias when outcome prevalence is less rare. The exact estimator is a more appropriate choice for a common outcome.
Appendix a: instructions for the mediateP package
The mediateP package calculates the point and interval estimates for the NIE, TE and MP, based on the product method, as described in this paper. The source files for the mediateP package was provided on CRAN https://cran.rproject.org.
First, use the following statements to install the mediateP package
The main function of the "mediateP" package is mediate(), which provides the mediation analysis results. It can be called with,
The function has 14 arguments. These are

data= (Required) The name of the dataset.

outcome= (Required) Name of the outcome variable, which should be either a continuous or binary datatype.

mediator= (Required) Name of the mediator variable, which should be either a continuous or binary datatype.

exposure= (Required) Name of the exposure variable, which should be either a continuous or binary datatype.

binary.outcome= (Required) If the outcome is binary, set to 1. If the outcome is continuous, set to 0. The default value is 0.

binary.mediator= (Required) If the mediator is binary, set to 1. If the mediator is continuous, set to 0. The default value is 0.

covariate.outcome= A vector of names showing the confounding variables used in the outcome model. The default value is NULL, which represents no confounding variables. We only accepted continuous and binary confounding variables, if one confounding variable is categorical, please set it to a series of binary variables in advance.

covariate.mediator= A vector of names showing the confounding variables used in the mediator model. The default value is NULL, which represents no confounding variables. We only accepted continuous and binary confounding variables, if one confounding variable is categorical, please set it to a series of binary variables in advance.

x0= (Required) The baseline exposure level (i.e., x^{∗}). The default value is 0.

x1= (Required) The new exposure level (i.e., x). The default value is 1.

c.outcome= A vector of numbers representing the conditional level of the confounding variables in the outcome model. The default value is a vector of 0.

c.mediator= A vector of numbers representing the conditional level of the confounding variables in the mediator model. The default value is a vector of 0.

boot= If a percentile bootstrap confidence interval needed to be added, set to 1. Otherwise, set to 0. The default value is 0.

R= (Required if boot=1) The number of replications when apply the percentile bootstrap method to calculate the confidence interval. The default value is 2,000.
We now illustrate the usage of the mediate function. First, using the following statements to simulate a dataset
This dataset, named mydata, includes a continuous exposure (X), a binary mediator (M), a binary outcome (Y), as well as two confounding variables (C1 and C2). mydata has 2,000 observations, where the first 6 observations are shown as follows
We conducted a mediation analysis using mediate(). In the outcome model, we adjusted for C1 and C2. In the mediator model, we only adjusted for C1. We calculated the NIE, TE and MP for exposure in change from 0 to 1, conditional on C1=0 and C2=1, as follows
Finally, we got the following mediation analysis output
More illustrative examples under other datatypes can be found by using the syntax help(mediate).
Availability of data and materials
The R codes for simulation analyses are available from the first author on reasonable request. The MaxART dataset is available upon request from the second author.
Abbreviations
 NDE:

Natural direct effect
 NIE:

Natural indirect effect
 MP:

Mediation proportion
 TE:

Total effect
References
 1
Baron RM, Kenny DA. The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations. J Personal Soc Psychol. 1986; 51(6):1173.
 2
VanderWeele T. Explanation in Causal Inference: Methods for Mediation and Interaction. Oxford: Oxford University Press; 2015.
 3
Li H, Stampfer MJ, Mucci L, Rifai N, Qiu W, Kurth T, Ma J. A 25year prospective study of plasma adiponectin and leptin concentrations and prostate cancer risk and survival. Clin Chem. 2010; 56(1):34–43.
 4
Khan S, Spiegelman D, Walsh F, Mazibuko S, Pasipamire M, Chai B, Reis R, Mlambo K, Delva W, Khumalo G, et al.Early access to antiretroviral therapy versus standard of care among HIVpositive participants in Eswatini in the public health sector: the MaxART steppedwedge randomized controlled trial. J Int AIDS Soc. 2020; 23(9):25610.
 5
Robins JM, Greenland S. Identifiability and exchangeability for direct and indirect effects. Epidemiology. 1992:143–55.
 6
Pearl J. Direct and indirect effects. In: Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence, UAI’01. San Francisco, CA, USA: Morgan Kaufmann Publishers Inc.: 2001. p. 411–20.
 7
VanderWeele TJ. Policyrelevant proportions for direct effects. Epidemiology. 2013; 24(1):175.
 8
Bebu I, Braffett BH, PopBusui R, Orchard TJ, Nathan DM, Lachin JM. The relationship of blood glucose with cardiovascular disease is mediated over time by traditional risk factors in type 1 diabetes: the dcct/edic study. Diabetologia. 2017; 60(10):2084–91.
 9
Bowe B, Xie Y, Yan Y, Xian H, AlAly Z. Diabetes minimally mediated the association between pm 2.5 air pollution and kidney outcomes. Sci Rep. 2020; 10(1):1–9.
 10
Chang CH, Huang YF, Wang PW, Lai CH, Huang LW, Chen HC, Lin MH, Yang W, Mao IF, Chen ML. Associations between prenatal exposure to bisphenol a and neonatal outcomes in a taiwanese cohort study: Mediated through oxidative stress?Chemosphere. 2019; 226:290–7.
 11
Huang L, Wei Y, Shen S, Shi Q, Bai J, Li J, Qin S, Yu H, Chen F. Therapeutic effect of apatinib on overall survival is mediated by prolonged progressionfree survival in advanced gastric cancer patients. Oncotarget. 2017; 8(17):29346.
 12
Inzaule SC, Kityo CM, Siwale M, Akanmu AS, Wellington M, de Jager M, Ive P, Mandaliya K, Stevens W, Boender TS, et al.Previous antiretroviral drug use compromises standard firstline hiv therapy and is mediated through drugresistance. Sci Rep. 2018; 8(1):1–7.
 13
Parker SE, Collett BR, Speltz ML, Werler MM. Prenatal smoking and childhood behavior problems: is the association mediated by birth weight?J Dev Orig Health Dis. 2016; 7(3):273.
 14
Imai K, Keele L, Tingley D. A general approach to causal mediation analysis. Psychol Methods. 2010; 15(4):309.
 15
Imai K, Keele L, Yamamoto T. Identification, inference and sensitivity analysis for causal mediation effects. Stat Sci. 2010:51–71.
 16
Tchetgen EJT, Shpitser I. Semiparametric theory for causal mediation analysis: efficiency bounds, multiple robustness, and sensitivity analysis. Ann Stat. 2012; 40(3):1816.
 17
Tchetgen Tchetgen EJ, Shpitser I. Estimation of a semiparametric natural direct effect model incorporating baseline covariates. Biometrika. 2014; 101(4):849–64.
 18
Huber M, Lechner M, Strittmatter A. Direct and indirect effects of training vouchers for the unemployed. J R Stat Soc Ser A Ser A (Stat Soc). 2018; 181(2):441–63.
 19
MacKinnon DP, Warsi G, Dwyer JH. A simulation study of mediated effect measures. Multivar Behav Res. 1995; 30(1):41–62.
 20
VanderWeele TJ. Mediation analysis: a practitioner’s guide. Annu Rev Public Health. 2016; 37:17–32.
 21
Nevo D, Liao X, Spiegelman D. Estimation and inference for the mediation proportion. Int J Biostat. 2017; 13(2):20170006.
 22
Jiang Z, VanderWeele TJ. When is the difference method conservative for assessing mediation?Am J Epidemiol. 2015; 182(2):105–8.
 23
VanderWeele TJ, Vansteelandt S. Odds ratios for mediation analysis for a dichotomous outcome. Am J Epidemiol. 2010; 172(12):1339–48.
 24
Bollen KA, Stine R. Direct and indirect effects: Classical and bootstrap estimates of variability. Sociol Methodol. 1990:115–40.
 25
MacKinnon DP, Fairchild AJ, Fritz MS. Mediation analysis. Annu Rev Psychol. 2007; 58:593–614.
 26
MacKinnon DP, Lockwood CM, Williams J. Confidence limits for the indirect effect: Distribution of the product and resampling methods. Multivar Behav Res. 2004; 39(1):99–128.
 27
Shrout PE, Bolger N. Mediation in experimental and nonexperimental studies: new procedures and recommendations. Psychol Methods. 2002; 7(4):422.
 28
Valeri L, VanderWeele TJ. Mediation analysis allowing for exposure?mediator interactions and causal interpretation: theoretical assumptions and implementation with SAS and SPSS macros. Psychol Methods. 2013; 18(2):137.
 29
Gaynor SM, Schwartz J, Lin X. Mediation analysis for common binary outcomes. Stat Med. 2019; 38(4):512–29.
 30
Samoilenko M, Blais L, Lefebvre G. Comparing logistic and logbinomial models for causal mediation analyses of binary mediators and rare binary outcomes: evidence to support crosschecking of mediation results in practice. Obs Stud. 2018; 4:193–216.
 31
Doretti M, Raggi M, Stanghellini E. Exact parametric causal mediation analysis for a binary outcome with a binary mediator. Stat Methods Appl. 2021:1–22.
 32
Barfield R, Shen J, Just AC, Vokonas PS, Schwartz J, Baccarelli AA, VanderWeele TJ, Lin X. Testing for the indirect effect under the null for genomewide mediation analyses. Genet Epidemiol. 2017; 41(8):824–33.
 33
Biesanz JC, Falk CF, Savalei V. Assessing mediational models: Testing and interval estimation for indirect effects. Multivar Behav Res. 2010; 45(4):661–701.
 34
Fritz MS, MacKinnon DP. Required sample size to detect the mediated effect. Psychol Sci. 2007; 18(3):233–9.
 35
MacKinnon DP, Lockwood CM, Hoffman JM, West SG, Sheets V. A comparison of methods to test mediation and other intervening variable effects. Psychol Methods. 2002; 7(1):83.
 36
Rijnhart JJ, Twisk JW, Eekhout I, Heymans MW. Comparison of logisticregression based methods for simple mediation analysis with a dichotomous outcome variable. BMC Med Res Methodol. 2019; 19(1):1–10.
 37
VanderWeele TJ, Vansteelandt S. Conceptual issues concerning mediation, interventions and composition. Stat Interface. 2009; 2(4):457–68.
 38
Cole SR, Frangakis CE. The consistency statement in causal inference: a definition or an assumption?Epidemiology. 2009; 20(1):3–5.
 39
Pearl J. The causal mediation formula–a guide to the assessment of pathways and mechanisms. Prev Sci. 2012; 13(4):426–36.
 40
Spiegelman D, VanderWeele TJ. Evaluating public health interventions: 6. modeling ratios or differences? let the data tell us. Am J Public Health. 2017; 107(7):1087–91.
 41
Pandey D, Shrivastava P. Mediation effect of social support on the association between hardiness and immune response. Asian J Psychiatr. 2017; 26:52–55.
 42
Pandey D, Shrivastava P. Psychometric properties and confirmatory factor analysis of the social support scale. Int J Indian Psychol. 2016; 3(4):191–8.
 43
VanderWeele TJ, Asomaning K, Tchetgen Tchetgen EJ, Han Y, Spitz MR, Shete S, Wu X, Gaborieau V, Wang Y, McLaughlin J, et al.Genetic variants on 15q25. 1, smoking, and lung cancer: an assessment of mediation and interaction. Am J Epidemiol. 2012; 175(10):1013–20.
 44
Liu Q, Pierce DA. A note on gauss–hermite quadrature. Biometrika. 1994; 81(3):624–9.
 45
Lin D, Zeng D. Proper analysis of secondary phenotype data in casecontrol association studies. Genet Epidemiol: Off Publ Int Genet Epidemiol Soc. 2009; 33(3):256–65.
 46
Li Y, Schneider JA, Bennett DA. Estimation of the mediation effect with a binary mediator. Stat Med. 2007; 26(18):3398–414.
 47
Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986; 73(1):13–22.
 48
Wang YG, Carey V. Working correlation structure misspecification, estimation and covariate design: implications for generalised estimating equations performance. Biometrika. 2003; 90(1):29–41.
 49
Oehlert GW. A note on the delta method. Am Stat. 1992; 46(1):27–9.
 50
Davison AC, Hinkley DV, Vol. 1. Bootstrap Methods and Their Application. Cambridge: Cambridge university press; 1997.
 51
Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Boca Raton: CRC press; 1994.
 52
Walsh FJ, Bärnighausen T, Delva W, Fleming Y, Khumalo G, Lejeune CL, Mazibuko S, Mlambo CK, Reis R, Spiegelman D, et al.Impact of early initiation versus national standard of care of antiretroviral therapy in swaziland’s public sector health system: study protocol for a steppedwedge randomized trial. Trials. 2017; 18(1):1–10.
 53
Groenwold RH, White IR, Donders ART, Carpenter JR, Altman DG, Moons KG. Missing covariate data in clinical research: when and when not to use the missingindicator method for analysis. CMAJ. 2012; 184(11):1265–9.
 54
Agerbo E, Sullivan PF, Vilhjalmsson BJ, Pedersen CB, Mors O, Børglum AD, Hougaard DM, Hollegaard MV, Meier S, Mattheisen M, et al.Polygenic risk score, parental socioeconomic status, family history of psychiatric disorders, and the risk for schizophrenia: a danish populationbased study and metaanalysis. JAMA Psychiatry. 2015; 72(7):635–41.
 55
Dadvand P, Ostro B, Figueras F, Foraster M, Basagaña X, Valentín A, Martinez D, Beelen R, Cirach M, Hoek G, et al.Residential proximity to major roads and term low birth weight: the roles of air pollution, heat, noise, and roadadjacent trees. Epidemiology. 2014:518–25.
 56
Scott R, Langenberg C, Sharp S, Franks P, Rolandsson O, Drogan D, van der Schouw Y, Ekelund U, Kerrison N, Ardanaz E, et al.The link between family history and risk of type 2 diabetes is not explained by anthropometric, lifestyle or genetic risk factors: the epicinteract study. Diabetologia. 2013; 56:60–9.
Acknowledgements
Not applicable.
Funding
This work was funded by NIH grant DP1ES025459. The funding body had no role in the design, collection, analysis and interpretation of the data, and writing of the manuscript.
Author information
Affiliations
Contributions
CC: programming, interpretation of results, writing of manuscript. DS and FL: results checking, interpretation of results, writing of manuscript. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
The MaxART study was approved by the Swaziland National Health Research Review Board in July 2014 (Reference Number: MH/599C/FWA 000 15267). All methods were carried out in accordance with the relevant guidelines and regulations.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1
“Web Material.pdf”: online supplementary material for “Estimating the natural indirect effect and the mediation proportion via the product method”, including Web Appendices A–D, Web Tables 1–8 and Web Figures 1–4 not shown in the manuscript.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Cheng, C., Spiegelman, D. & Li, F. Estimating the natural indirect effect and the mediation proportion via the product method. BMC Med Res Methodol 21, 253 (2021). https://0doiorg.brum.beds.ac.uk/10.1186/s12874021014254
Received:
Accepted:
Published:
DOI: https://0doiorg.brum.beds.ac.uk/10.1186/s12874021014254
Keywords
 Estimating equations
 Mediation analysis
 Natural indirect effect
 Total effect
 Product method
 Asymptotically uncorrelated