Estimating the natural indirect effect and the mediation proportion via the product method.

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 finite-sample 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 closed-form 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 non-negligible 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 exposure-disease 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 anti-retroviral 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 HIV-positive 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 exposure-disease 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 finite-sample 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 weighting-based 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 non-continuous 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 closed-form 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 closed-form variance and interval estimators through the method of generalized estimating equations and the multivariate delta method. We provided a comparison between the closed-form 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 closed-form 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.

Table 1

A comparison of the current work with several previous literature evaluating the empirical performance of the product method in mediation analysis under four data types: Case #1, continuous outcome and continuous mediator; Case #2, continuous outcome and binary mediator; Case #3, binary outcome and continuous mediator; and Case #4, binary outcome and binary mediator

	Natural Indirect Effect							Mediation Proportion
Literatures	Case #1	Case #2	Case #3			Case #4		Case #1	Case #2	Case #3			Case #4
			Approx.	Exact	P.A.	Approx.	Exact			Approx.	Exact	P.A.	Approx.	Exact
Current work	B.V.I	B.V.I	B.V.I	B.V.I		B.V.I	B.V.I	B.V.I	B.V.I	B.V.I	B.V.I		B.V.I	B.V.I
Barfield et al. (2017)	T	T	T			T
Biesanz, Falk, and Savalei (2010)	I.T
Fritz and MacKinnon (2007)	T
Gaynor et al. (2019)					B.I		B.I
MacKinnon, Warsi, and Dwyer (1995)	B.V							B.V
MacKinnon et al. (2002)	T
MacKinnon, Lockwood, and Williams (2004)	I
Rijnhart et al. (2019)			B.V			B.V				B.V			B.V
Samoilenko, Blais, and Lefebvre (2018)						B.V.I	B.V.I

1Note: The B, V, I, and T denote the bias, variance, confidence interval, and hypothesis testing, respectively. If one of those appears in one cell, it indicates that this operating characteristic has been covered in this literature. In Cases #3 and #4, Approx., Exact, and P.A. denote the approximate expression, exact expression, and the probit approximation expression, respectively (See Table 2 for their specific formulas)

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 12-month retention is mediated by the 6-month 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 , a vector of covariates, associated with outcome and measured before the exposure, some of which may be confounders of the estimated exposure-outcome association and/or the mediator-outcome 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 and M as the outcome Y and mediator M, respectively, that would have been observed when setting X=x. Similarly, let Y be the value of outcome Y that would have been observed when setting X=x and M=m.

Fig. 1

Mediation directed acyclic graph, where Y, X, M and denote the outcome, exposure, mediator, and confounders of the exposure-outcome and exposure-mediator relationships. The NIE of exposure X on outcome Y through mediator M is highlighted in blue and the NDE of exposure X on outcome Y is highlighted in red

Based on the classical counterfactual framework, we made several standard assumptions. The first assumption is consistency, which assumes that (i) the potential outcome Y and potential mediator M equal their respective observed variables Y and M when we set X=x, and (ii) the potential outcome Y 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 well-defined 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 , 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|, (A.2) Y⊥M|,X, (A.3) M⊥X|, and (A.4) for all x, x∗, and m. That is, (A.1)-(A.3) assume that the exposure-outcome, mediator-outcome, and exposure-mediator relationships are not confounded conditional on covariates . (A.4) is sometimes termed cross-world independence, which stresses that none of the confounders in the mediator-outcome relationship can be affected by exposure X [31]. Intuitively, (A.1)-(A.4) requires the investigators to measure all pre-exposure covariates that may confound the exposure-mediator-outcome mechanisms.

Under this framework, when changing the exposure levels from x∗ to x conditional on =, Nevo et al. [21] defined the NIE and NDE on a g-function scale:

where g(.) is a pre-specified 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 to 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 . As discussed in [14, 15], under assumptions (A.1)–(A.4), the NIE, NDE, and TE are given by , , and , 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 =, 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+ T-lymphocyte 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 . 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 M|X, is not normally distributed or Var(M|X,) 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 M|X, in model (4) follows a normal distribution with a constant variance σ2. Specifically, define an integrand we have that

both of which involve one-dimensional logistic-normal integrals that do not have closed-form solutions. Gaynor et al. (2019) [29] uses a probit function to approximate the logistic function in the integral and obtained closed-form 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 logistic-normal integrals, we consider in this paper the Gauss-Hermite Quadrature (GHQ) approach [44] to numerically calculate the integrals.

Table 2

Expressions of mediation measures under four different datatypes of the outcome and mediator. Case #1, continuous outcome and continuous mediator; Case #2, continuous outcome and binary mediator; Case #3, binary outcome and continuous mediator; and Case #4, binary outcome and binary mediator

Datatypes		NIE	NDE	Depend on W	Reference(s)
Case #1		β2γ1(x−x∗)	β1(x−x∗)	No	[28]
Case #2		β_2κ(x,w)−κ(x^∗,w)(1+κ(x,w))(1+κ(x^∗,w))	β1(x−x∗)	Yes	[32]
	Approx.	β2γ1(x−x∗)	β1(x−x∗)	No	[23]
Case #3	Exact	log∫mexp(β_2m)τ(x,x,m,w)dm∫mτ(x,x,m,w)dm−log∫mexp(β_2m)τ(x^∗,x,m,w)dm∫mτ(x^∗,x,m,w)dm	β_1(x−x^∗)+log∫mexp(β_2m)τ(x^∗,x,m,w)dm∫mτ(x^∗,x,m,w)dm−log∫mexp(β_2m)τ(x^∗,x^∗,m,w)dm∫mτ(x^∗,x^∗,m,w)dm	Yes	Web Appendix A
	Probit Approx.	logitΦsβ_0+sβ_1x+sβ3Tw+sβ_2(γ_0+γ_1x+γ2Tw)1+s^2β22σ^2−logitΦsβ_0+sβ_1x+sβ3Tw+sβ_2(γ_0+γ_1x^∗+γ2Tw)1+s^2β22σ^2	logitΦsβ_0+sβ_1x+sβ3Tw+sβ_2(γ_0+γ_1x^∗+γ2Tw)1+s^2β22σ^2−logitΦsβ_0+sβ_1x^∗+sβ3Tw+sβ_2(γ_0+γ_1x^∗+γ2Tw)1+s^2β22σ^2	Yes	[29]
Case #4	Approx.	β_2κ(x,w)−κ(x^∗,w)(1+κ(x,w))(1+κ(x^∗,w))	β1(x−x∗)	Yes	[2]
	Exact	log1+e^β_2η(x,w)+κ(x^∗,w)(1+η(x,w))1+e^β_2η(x,w)+κ(x,w)(1+η(x,w))−log1+e^β_2η(x,w)+e^β_2κ(x,w)(1+η(x,w))1+e^β_2η(x,w)+e^β_2κ(x^∗,w)(1+η(x,w))	β_1(x−x^∗)+log1+e^β_2η(x^∗,w)+κ(x^∗,w)(1+η(x^∗,w))1+e^β_2η(x,w)+κ(x^∗,w)(1+η(x,w))+log1+e^β_2η(x,w)+e^β_2κ(x^∗,w)(1+η(x,w))1+e^β_2η(x^∗,w)+e^β_2κ(x^∗,w)(1+η(x^∗,w))	Yes	[29] and Web Appendix B

1Note: NIE and NDE denote the natural indirect effect and natural direct effect, respectively, which are defined for X in change from x∗ to x conditional on =, on an identity scale in Cases #1 and #2 and a log odds ratio scale in Cases #3 and #4. Given NIE and NDE, the mediation proportion (MP) can be obtained by . In the probit approximation method, s=1/1.6, , and Φ(.) is the cumulative density function for the standard normal distribution

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 . Here we can obtain simple and closed-form 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 logistic-normal integrals in the mediation expressions can be approximated by the log-normal integrals that have closed-form 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(, NDE(, TE(, and MP(, 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 . 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 .

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 12-month retention is mediated by 6-month visit adherence in the MaxART study. Here, the 12-month retention and 6-month 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 . Given NIE and NDE, the MP is given by . If the outcome is rare, the following approximate NIE and NDE have been widely used [2]:

and NDE(=β1(x−x∗). As a result, . 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 , and unknown parameters in the mediator model, . 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 [,,σ2] for the exact expressions in Case #3 and [,] for other cases. Hereafter, we will rewrite the expressions of NIE and MP as and , respectively, to emphasize that those expressions are functions of .

In practice, the coefficients in the outcome model, , can be obtained by solving the following generalized estimating equation (GEE) [47]:

where are n observations of {Y,X,M,} and V is a working variance term for the outcome Y. The optimal asymptotic efficiency of will be obtained when V=Var(Y). When V is misspecified, the resultant is still consistent under mild regularity conditions but could be less efficient [48]. Similarly, the coefficients in the mediator model, , can be obtained by solving the GEE below

where is the working variance for M. Also, misspecification of only impacts the efficiency of . When evaluating the exact expressions in Case #3, is needed and can be estimated by solving the GEE After obtaining through the above estimating equations, we can calculate and by plugging those parameter estimates in their expressions introduced in the previous sections; i.e., and . In what follows, we develop the closed-form asymptotic variance expressions of and 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 variance-covariance matrix of , abbreviated by , based on the theory of estimating equations [47]. Specifically, , , and can be estimated by the robust sandwich variance estimators [47] based on their respective estimating equations. If the working variance terms, V and , are correctly specified when estimating and , one can also approximate and by the negative inverse information matrix, and , respectively. In general, to obtain , we also need to estimate the covariances between , , and . 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 , , and obtained by solving U()=, U()=, 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 for the exact expressions in Case #3 and for the remaining Cases. Finally, the variances of and 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 and , respectively.

Alternatively, the nonparametric percentile bootstrap approach [50, 51] can be used to approximate the empirical distributions of by resampling the dataset with replacement and re-estimating 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 finite-sample 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 closed-form 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+γ1X,1), where γ0=0. γ1 was chosen to 0.408, corresponding to the exposure-mediator correlation, Corr(X,M)=0.2. Finally, given X and M, we simulated Y∼N(β0+β1X+β2M,1), where β0 was fixed as 0. We let β1=(1−MP)×TE and 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=1|X))=γ0+γ1X. Noting that , we chose the values of γ0 such that the baseline prevalence of the mediator P(M=1|X=0)=0.2. Similar to Case #1, we chose γ1=0.903 to give an exposure-mediator correlation 0.2. Finally, we simulated Y∼N(β0+β1X+β2M,1), where β0=0, β1=(1−MP)×TE from the definition NDE=(1−MP)×TE=β1, and by noticing .

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 , we selected β0 so the baseline outcome prevalence was 3%. Then, we selected β1 and β2 by numerically solving the system of equations and , where and 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 and , where and refer to the exact expressions of NIE and NDE, given in (8) and (9).

For each data type, we obtained and , the variance estimates of and 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( and MP(), 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.

Table 3

Simulation results for Case #1: continuous outcome and continuous mediator

			NIE^				MP^
N	MP	TE	Bias(%)	CR (d)	CR (b)	VR	Bias(%)	CR (d)	CR (b)	VR
150	0.05	0.25	-16.4	99.2	96.4	0.998	-27.5	99.2	98.8	0.000
		0.5	-15.1	96.5	95.7	0.992	-14.1	96.9	96.5	0.058
		1	-10.5	92.2	94.7	0.993	-9.8	92.6	95.1	0.949
	0.2	0.25	-10.5	92.2	94.7	0.993	-19.7	90.6	97.3	0.000
		0.5	-6.2	92.5	94.8	1.001	-2.8	92.9	96.4	0.092
		1	-2.1	94.8	94.9	1.012	-1.7	95.7	95.4	0.944
	0.5	0.25	-4.7	93.2	94.8	1.007	-14.1	89.1	97.0	0.000
		0.5	-1.6	95.1	95.0	1.012	-0.4	96.0	97.0	0.026
		1	-0.1	95.5	94.9	1.013	0.2	97.5	95.3	0.710
500	0.05	0.25	-4.4	96.0	94.4	0.997	-3.5	97.8	96.6	0.002
		0.5	-3.5	94.7	94.2	0.988	-3.4	95.3	94.5	0.855
		1	-2.6	93.1	94.2	0.983	-2.8	93.4	94.4	0.948
	0.2	0.25	-2.6	93.1	94.2	0.983	-1.5	92.1	97.0	0.001
		0.5	-1.5	93.6	94.5	0.982	-0.8	94.4	94.8	0.801
		1	-0.2	94.5	94.5	0.981	-0.1	94.7	94.6	0.945
	0.5	0.25	-0.9	94.0	94.5	0.980	-0.4	92.3	97.5	0.000
		0.5	-0.1	94.7	94.3	0.982	0.1	95.5	94.9	0.770
		1	0.6	94.5	94.5	0.983	0.0	95.6	94.7	0.915
1000	0.05	0.25	-4.0	95.4	94.7	0.969	-3.2	97.0	95.4	0.724
		0.5	-2.2	94.8	94.7	0.969	-2.8	95.2	95.1	0.924
		1	-2.1	94.1	94.8	0.971	-1.9	94.0	94.7	0.967
	0.2	0.25	-2.1	94.1	94.8	0.971	-1.1	93.5	96.3	0.580
		0.5	-1.1	94.2	94.6	0.981	-1.2	94.2	94.8	0.923
		1	-0.9	94.9	94.8	0.986	-0.8	94.6	94.4	0.981
	0.5	0.25	-1.1	94.4	94.6	0.983	-0.5	92.8	96.6	0.441
		0.5	-0.5	95.0	94.9	0.987	-0.5	94.9	94.6	0.916
		1	-0.3	94.8	94.7	0.988	-0.4	95.0	94.5	0.966
5000	0.05	0.25	-0.5	94.8	94.8	0.990	0.3	95.5	94.8	0.961
		0.5	-0.2	94.9	95.0	0.995	-0.3	95.4	95.0	0.990
		1	-0.2	95.0	95.1	1.007	-0.3	95.0	95.1	1.006
	0.2	0.25	-0.2	95.0	95.1	1.007	0.1	95.6	95.1	0.936
		0.5	-0.1	95.0	95.2	1.027	-0.1	95.1	95.3	1.004
		1	0.0	95.3	95.2	1.038	0.0	95.1	95.1	1.031
	0.5	0.25	0.0	95.3	95.5	1.032	0.1	95.0	94.7	0.920
		0.5	-0.1	95.3	95.2	1.040	0.0	95.2	95.0	0.995
		1	0.0	95.4	95.2	1.041	0.0	95.3	95.4	1.022

1Note: Bias(%), CR (, CR (, and VR denote the median percent bias, 95% confidence interval coverage rate of multivariate delta method, 95% confidence interval coverage rate of percentile bootstrap method, and mediation variance ratio, respectively. The coverage rates outside the 95% confidence boundary, i.e., , were highlighted in bold, where q denotes the nominal confidence interval threshold (95%) and B denotes the number of replication (5,000). The median percent bias was calculated as the median of the ratio of bias to the true value over 5,000 replications, i.e., , where p denotes the true value of the causal mediation measure, and is the point estimate of the simulated causal mediation measure. The median variance ratio is defined by the ratio of median delta-method variance estimators across 5,000 replications to the empirical variance of causal mediation measure estimates from the 5,000 replications

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 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 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., ). 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 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 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 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 ( v.s. and v.s. ) 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 small-sample 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, and , provided accurate point and interval estimates when the outcome prevalence >1%. When TE= log(1.2) and the baseline outcome prevalence was 1%, and 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 diverged from 0 and the coverage rate of by the delta method substantially decreased from 95%. Similarly, was also quite different from its true value when the baseline outcome prevalence ≥10% (See Web Figure 1 in additional file 1).

Fig. 2

Performance of MP ( estimates (black line) and MP estimates (blue dotted line) when changing baseline outcome prevalence from 1% to 50% in Case #4, where sample size is 20,000. Bias(%), CR (, and VR denote the percent bias, coverage rate by the multivariate delta method, and variance ratio. Upper row: results for TE=log(1.2) and MP=0.1; second row: results for TE=log(1.2) and MP=0.5; third row: results for TE=log(2) and MP=0.1; bottom row: results for TE=log(2) and MP=0.5

With a binary outcome and continuous mediator (Case #3), and 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 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 was sensitive to the baseline prevalence and the magnitude of MP. For example, when MP=0.5, the percent bias of increased as baseline outcome prevalence increased, and the confidence interval coverage rates rapidly declined. When the true MP=0.1, the percent bias of was negligible over the range of outcome prevalences considered.

Homoscedasticity and normality assumptions in case #3

In the above simulations, we simulated M|X 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 , where . Here, η1 measures the extent to which the variance depends upon the exposure and η0 was chosen such that the expectation of 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+γ1X+ε, where ε follows a gamma distribution. Specifically, we let , where b follows a gamma distribution with density (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 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 confounder-exposure, confounder-mediator, and confounder-outcome 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 stepped-wedge cluster randomized trial among HIV-positive 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 pre-specified 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 intervention-retention relationship is unknown. In this illustrative example, we investigated the extent to which the effect of the intervention (SoC v.s. EAAA) on 12-month 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 12-month 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 6-month 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 follow-up 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 6-month 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., ) 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 intervention-retention relationship, visit adherence-retention relationship, or intervention-visit 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/m2), 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, 1-3 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 non-retention 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 steptime-adjusted model, we found that the intervention was protective against 12-month non-retention, with odds ratios of 0.23 and 0.55 for and 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 steptime-adjusted was 40.8% (95% CI by delta method: (0.27, 0.54)), implying that over 40% of the intervention effect was mediated by 6-month visit adherence. In multivariate-adjusted analyses (Scenario II), stronger NIE and TE effects were obtained, corresponding to odds ratios of 0.08 and 0.38, respectively and the multivariate-adjusted 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.

Table 4

Mediation analysis of MaxART [4]. (n=1731)

Expression	Scenario	Parameter	Point	S.E.	Delta 95% CI	Bootstrap 95% CI
Steptime adjusted	Approximate	NIE (a)	-0.601	0.091	(-0.779,-0.424)	(-0.782,-0.445)
		TE (a)	-1.472	0.226	(-1.915,-1.030)	(-1.973,-1.031)
		MP (a)	0.408	0.069	(0.273,0.544)	(0.305,0.559)
	Exact	NIE	-0.630	0.093	(-0.813,-0.448)	(-0.816,-0.474)
		TE	-1.444	0.213	(-1.862,-1.027)	(-1.922,-1.023)
		MP	0.437	0.073	(0.293,0.580)	(0.327,0.597)
Multivariate adjusted	Approximate	NIE (a)	-0.972	0.121	(-1.208,-0.735)	(-1.287,-0.775)
		TE (a)	-2.520	0.287	(-3.082,-1.958)	(-3.282,-1.975)
		MP (a)	0.386	0.050	(0.288,0.483)	(0.292,0.494)
	Exact	NIE	-0.970	0.120	(-1.205,-0.735)	(-1.282,-0.772)
		TE	-2.316	0.267	(-2.841,-1.792)	(-3.033,-1.819)
		MP	0.419	0.054	(0.313,0.525)	(0.322,0.539)

1Note: All the mediation measures, including NIE, TE, and MP, are defined on a log odds ratio scale for the intervention in change from SoC to EAAA, conditional on the most frequent level of the confounding variables. S.E. denotes the standard error of the point estimates, which is calculated by the multivariate delta method. We implemented the delta method and bootstrap method with 1,000 replications to calculate the 95% confidence interval (95% CI) of each mediation measure

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 steptime-adjusted model and the multivariate-adjusted model, the adjusted were slightly weaker than previous results assuming a rare outcome (steptime-adjusted and multivariate-adjusted on the odds ratio scale, compared to 0.23 and 0.08, respectively). As a result, slightly increased in the adjusted analyses (steptime-adjust ed , multivariate-adjusted ), 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 12-month retention in care was mediated by 6-month 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 and obtained by the product method under various scenarios likely to be encountered in practice. We also provided the and 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.r-project.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 had very little bias and the variance estimate for 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 or 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 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 under-coverage, 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( and MP( 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 and . When the outcome prevalence was greater than 5% and with a binary mediator (Case #4), the percent bias of and 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 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.r-project.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).

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.