Skip to content
Publicly Available Published by De Gruyter September 20, 2017

Estimation and Inference for the Mediation Proportion

  • Daniel Nevo EMAIL logo , Xiaomei Liao and Donna Spiegelman

Abstract

In epidemiology, public health and social science, mediation analysis is often undertaken to investigate the extent to which the effect of a risk factor on an outcome of interest is mediated by other covariates. A pivotal quantity of interest in such an analysis is the mediation proportion. A common method for estimating it, termed the “difference method”, compares estimates from models with and without the hypothesized mediator. However, rigorous methodology for estimation and statistical inference for this quantity has not previously been available. We formulated the problem for the Cox model and generalized linear models, and utilize a data duplication algorithm together with a generalized estimation equations approach for estimating the mediation proportion and its variance. We further considered the assumption that the same link function hold for the marginal and conditional models, a property which we term “g-linkability”. We show that our approach is valid whenever g-linkability holds, exactly or approximately, and present results from an extensive simulation study to explore finite sample properties. The methodology is illustrated by an analysis of pre-menopausal breast cancer incidence in the Nurses’ Health Study. User-friendly publicly available software implementing those methods can be downloaded from the last author’s website (SAS) or from CRAN (R).

1 Introduction

In many public health, biological, and biomedical systems, the mechanism that explains how an intervention or exposure affects the outcome of interest is unknown, even after a causal association between the exposure and the outcome is established. It is sometimes hypothesized that there exists a mediator that connects the exposure and the outcome, sitting on the causal pathway between the exposure and the outcome. In observational studies, identifying a plausible ideally pre-specified, mediator can strengthen the casual inference of the findings. For example, in an evaluation of the effectiveness of the ongoing, trillion dollar President’s Emergency Plan for AIDS Relief (PEPFAR) in reducing HIV incidence and prevention in sub-Saharan Africa, it would strengthen the evidence of a causal inference if it could be shown that a substantial proportion of the reduction in disease incidence in time was mediated by increased programmatic coverage in the region, thus diminishing exogenous time trends as the best explanation for any observed decline.

Several methods have been proposed to assess whether mediation exists and to quantify its magnitude [1, 2, 3, 4]. [5] described a sequence of hypothesis tests to assess the evidence in the data for mediation by a specific covariate. They assumed a linear model for the relationship between the outcome and the exposure, both marginally and conditionally on the mediator. They also assumed a linear model for the relationship between the exposure and the mediator. Within the counterfactual framework, the building blocks of mediation analysis are the natural direct effect (NDE), defined as the effect on the outcome when increasing the value of the exposure in one unit while holding the mediator at a fixed level, and the natural indirect effect (NIE), which is the effect on the outcome when the exposure is held fixed but the mediator value is changed as it would have been changed if the exposure value were increased by one unit [6, 7, 8]. The sum of the NDE and NIE is the total effect (TE). Under this framework, estimation methods for the TE, NDE and NIE were developed for various statistical models. Other examples include logistic regression [9], zero-inflated regression models [10] and high-dimensional mediators in linear regression with normal errors [11].

One way to estimate mediation is through the product method [5]. Another widely used method for assessing mediation is the difference method [5, 12, 13]. It quantifies the difference in estimates obtained from separate exposure-outcome relationship models, with and without the mediator. The mediation proportion, defined as the change in the effect of the exposure due to mediation by the mediator relative to the total effect, is a main parameter of interest when performing mediation analysis. An analogous measure in surrogacy analysis, termed proportion of treatment effect (PTE), aids researchers in deciding whether an intermediate marker can be used as a surrogate for a final outcome of interest. Quantifying PTE entails statistical questions relevant to those that arise in studying the mediation proportion. When the intermediate and the final outcomes are both binary, confidence intervals for the PTE were developed [14]. A time-to-event final outcome with surrogate biomarkers was also considered by [15], who used a data duplication algorithm in order to estimate the covariance between estimators obtained by separate models. The PTE measure in surrogacy research is still actively used and researched [16, e.g.,].

Methods for variance estimation, statistical testing and confidence interval construction of mediation parameters have been suggested by past authors. For the difference method, [14] suggested to use results from the linear model in binary outcome setup to approximate the covariance between the two estimators. Other approximation were described and compared in [4]. For the product method, variance can calculated either using delta method [17] or Goodmans exact variance-of-product formula [18]. However, since the finite samples behavior of the product method estimator was found to be nonnormal, the bootstrap was generally recommended [2, 19], at least for medium or small sample size. For the simulation-based mediation approach proposed by [20], a quasi-Bayesian Monte Carlo method or the bootstrap can be used [20, 21].

While the mediation proportion is a popular measure in mediation analysis [22, 23, 24, 25, e.g.,], statistical inference for this parameter is not sufficiently developed. The NIE and NDE are well-defined concepts, however, in practice, researchers are often primarily interested in the mediation proportion, as exemplified by the aforementioned papers. In this paper, we provide a framework for mediation analysis in generalized linear models (GLMs). We combine a generalized estimation equations (GEE) approach together with a data duplication algorithm to formulate valid statistical inference under minimal assumptions on the marginal and conditional distribution of the outcome. We discuss situations in which these assumptions should hold, and assess robustness to departures from these assumptions in extensive simulation studies. This paper further provides methods for statistical inference in mediation analysis using the difference method, including studying confidence intervals for the mediation proportion and hypothesis tests. Our investigation of these aspects is expanded beyond GLMs to inference about the mediation proportion for Cox model.

The reminder of this paper is organized as follows. In Section 2, we formulate the models needed for the estimation of the mediation proportion in GLMs. In Section 3, we consider the g-linkability property for common link functions. In Section 4, we present methods for inference for this parameter using the multivariate delta method and a data duplication method that enables consistent variance estimation. In Section 5, we present results from a simulation study. In Section 6, we illustrate the use of the methodology developed in studying mediation of the effect of risk factors for pre-menopausal breast cancer incidence by mammographic density in the Nurses’ Health Studies NHSI and NHSII. In Section 7, we discuss results and related issues. We describe the software we have made publicly available in the Appendix.

2 The models

Assume Y1,...,Yn is a sample of results of an outcome of interest, and that for each subject i we also observe a vector of factors Zi=(Xi,Mi,Wi) where Xi,Mi and Wi are an exposure of interest, a mediator and a vector of confounders, respectively. We assume the conditional mean function for the outcome is E(Yi|Zi)=g1(ZiTβ) with g being the link function and where β is an unknown parameter vector. A consistent estimator, βˆ, for β is obtained as the solution to the estimating equations

(1)U(β)=i=1nDivi1[yiE(Yi|Zi)]=0

where Di=E(Yi|Zi)/β and vi is the working variance of yi. By GEE theory, the variance of βˆ can be consistently estimated by the robust sandwich estimator [26, 27].

In this paper, we consider a mediator, M, which may be binary, multilevel, continuous or a vector of any one of these. For convenience, we henceforth treat M as a scalar. Nevertheless, our framework also applies to the investigation of the joint mediation effect of multiple covariates. Consider the following conditional and marginal mean models for Y, with respect to M

(2)E(Y|X,M,W)=g1(β0+β1X+β2M+β3TW)
(3)E(Y|X,W)=g1(β0+β1X+β3TW).

Let B=(β0,β1,β2,β3) and B=(β0,β1,β3) be the vectors of conditional and marginal regression model parameters, and denote Bˆ and Bˆ for their estimators obtained by solving eq. (1) under models (2) and (3), separately, respectively. When the two models (2) and (3) both hold simultaneously, we say we have g-linkability.

The definitions of the NDE and NIE, as given by [7], use the counterfactual framework. Let Y(x,m) be the counterfactual outcome value when setting X=x and M=m, and let M(x) be the counterfactual mediator value when setting X=x. When comparing two exposure or treatment levels x and x, the TE, NDE and NIE can be then written [6, 7] as

TE(x,x)=E(Y(x))E(Y(x))=E(Y(x,M(x))E(Y(x,M(x)))NDE(x,x)=E(Y(x,M(x)))E(Y(x,M(x)))NIE(x,x)=E(Y(x,M(x)))E(Y(x,M(x))).

Under certain identifiability conditions, to estimate these effects, the mediation formula, given by [28, 29] and [6], can be used. Extensive work has been published on nonparametric and non-linear models [6, 20, 21, e.g., ].

In GLMs, alternative definitions for the NDE, NIE and TE have been proposed. For example, for a binary outcome following a logistic model, [30] defined the TE on the odds ratio scale as

ORTE(x,x|W)=P(Y(x)=1|W)/(1P(Y(x)=1|W))P(Y(x)=1|W)/(1P(Y(x)=1|W)),

and have shown that it can be decomposed to a product of a NDE and a NIE. Therefore, on the log odds ratio scale, the TE decomposes to the sum of the NDE and the NIE. See [30] and [31] for further details. While these definitions depend on the value of the confounders, they show that, for a rare outcome, under the logistic regression outcome model and a linear mediator-exposure regression model, estimates of the NDE and NIE can be obtained using estimates of the regression model parameters. This result also extends to the log link function in GLMs, without a rare outcome assumption.

In this paper, we consider the causal effects on the link function scale, and assuming no exposure-mediator interaction. That is, as discussed by [32], under the assumptions given below, the TE, NDE and NIE are

(4)TE(x,x|W)=g[E(Y(x)|W)]g[E(Y(x)|W)]NDE(x,x|W)=g[E(Y(x,M(x))|W)]g[E(Y(x,M(x))|W)]NIE(x,x|W)=g[E(Y(x,M(x))|W)]g[E(Y(x,M(x))|W)].

Throughout this paper, we assume that, after adjusting for measured confounders, there is no unmeasured confounding of the estimates of the exposure-outcome relationship, the mediator-outcome relationship or the exposure-mediator relationship. We also assume that confounders of the mediator-outcome relationship are unaffected by the exposure. Alternative identifiability assumptions have been given, e.g., in [6], which we will not consider further in this paper.

The mediation proportion, defined as NIE/TE, can be estimated using either the difference or the product method. The product method fits, in addition to model (2), a model for mediator-exposure relationship,

(5)E(M|X,W)=γ0+γ1X+γ2TW.

Under the aforementioned no-confounding assumptions, considering a unit change (i.e., xx=1) and when models (2) and (5) hold, the NIE is γ1β2, the NDE is β1 and the mediation proportion is γ1β2/(γ1β2+β1). The product method estimator for the mediation proportion is obtained by plugging in estimates of β2 and γ1 in this expression. Variance estimation and confidence intervals construction have been typically conducted by the bootstrap, due to the skewness of the finite sample distribution [2, 4].

In this paper, we focus on the difference method, for which we will develop asymptotic properties. Under g-linkability, the TE equals β1 and the NIE equals β1β1. Therefore, the mediation proportion equals to

p=β1β1β1=1β1β1.

If p(0,1], p can be interpreted as a proportion. The situation where p=0 corresponds to β1=β1, hence in this case M does not mediate the effect of X at all. On the other hand, if p=1 then the effect of X is fully mediated by M. Finally, if p[0,1], the NDE and NIE are in opposite directions. One can get a consistent estimate for p by simply plugging in the appropriate estimator from each model. That is, pˆ=1βˆ1βˆ1, where βˆ1 and βˆ1 are the appropriate components of Bˆ and Bˆ. Under g-linkability, this estimator is consistent by standard GEE theory and the general mapping theorem.

The question of mediation can also be investigated when the available data is survival data. A counterfactual framework for mediation analysis of survival data has been previously provided [33, 34, 35, 36]. [15] considered this question for the Cox model in the context of the PTE. First, as in [15], we define the following two models for the hazard function at time t, h(t), conditionally and marginally, with respect to M

(6)h(t|X,M,W)=λ0(t)exp(β0+β1X+β2M+β3TW)h(t|X,W)=λ0(t)exp(β0+β1X+β3TW),

where λ0(t) and λ0(t) are baseline hazard functions. [15] have shown that these two models cannot hold at same time. However, they claimed that if either β3 or Λ0(t)=0tλ0(s)ds are small, then model (6) is a good approximation to the true conditional model. The assumption that Λ0(t) is small is the rare outcome assumption. They confirmed this claim using a small scale simulation study. When (6) holds, approximately, the Cox model is approximately g-linkable. Thus, in addition to GLMs, we investigate in this paper estimation and inference for p in approximately g-linkable Cox models.

3 Further results on g-linkability

In this section, we consider the issue of when the full model (2) and the marginal model (3) both hold with the same function g exactly or approximately. Recall that g-linkability is sufficient for ensuring that pˆ, the point estimator of p, is consistent. This subject was also discussed in the context of random effects models [37], in which the authors showed, for each common statistical model, what random effect’s distribution would provide a g-linkable conditional mean model condition on, and marginal over, the random effect. If g-linkability does not hold, then pˆ converges to pp. However, if g-linkability holds approximately, as in the case of logistic regression under rare outcome assumption (see below, and [38]), then one may expect p to be close to p, as discussed in [15] for the Cox model.

We consider the three common link functions: identity, log and logit. For each of these functions we give a general condition for the distribution of M given X and W that ensures g-linkability, where in the logit link function a rare outcome assumption is also needed. Numerous detailed and practical examples that fulfill these conditions can be constructed. In practice, the difference method does not require fitting of the mediator-exposure relationship model, as noted by [39]. For the validity of the product method, however, this model has to be correctly specified.

3.1 Identity link function

Under the identity link function, models (2) and (3) simplify to

E(Y|X,M,W)=β0+β1X+β2M+β3TWE(Y|X,W)=β0+β1X+β3TW.

We now show that g-linkability holds whenever E(M|X,W) is a linear function of X and W. To see that, let E(M|X,W)=a+b1X+b3TW, for some a,b1 and b3. Then,

E(Y|X,W)=E(E(Y|X,M,W)|X,W)=β0+β1X+β2E(M|X,W)+β3TW=β0+β1X+β3TW

where β0=β0+β2a, β1=β1+β2b1 and β3=β3+β2b3.

3.2 Log link function

Under the log link function, g(u)=log(u), the mean models become

E(Y|X,M,W)=exp(β0+β1X+β2M+β3TW)E(Y|X,W)=exp(β0+β1X+β3TW)

and we have

E(Y|X,W)=E(E(Y|X,M,W)|X,W)=exp(β0+β1X+β3TW)×E[exp(β2M)|X,W].

Therefore, in the log link case, glinkability holds if the log of the moment generating function of M|X,W can be written as a linear function of X and W. That is, logE[exp(β2M)|X,W]=a+b1X+b3TW.

3.3 Logit link function

The issue of whether the logistic regression model holds for both the conditional and marginal models has been discussed in the literature [37, 38, 39]. The logit link function, defined as logit(p)=log(p/(1p)), is typically used when Y is binary. It is well known and readily seen that when the outcome is rare, the logit function is similar to the log function. Thus, under rare outcome scenario, one may expect that g-linkability holds approximately for the logit link function, as is typical in many epidemiologic and public health studies [30]. We empirically investigate the limits of the rare outcome assumption in Section 5.

4 Inference for the mediation proportion

For simplicity of presentation, we assume throughout this section that g-linkability holds. Then, pˆ=1βˆ1βˆ1, where βˆ1 and βˆ1 are the appropriate components of Bˆ and Bˆ defined in Section 2. By the aforementioned GEE theory together with the general mapping theorem, this estimator is consistent.

Asymptotic confidence intervals for p have been constructed previously using either Fieller’s theorem or the delta method [14, 15, 40]. Here, we consider the latter. As written in [15], by the delta method pˆ has an asymptotic normal distribution with variance equals to

(7)σpˆ2=σβˆ12(β1)2+β12σβˆ12(β1)42β1σβˆ1,βˆ1(β1)3,

where

σβˆ12=Var(βˆ1),σβˆ12=Var(βˆ1)andσβˆ1,βˆ1=Cov(βˆ1,βˆ1).

While σβˆ12 and σβˆ12 can be consistently estimated using the robust sandwich estimator [26] for each of the models (2) and (3) separately, it is not obvious how to estimate σβˆ1,βˆ1, the covariance of estimators obtained from two separate models. In Section 4.1, we propose a data duplication algorithm to estimate this quantity.

Assume now we have estimates σˆβˆ12,σˆβˆ12 and σˆβˆ1,βˆ1 for σβˆ12,σβˆ12 and σβˆ1,βˆ1, respectively. These estimates are plugged in (7) in order to get an estimate σˆp2 and a (1α) level confidence interval for p may be obtained as

(8)pˆ±z1α/2σˆpˆ

with z1α/2 being the appropriate quantile of the normal distribution.

Past authors concentrated on methods for testing that the mediation proportion is at least some fraction f, with f typically being 0.5 or more [14, 15, 40]. In the context of PTE, where the validation of intermediate biomarkers for outcome is of interest, this may be reasonable. However, when considering a mediator, the more relevant question is whether M is indeed a mediator. Then, the hypothesis is H0:p=0 vs. H1:p0. Let Zp=σp1pˆ be the scaled estimate. By the delta method, the distribution of Zp under the null converges to a standard normal distribution and the null is rejected at significance level α if |Zp|>z1α/2.

An alternative test statistic is based upon a test for the difference between the effect estimates in the marginal and conditional models. That is, on dˆ=βˆ1βˆ1. Under the assumptions in this paper, dˆ is a consistent estimate for the NIE. A test statistic based on dˆ is based on Zd=σdˆ1dˆ, where

σdˆ2=Var(dˆ)=σβˆ12+σβˆ122σβˆ1,βˆ1.

Then, the null hypothesis is rejected if |Zd|>z1α/2.

4.1 The data duplication algorithm

A main challenge when conducting inference for the mediation proportion p is to estimate the covariance of estimators obtained from the two models (2) and (3). It turns out that the covariance between Bˆ and Bˆ can be estimated by fitting both models by stacking the estimating equations for the two models using a data duplication algorithm. A similar method was presented in [15] for the Cox model in survival data. Here, we extend it to GEE for GLMs. First, the data are augmented with additional pseudo-variables and pseudo-observations. Each variable, including the intercept, the exposure, the confounders, but not the mediator, appears twice, and each of the original observations is included as two pseudo-observations in the new data set. See Table 1 for an illustration of the duplicated data structure.

Table 1:

The augmented data used by the data duplication algorithm. For each original observation i, two rows j=1,2 are created. The duplicated data is used for the pseudo model presented in eq. (9).

ijInterceptInterceptXXMWWY&% Model
1110x10m1w10y1
12010x100w1y1&%2
2110x20m2w20y2
22010x200w2y2
&%

The following pseudo model is fitted to the duplicated data using GEE [27],

(9)E(Yij|Xi,Xi,Mi,Wi,Wi)=g1(β0I{j=1}+β1Xi+β2Mi+β3TWi+β0I{j=2}+β1Xi+β3TWi),

where j=1,2 are the rows created from duplicating each observation and are treated as repeated measures. Model (9) implies that we can write E(Yi1|Xi,Xi,Mi,Wi,Wi)=E(Yi1|Xi,Mi,Wi) and E(Yi2|Xi,Xi,Mi,Wi,Wi)=E(Yi2|Xi,Wi). Let R be a 2×2 working correlation matrix and denote Bi=diag(vi1,vi2), where vij=Var(Yij). Let also Vi=Bi1/2RBi1/2 be a 2×2 working variance for the vector (Yi1,Yi2). Here, the GEE are defined as

(10)UGEE(B)=i=1n(Di,Di)Vi1yi1E(Yi1|Xi,Mi,Wi)yi2E(Yi2|Xi,Wi),

where Di=E(Yi1|Xi,Mi,Wi)/B and Di=E(Yi2|Xi,Wi)/B are two column vectors. If R is taken to be the identity matrix, then Vi=Bi and eq. (10) simplifies to the following estimating equations

(11)UIEE(B)=UIEE(1)(β)UIEE(2)(β)=i=1nDivi11[yiE(Yi1|Xi,Mi,Wi)]i=1nDivi21[yiE(Yi2|Xi,Wi)]=0.

Then, the estimating equations given by eq. (11) are identical to the estimating equations for fitting models (2) and (3) separately, because Di, vi and Zi in eq. (1) are equal to Di,vi1 and (Xi,Mi,Wi), respectively, under model (2), and they are equal to Di,vi2 and (Xi,Wi), respectively, under model (3). The major advantage of the data duplication algorithm is that it provides an estimator for σβ1,β1 in a straightforward manner. Taking a working correlation matrix other than the identity may result in more efficient estimators of Bˆ, but would not have the desirable property that the duplicated data estimating equations are identical to the two separate estimating equations from the two separate models.

5 Simulation study

The simulation studies and data analysis were conducted in R. The code is available upon request from the first author. In addition, we have developed a SAS macro and an R package that are publicly available (Appendix A.1). In the simulation studies, we considered several issues regarding the performance of the methodology we presented throughout the paper. We first present results concerning g-linkability for the logit link function and the Cox model. Then, we turn to the performance of the mediation proportion estimator, studying its bias, the coverage rate of the accompanied confidence interval and the type I error and the power of the statistical tests described in Section 4.

Throughout these simulation studies, we assume that there are no confounders in the model. X and M were generated using a bivariate normal with mean (0,0)T and covariance matrix 1ρρ1. Then, we have that β2=pρβ1 for the identity, log and logit link functions (the latter under the rare outcome assumption); see Web Appendix A. In these scenarios, glinkability holds for all three link functions. The estimation and inference procedures apply to any bivariate distribution of X and M that satisfies the simple moment conditions given in Section 3, and here we used the bivariate normal distribution for generating the data merely for convenience. The estimation and inference procedures do not use the bivariate normal distribution of (X,M).

5.1 g-linkability for the logit link function and of the Cox model

In order to assess the magnitude of the bias when assuming g-linkability of the logit link function and the Cox model, we conducted a simulation study under various conditions and inspected the resulting bias in pˆ, as estimated using the data duplication algorithm described in Section 4.1 while taking the working correlation matrix to be the identity. First we describe the logit link function model. We simulate Y under the logistic regression model

logit(P(Y=1|X,M))=β0+β1X+β2M.

We chose the model parameter values in the following way. First, we chose ρ=corr(X,M),p and β1. Then, by definition we had β1=(1p)β1, and we took β2 as if g-linkability exactly holds. That is, β2=pρβ1. Then, we fixed the unconditional case probability P(Y=1) and found the appropriate β0 value by solving for β0 in the equation

P(Y=1)=E(expit(β0+β1X+β2M)),

where expit(u)=exp(u)/(1+exp(u)). Finally, the sample size was given as n=E(Ncases)/P(Y=1) where E(Ncases) is number of expected cases. We considered the following values for the parameters. p=0.1,0.2,...,0.8; ρ=p,p+0.1,...,0.8, with ρp to satisfy that β2β1 or in words, to ensure that the total effect of X is larger than effect of M; β1=log(1.25),log(1.5),log(2); P(Y=1)=0.005,0.01,0.1,0.25; E(Ncases)=100,500,1000. The number of simulation iterations per scenario was 1000.

For the Cox model, we simulated the data similarly to the logit link function simulations. First, we simulated X and M as before. Then, given fixed ρ,p and β1, β2=pρβ1. We took a Weibull distribution for the baseline hazard and used Exponential distribution for the censoring (mean=50), with additional cutoff at age 90. Given the desired proportion number of cases in the population, we used simulations to find the appropriate values for the Weibull distribution shape parameter, while fixing the scale parameter at 200. As in the logit link case, we chose the sample size as the number of expected cases (E(Ncases)) divided by the expected proportion of cases (P(δ=1)), where δ is the event indicator.

In order to assess g-linkability, and the finite sample performance of pˆ, we calculated the relative bias, defined as 100×|mean(pˆ)pp|. Ideally, this quantity should be close to zero. We note that bias may arise either because g-linkability fails to hold, or because of a sample size not large enough. Figure 1 presents bias for β1=log(1.5) as a function of the parameters. First, it is of note that whenever the overall prevalence or cumulative incidence of Y was small, as in the rare disease scenario, and the number of cases was sufficiently large, bias was minimal. Even when the disease was not as rare, e.g., P(Y=1)=0.25, when there were enough cases, and when p was large enough (e.g., p>0.2 in this case), the bias was minimal.

Figure 1: Relative bias of the mediation proportion estimator under the logistic model as a function of the mediation proportion (p)$(p)$, the correlation between the exposure and the mediator (ρ)$(\rho)$, the number of expected cases (E(Ncases))$(E(N_{cases}))$ and the outcome rate (P(Y=1))$(P(Y=1))$. The value of β1⋆$\beta^{\star}_{1}$ was taken to be log(1.5)$\log(1.5)$.
Figure 1:

Relative bias of the mediation proportion estimator under the logistic model as a function of the mediation proportion (p), the correlation between the exposure and the mediator (ρ), the number of expected cases (E(Ncases)) and the outcome rate (P(Y=1)). The value of β1 was taken to be log(1.5).

Figure 2: Relative bias of the mediation proportion estimator under the Cox model as a function of the mediation proportion (p)$(p)$, the correlation between the exposure and the mediator (ρ)$(\rho)$, the number of expected cases (E(Ncases))$(E(N_{cases}))$ and the event rate (P(δ=1))$(P(\delta=1))$. The value of β1⋆$\beta^{\star}_{1}$ was taken to be log(1.5)$\log(1.5)$.
Figure 2:

Relative bias of the mediation proportion estimator under the Cox model as a function of the mediation proportion (p), the correlation between the exposure and the mediator (ρ), the number of expected cases (E(Ncases)) and the event rate (P(δ=1)). The value of β1 was taken to be log(1.5).

Considering the g-linkability of the Cox model, presented for β1=log(1.5) in Figure 2, the results were similar to the results obtained for the logit link function. That is, when the outcome was rare (P(δ=1) was small) then the corresponding marginal Cox model for the hazard function approximately holds and the bias in mediation proportion estimation was minimal. Figures similar to Figures 1 and 2 are presented for β1=log(1.25),log(2) in Web Appendix B. The overall trends were similar.

5.2 Estimation and inference performance

For the Cox model and the logit link function, data were simulated as described above. For the identity link function, data were simulated from the model Y=E(Y|X,M)+ϵ=β0+β1X+β2M+ϵ. As before, (X,M) were simulated from a bivariate normal distribution with zero mean, unit variance, and correlation ρ. The error, ϵ, was a vector of iid standard normal random variables. We also considered other distributions for ϵ; we will expand on this matter later on. As in the logit link case, we fixed β1, ρ, and p and took β1=(1p)β1 and β2=pρβ1. The intercept, β0, was chosen arbitrarily to be equal to 2. We considered various values for β1,p and ρ, where as before we were only interested in scenarios where ρp, since then β2β1. We present results for β1=0.1,0.3,0.5, which imply multiple correlations between (X,M) and Y of about 0.1,0.3 and 0.5, respectively.

Under the simple linear model (identity link), estimates obtained from the difference and product methods are algebraically identical [41]. This result does not extend to the logistic link function or to the Cox model [3, 42]. We therefore compared between the difference and the product methods. Let pˆdiff and pˆprod be the difference and product method estimators, respectively. Table 2 presents %RBiasDiff=100×[|mean(pˆdiff)pp||mean(pˆprod)pp|], and the empirical variance ratio Var(pˆdiff)Var(pˆratio). The two methods are generally comparable when g-linkability holds.

Table 2:

Percent relative bias and efficiency of estimators of the mediation proportion under the logistic link function. RBiasDiff is the difference between the relative bias of difference method compared to that of the product method times 100. Negative values indicate that the difference method is preferable. Vratio is the ratio of the empirical variances of the estimates of the difference and product methods. Vratio>1 indicates that the product method is preferable.

P(Y=1)pρexp(β1)=1.25exp(β1)=1.5exp(β1)=2
%RBiasDiffVratio%RBiasDiffVratio%RBiasDiffVratio
0.0050.10.1-0.231.00-0.821.002.841.14
0.50.011.000.041.00-0.141.00
0.7-0.011.000.021.00-0.041.00
0.30.3-0.071.00-0.201.000.120.99
0.5-0.021.00-0.071.00-0.260.99
0.7-0.001.00-0.021.00-0.101.00
0.50.5-0.021.00-0.081.00-0.331.00
0.7-0.011.00-0.021.00-0.111.00
0.010.10.1-0.550.99-1.001.026.001.13
0.50.021.00-0.081.000.220.99
0.70.021.000.041.00-0.121.00
0.30.3-0.111.00-0.420.990.881.01
0.5-0.031.00-0.121.00-0.040.99
0.7-0.011.00-0.041.000.151.00
0.50.5-0.041.00-0.141.00-0.351.00
0.7-0.011.00-0.041.00-0.201.00
0.10.10.1-1.440.9711.520.9746.421.33
0.50.271.00-0.630.991.450.96
0.7-0.181.00-0.321.00-0.240.98
0.30.3-0.930.99-0.860.9711.260.97
0.5-0.301.00-1.050.982.990.94
0.7-0.101.00-0.360.99-1.200.97
0.50.5-0.301.00-0.351.002.711.01
0.7-0.091.00-0.411.000.620.99
  1. %Rbiasdiff=100×Mean(pˆDiff)ppMean(pˆProd)pp

Table 3:

Type I error and power for tests for mediation under the identity (n=1000) and logit link functions and the Cox model (E(Ncases)=500).

Identity link function (YNormal)
β1=0.1β1=0.3β1=0.5
pˆtestdˆtestpˆtestdˆtestpˆtestdˆtest
p=0.0ρ=0.10.000.020.010.010.010.02
ρ=0.50.010.040.050.050.050.05
ρ=0.70.010.050.040.050.060.06
p=0.1ρ=0.10.280.600.910.900.900.88
ρ=0.50.020.070.370.380.790.79
ρ=0.70.020.050.140.150.350.36
p=0.2ρ=0.30.280.521.001.001.001.00
ρ=0.50.070.180.900.911.001.00
ρ=0.70.030.100.470.490.890.89
p=0.3ρ=0.30.540.831.001.001.001.00
ρ=0.50.170.361.001.001.001.00
ρ=0.70.060.180.830.840.991.00
Logit link function (YBer) with P(Y=1)=0.01
β1=log(1.25)β1=log(1.5)β1=log(2)
pˆtestdˆtestpˆtestdˆtestpˆtestdˆtest
p=0.0ρ=0.10.030.060.040.050.040.04
ρ=0.50.030.050.050.060.050.05
ρ=0.70.040.050.040.040.050.05
p=0.1ρ=0.11.001.001.001.001.001.00
ρ=0.50.090.130.340.360.740.74
ρ=0.70.060.090.150.160.350.36
p=0.2ρ=0.30.850.891.001.001.001.00
ρ=0.50.310.400.860.871.001.00
ρ=0.70.130.170.420.440.870.87
p=0.3ρ=0.30.991.001.001.001.001.00
ρ=0.50.660.731.001.001.001.00
ρ=0.70.250.320.740.760.990.99
Cox model with P(δ=1)=0.01
β1=log(1.25)β1=log(1.5)β1=log(2)
pˆtestdˆtestpˆtestdˆtestpˆtestdˆtest
p=0.0ρ=0.10.030.050.040.040.040.04
ρ=0.50.030.050.040.050.060.06
ρ=0.70.030.050.060.060.060.06
p=0.1ρ=0.10.991.001.001.001.001.00
ρ=0.50.090.140.330.350.750.76
ρ=0.70.050.090.150.160.320.33
p=0.2ρ=0.30.840.881.001.001.001.00
ρ=0.50.330.430.850.861.001.00
ρ=0.70.150.190.460.470.870.88
p=0.3ρ=0.30.991.001.001.001.001.00
ρ=0.50.670.751.001.001.001.00
ρ=0.70.250.330.780.801.001.00

From estimation we move to hypothesis testing. The two test statistics compared were described in Section 4, where the variance estimators used in the test statistics were obtained by the data duplication algorithm described in Section 4.1. Results are presented in Table 3. In terms of type I error, both tests were adequate, with a conservative type I error when the correlation between the exposure and the mediator was low. When the total effect was low, the test based on dˆ had greater power, usually by 5%10%, compared to the test based on pˆ. The power of both tests was highly affected by the effect size (β1) and the correlation between the exposure and the mediator (ρ). High correlation between X and M decreased the power. The power of both tests was lower when the total effect β1 is low. It should be noted that mediation analysis is performed only after a risk factor was found to be significant, which, in general, is less likely to happen if both β1 and the sample size are low. We further address this point in Section 7. We next consider the finite sample properties of the confidence interval presented in Section 4. Table 4 presents empirical coverage rates and confidence interval widths. Coverage rates were generally adequate. When ρ and p were of similar size, and the total effect was large, the new method did not produce confidence intervals with nominal coverage. For the logistic link function and the Cox model, the worst results were obtained for the combination of p=ρ=0.1 and log(β1)=2. As can be seen from Figures 1 and 2 (and from the figures in Web Appendix B) in these cases, g-linkability does not hold because the outcome is not that rare and the relative risk is strong.

Table 4:

Empirical coverage rates (CR) and lengths (LEN) of confidence intervals for the mediation proportion under the identity and logit link functions and the Cox model

Identity link function (YNormal)
n=1000n=5000
β1=0.1β1=0.3β1=0.5β1=0.1β1=0.3β1=0.5
p=0.1ρ=0.1CR0.900.960.950.940.950.94
LEN0.240.130.120.100.060.05
ρ=0.5CR0.980.950.960.960.950.96
LEN0.850.250.150.330.110.07
ρ=0.7CR0.980.960.950.960.950.96
LEN1.420.410.240.560.180.11
p=0.3ρ=0.3CP0.930.950.940.960.970.95
LEN0.650.200.140.250.090.06
ρ=0.5CR0.960.960.940.960.950.96
LEN0.930.280.170.370.120.07
ρ=0.7CR0.970.960.960.950.950.97
LEN1.510.430.260.590.190.11
p=0.5ρ=0.5CR0.940.940.960.960.960.95
LEN1.110.320.200.430.140.09
ρ=0.7CR0.960.950.950.960.960.95
LEN1.650.470.270.630.200.12
Logit link function (YBer) with P(Y=1)=0.01
E(Ncases)=500E(Ncases)=1000
β1β1
log(1.25)log(1.5)log(2)log(1.25)log(1.5)log(2)
p=0.1ρ=0.1CR0.950.950.870.960.930.80
LEN0.130.070.040.080.050.03
ρ=0.5CR0.970.960.940.960.950.96
LEN0.490.260.150.330.190.11
ρ=0.7CR0.960.950.960.950.940.95
LEN0.840.430.250.570.310.18
p=0.3ρ=0.3CR0.960.950.920.950.950.94
LEN0.390.200.110.260.130.08
ρ=0.5CR0.960.950.930.960.960.95
LEN0.550.290.170.370.200.12
ρ=0.7CR0.960.940.960.950.950.94
LEN0.880.470.260.590.320.19
p=0.5ρ=0.5CR0.950.950.940.960.960.94
LEN0.650.340.200.450.240.14
ρ=0.7CR0.970.940.950.960.950.93
LEN0.950.490.290.660.340.20
Cox model with P(δ=1)=0.01
E(Ncases)=500E(Ncases)=1000
β1β1
log(1.25)log(1.5)log(2)log(1.25)log(1.5)log(2)
p=0.1ρ=0.1CR0.930.950.920.930.810.73
LEN0.120.080.070.050.060.04
ρ=0.5CR0.960.950.960.940.940.95
LEN0.480.340.260.180.150.10
ρ=0.7CR0.980.950.950.960.970.94
LEN0.830.560.430.300.250.18
p=0.3ρ=0.3CR0.950.960.950.950.930.91
LEN0.380.260.190.130.110.08
ρ=0.5CR0.970.960.950.950.950.94
LEN0.570.380.290.200.170.12
ρ=0.7CR0.960.940.950.960.950.95
LEN0.870.600.450.320.260.18
p=0.5ρ=0.5CR0.960.950.940.960.940.94
LEN0.660.440.340.240.200.14
ρ=0.7CR0.960.960.950.950.950.95
LEN0.940.630.490.340.280.20
Table 5:

Ratio between mean estimated Covˆ(βˆ1,β1ˆ) and empirical Cov(βˆ1,β1ˆ) under the identity (n=1000) and logit link functions and the Cox model E(Ncases)=100

Identity link function (YNormal)
β1=0.1β1=0.3β1=0.5
p=0.1ρ=0.10.9890.9940.985
ρ=0.51.0361.0140.963
ρ=0.70.9600.9621.098
p=0.3ρ=0.31.0310.9590.967
ρ=0.51.1030.8890.989
ρ=0.70.9880.9750.955
p=0.5ρ=0.51.0640.9591.008
ρ=0.71.0100.9801.100
Logit link function (YBer) with P(Y=1)=0.1
β1=log(1.25)β1=log(1.5)β1=log(2)
p=0.1ρ=0.10.9960.9681.031
ρ=0.51.0700.9771.023
ρ=0.71.0120.9760.940
p=0.3ρ=0.31.0020.9541.057
ρ=0.50.9541.0161.080
ρ=0.70.9580.9020.995
p=0.5ρ=0.50.9930.9831.008
ρ=0.50.9961.0110.959
Cox model with P(δ=1)=0.1
β1=log(1.25)β1=log(1.5)β1=log(2)
p=0.1ρ=0.11.0891.0350.951
ρ=0.51.0111.0130.945
ρ=0.71.0991.0441.068
p=0.3ρ=0.31.0530.8920.906
ρ=0.50.9570.9900.887
ρ=0.70.9990.9820.966
p=0.5ρ=0.50.9450.9381.036
ρ=0.50.9611.0570.963

For the difference method, we also compared the confidence intervals using the asymptotic variance to confidence intervals constructed using the bootstrap, both by estimating the bootstrapped variance and assuming normality, and by using the quantiles of the bootstrap samples. The results, presented in Web Appendix C, show that the asymptotic approach is comparable to both bootstrap procedures in terms of nominal coverage. Both versions of the bootstrap confidence intervals were wider than the asymptotic confidence intervals. Furthermore, the bootstrap is time consuming, especially for large data sets, where our method, implemented by publicly-available software, is as fast as a single GEE model fit. For n=5000, the computation time of 500 bootstrap replications, for a single data set, was about 115 seconds; For n=10000, it was 250 seconds, and for n=50000 it was about 1300 seconds (more than 20 minutes). For comparison, calculation of confidence intervals in our method was less than one second for n=5000,10000 and about 4 seconds for n=50000.

Throughout this section, we presented in parallel results for the identity and logit link function and the Cox model. There was a very strong agreement between the results for the logit link function for binary data and the Cox model, as one may have expect given the close relationship between the logistic regression model and the Cox model in epidemiology and public health evaluations.

In addition to the scenarios we described above, we conducted simulations for the identity link function with error distributions other than the normal one. We considered symmetric distribution with tails heavier than the normal distribution as well as skewed distributions. As predicted by GEE theory, the performance of the mediation proportion estimator, the statistical tests and the confidence interval was only slightly changed. Details are given in Web Appendix D.

6 Illustrative example

We illustrate the use of our methodology in an analysis of the etiology of pre-menopausal breast cancer data from the Nurses Health’s Studies (NHS and NHSII) [43, 44]. It was previously found that high mammographic density (MD) is a risk factor for breast cancer [45]. The goal here is to investigate whether, and to what extent, the effects of more distal risk factors for pre-menopausal breast cancer are mediated by high MD. A detailed description of this study is given in [46]. In this nested case-control study, controls were matched to cases by current age, menopausal status, current hormone use, month, time of day, fasting status and time of the day at blood collection and luteal day (for NHSII samples only). There were 559 pre-menopausal cases and 1727 controls. Since the disease is rare, as shown in the previous section, g-linkability should hold. Mediation by percent MD, a single mediator, was considered for a number of well-established breast cancer risk factors. Following [46], the mediation analysis was conducted for each risk factor separately.

We only considered risk factors with significant total effects: personal history of benign breast disease (HBBD), family history of breast cancer (FH), adolescent somatotype (ASM), body mass index at age 18 (BMI18), age at first birth (AFB), age at menarche (AM) and height (HT). In each of the analyses, we included potential confounders for the risk factor-MD, MD-outcome and risk factor-outcome relationships, based upon subject matter considerations following [46]. For example, we did not adjust for current (adult) BMI in the analysis of BMI18, as the latter affects the former. In addition, some of the risk factors studied may have been confounders in analysis of mediation via MD of another risk factor. The set of confounders used in at least one analysis included current age, fasting status, blood collection time of the day, mammography batch (NHS batch 1, NHS batch 2 or NHSII), current BMI, BMI18, ASM, HBBD, parity, AFB, and AM. As in most observational studies, residual confounding may bias our results.

Since our method assumes no exposure-mediator interaction, we fitted, for each risk factor, a logistic regression that included the risk factor, mediator, potential confounders and an interaction term involving the risk factor and the mediator. Then we tested for interaction using a standard Wald test.

Table 6 presents the estimated mediation proportions, confidence intervals, p-values, the estimated risk factor effects, and the p-value for risk factor-mediator interaction. MD was a significant mediator (at the 5% significance level) for HBBD, ASM and BMI18, regardless whether the test was based on pˆ or dˆ, although p-values corresponding to the latter test were much smaller. MD was significant as a mediator for AM according to the dˆ test but not according to the less powerful pˆ test. Confidence intervals were quite wide for the mediation proportions for ASM and BMI18. This may be due to the moderate sample size, and the relatively small effect.

There was no evidence for risk factor-mediator interaction for all but one risk factor studied here. For HT, the test for the interaction term was significant. However, the point estimate for the proportion of the effect of HT mediated through MD is very close to zero. Thus, the fact that this assumption is violated is unlikely to be of substantive importance. In the supplementary materials of [46], it was reported that when taking the interaction into account using the method of [30], the mediation proportion remained very small, although positive.

The results suggest that, if the non-confounding assumptions needed for causal interpretation of observed associations are met, MD mediates the effect of at least some pre-menopausal breast cancer risk factors, with evidence for a large mediation proportion for BMI18 and ASM and some mediation of HBBD, but not for the other risk factors.

Table 6:

Mediation analysis for pre-menopausal breast cancer incidence with mammographic density as the mediator in the NHS and NHSII studies (N= 559 cases and 1727 controls). RRˆtotal=exp(βˆ), RRˆdirect=exp(βˆ).

Risk factorpinterβˆ(RRˆtotal)p-valueβˆ(RRˆdirect)pˆ95% CIp-value, pˆ testp-value, dˆ test
Personal history of benign breast disease0.950.35 (1.42)0.0010.25 (1.28)0.300.10-0.510.004<106
Family history of breast cancer0.590.42 (1.52)0.010.42 (1.52)0.004-0.10-0.110.940.94
Adolescent somatotype Per 3 unit increase0.20-0.34 (0.72)0.02-0.12 (0.88)0.630.05-1.200.03<107
BMI at age 18 Per 5 unit increase0.20-0.23 (0.79)0.02-0.05 (0.95)0.780.06-1.500.03<107
Age at first birth Per 5 year increase0.170.15 (1.17)0.030.15 (1.16)0.03-0.09-0.150.310.30
Age at menarche Per 2 year increase0.22-0.16 (0.86)0.03-0.18 (0.84)-0.16-0.36-0.040.120.04
Height Per 3 inch increase0.030.13 (1.14)0.030.14 (1.14)-0.01-0.14-0.110.820.82
  1. Adjusted for age, fasting status, blood collection time of the day, mammography batch (NHS batch 1, NHS batch 2 or NHSII), current and at age 18 BMI, adolescent somatotype, history of BBD, parity, age at first birth, and age at menarche P-value of the test for interaction between percent MD and each risk factor Not adjusted for adolescent somatotype, BMI, current or at age 18 Among parous women only (478 cases, 1499 controls)

7 Discussion

In this paper, we have provided methodology for estimation and inference for the mediation proportion in GLMs and the Cox model using the difference method. Our methodology for GLMs uses a data duplication algorithm with GEE and allows for the consistent estimation of the covariance of the estimates.

Strictly speaking, the validity of the difference method relies on the assumption that the marginal model, the one that does not include the mediator, and the conditional model, the one that does, hold simultaneously. However, we demonstrate in this paper that g-linkability with respect to the mean functions ensures that the point estimator for the mediation proportion is consistent under standard assumptions for the identity and log link functions and under a rare outcome assumption for the logistic link function and Cox model. The rare outcome assumption is fulfilled in most chronic disease incidence studies in epidemiology, including the one motivating the present work. When the outcome is not rare, one may fit the log-binomial model instead, as noted in [38], which may be preferable anyway, as the odds ratio is typically not the parameter of interest [47]. Furthermore, the estimator is asymptotically normally distributed with a variance that can be consistently estimated using a robust sandwich estimator easily obtained by applying a data duplicated GEE. In some scenarios, g-linkability fails to hold, even approximately. A direction for future research is alternative definitions and estimation procedures that are based on nonparametric projections instead of exact generalized linear models [48].

Despite its popularity, the difference method for estimating the mediation proportion has been criticized due to what appeared to be undesirable finite samples properties [4, 40]. However, when considering binary outcomes, the covariance (or correlation) between estimates from the marginal and conditional models was typically estimated using approximations from the linear model [1]. We have now developed methodology for a valid covariance estimator and showed that testing for mediation using a test based on the difference yields a valid statistical test, even in finite samples.

An alternative to the difference method is the product method. In terms of finite sample properties, we have shown the two methods are comparable under g-linkability. A major advantage of using the difference method is that variance can be directly estimated by standard software, without relying on approximations or on the bootstrap, which are typically used for inference when working with the product method. Thus, we have provided a valid, simple alternative and provided software in SAS and R.

The causal structure and the underlying confounding assumptions are important to consider when our methods are used in applications. Confounding may occur due to exposure-mediator confounders, exposure-outcome confounders or mediator-outcome confounders. We refer the readers to [19], and references therein, for relevant discussions on assumptions needed and analysis conducted in order to avoid, or at least minimize, potential bias due to confounding when conducting mediation analysis. The difference method does not allow for mediator-exposure interaction, and alternative methods to allow for this interaction were previously developed [34, 38, 49].

In practice, mediation analysis is often conducted for well-established exposures or risk factors, or when the total effect is significant. As suggested by our simulation results, when the total effect was small, mediation analysis was less likely to provide adequate results. On the other hand, an analysis that only considers significant total effects should take into account that it was performed conditionally on the results of a first stage analysis. The properties of such conditional inference can be considered in future research.

In our implementation of the GEE methodology, we propose to use the independence working correlation matrix, which has the nice property of providing identical coefficient estimates when fitting the two models separately and when using the data duplication algorithm, fitting them together. Under other working correlation matrices, this property does not hold anymore, but efficiency may be gained.

In conclusion, the general framework for mediation analysis in GLMs developed in this paper along with the methodology established, will allow researchers to investigate mediation under various outcome scenarios and to quantify results based on rigorously derived and empirically studied estimators and hypothesis tests.

A Appendix

A.1 Software

The SAS macro %mediate implements the data duplication algorithm and reports point and interval estimates for the mediation proportion and the results for the mediation test using the difference method. It is available on the last author’s website http://www.hsph.harvard.edu/donna-spiegelman/software/mediate. The SAS macro supports GLMs and the Cox model. We also developed an R package named GEEmediate which is available on CRAN. The current version of the R package can be used for GLMs.

References

[1] Freedman LS, Schatzkin A. Sample size for studying intermediate endpoints within intervention trials or observational studies. Am. J. Epidemiol. 1992;136:1148–1159.10.1093/oxfordjournals.aje.a116581Search in Google Scholar PubMed

[2] MacKinnon DP. Introduction to statistical mediation analysis Routledge, 2008.Search in Google Scholar

[3] MacKinnon DP, Dwyer JH. Estimating mediated effects in prevention studies. Eval. Rev. 1993;17:144–158.10.1177/0193841X9301700202Search in Google Scholar

[4] MacKinnon DP, Lockwood CM, Hoffman JM, West SG, Sheets V. A comparison of methods to test mediation and other intervening variable effects. Psychol. Meth. 2002;7:83.10.1037/1082-989X.7.1.83Search in Google Scholar

[5] Baron RM, Kenny DA. The moderator–mediator variable distinction in social psychological research: Conceptual, strategic, and statistical considerations. J. Personality Social Psychol. 1986;51:1173.10.1037/0022-3514.51.6.1173Search in Google Scholar

[6] Imai K, Keele L, Identification Yamamoto T., inference and sensitivity analysis for causal mediation effects. Stat. Sci. 2010b:51–71.10.1214/10-STS321Search in Google Scholar

[7] Pearl J. Direct and indirect effects. In: Proceedings of the seventeenth conference on uncertainty in artificial intelligence. Morgan Kaufmann Publishers Inc., 2001:411–420.10.1145/3501714.3501736Search in Google Scholar

[8] Robins JM, Greenland S. Identifiability and exchangeability for direct and indirect effects. Epidemiology 1992:143–155.10.1097/00001648-199203000-00013Search in Google Scholar PubMed

[9] Huang B, Sivaganesan S, Succop P, Goodman E. Statistical assessment of mediational effects for logistic mediational models. Stat. Med. 2004;23:2713–2728.10.1002/sim.1847Search in Google Scholar PubMed

[10] Wang W., Albert JM. Estimation of mediation effects for zero-inflated regression models. Stat. Med. 2012;31:3118–3132.10.1002/sim.5380Search in Google Scholar PubMed PubMed Central

[11] Huang Y-T, Pan W-C. Hypothesis test of mediation effect in causal mediation model with high-dimensional continuous mediators. Biometrics 2015.10.1111/biom.12421Search in Google Scholar PubMed

[12] Alwin DF, Hauser RM. The decomposition of effects in path analysis. Am. Sociological Rev. 1975:37–47.10.2307/2094445Search in Google Scholar

[13] Judd CM. Kenny DA. Process analysis estimating mediation in treatment evaluations. Eval. Rev. 1981;5:602–619.10.1177/0193841X8100500502Search in Google Scholar

[14] Freedman LS, Graubard BI, Schatzkin A. Statistical validation of intermediate endpoints for chronic diseases. Stat. Med. 1992;11:167–178.10.1002/sim.4780110204Search in Google Scholar PubMed

[15] Lin D, Fleming T, De Gruttola V. Estimating the proportion of treatment effect explained by a surrogate marker. Stat. Med. 1997;16:1515–1527.10.1002/(SICI)1097-0258(19970715)16:13<1515::AID-SIM572>3.0.CO;2-1Search in Google Scholar PubMed

[16] Parast L, McDermott MM, Tian L. Robust estimation of the proportion of treatment effect explained by surrogate marker information. Stat Med. 2015.10.1002/sim.6820Search in Google Scholar PubMed

[17] Sobel ME. Asymptotic confidence intervals for indirect effects in structural equation models. Soc. Method. 1982;13:290–312.10.2307/270723Search in Google Scholar

[18] Goodman LA. On the exact variance of products. J. Am. Stat. Assoc. 1960;55:708–713.10.1080/01621459.1960.10483369Search in Google Scholar

[19] VanderWeele T. Explanation in causal inference: methods for mediation and interaction. Oxford University Press, 2015.10.1093/ije/dyw277Search in Google Scholar

[20] Imai K, Keele L, Tingley D. A general approach to causal mediation analysis. Psychol. Meth. 2010a;15:309.10.1037/a0020761Search in Google Scholar

[21] Tingley D, Yamamoto T, Hirose K, Keele L, Imai K. Mediation: R package for causal mediation analysis. 2014.10.18637/jss.v059.i05Search in Google Scholar

[22] Carney, RM, Howells WB, Blumenthal JA, Freedland KE, Stein PK, Berkman LF, Watkins LL, Czajkowski SM, Steinmeyer B, Hayano J, et al. Heart rate turbulence, depression, and survival after acute myocardial infarction. Psychosomatic Med. 2007;69:4–9.10.1097/01.psy.0000249733.33811.00Search in Google Scholar PubMed

[23] Lyall K, Ashwood P, Van de Water J, Hertz-Picciotto I. Maternal immune-mediated conditions, autism spectrum disorders, and developmental delay. J. Autism Dev. Disorders 2014;44, 1546–1555.10.1007/s10803-013-2017-2Search in Google Scholar PubMed PubMed Central

[24] Reisner SL, Greytak EA, Parsons JT, Ybarra ML. Gender minority social stress in adolescence: disparities in adolescent bullying and substance use by gender identity. J. Sex Res. 2015;52:243–256.10.1080/00224499.2014.886321Search in Google Scholar PubMed PubMed Central

[25] Roberts AL, Rosario M, Corliss HL, Koenen KC, Austin SB. Childhood gender nonconformity: A risk indicator for childhood abuse and posttraumatic stress in youth. Pediatrics 2012;129:410–417.10.1542/peds.2011-1804Search in Google Scholar PubMed PubMed Central

[26] Huber PJ. Robust statistics Springer, 2011.10.1007/978-3-642-04898-2_594Search in Google Scholar

[27] Liang K-Y, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika 1986:13–22.10.1093/biomet/73.1.13Search in Google Scholar

[28] Pearl J. Causal inference in statistics: An overview. Stat. Surv. 2009;3:96–146.10.1214/09-SS057Search in Google Scholar

[29] Pearl J. The causal mediation formula–a guide to the assessment of pathways and mechanisms. Prev. Sci. 2012;13:426–436.10.1007/s11121-011-0270-1Search in Google Scholar PubMed

[30] VanderWeele TJ, Vansteelandt S. Odds ratios for mediation analysis for a dichotomous outcome. Am. J. Epidemiol. 2010;172:1339–1348.10.1093/aje/kwq332Search in Google Scholar PubMed PubMed Central

[31] Valeri, L., Lin X, VanderWeele TJ. Mediation analysis when a continuous mediator is measured with error and the outcome follows a generalized linear model. Stat. Med. 2014;33:4875–4890.10.1002/sim.6295Search in Google Scholar PubMed PubMed Central

[32] Tchetgen Tchetgen EJ. Inverse odds ratio-weighted estimation for causal mediation analysis. Stat. Med. 2013;32:4567–4580.10.1002/sim.5864Search in Google Scholar PubMed

[33] Lange T, Hansen JV. Direct and indirect effects in a survival context. Epidemiology 2011;22:575–581.10.1097/EDE.0b013e31821c680cSearch in Google Scholar

[34] Lange T, Vansteelandt S, Bekaert M. A simple unified approach for estimating natural direct and indirect effects. Am. J. Epidemiol. 2012;176:190–195.10.1093/aje/kwr525Search in Google Scholar PubMed

[35] Tchetgen Tchetgen EJ. On causal mediation analysis with a survival outcome. Int. J. Biostat. 2011;7:1–38.10.2202/1557-4679.1351Search in Google Scholar

[36] VanderWeele TJ. Causal mediation analysis with survival data. Epidemiology (Cambridge, Mass.) 2011;22:582.10.1097/EDE.0b013e31821db37eSearch in Google Scholar PubMed

[37] Ritz J, Spiegelman D. Equivalence of conditional and marginal regression models for clustered and longitudinal data. Stat. Meth. Med. Res. 2004;13:309–323.10.1191/0962280204sm368raSearch in Google Scholar

[38] 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:137.10.1037/a0031034Search in Google Scholar PubMed

[39] Jiang Z, VanderWeele TJ. When is the difference method conservative for assessing mediation? Am. J. Epidemiol. 2015;182:105–8.10.1093/aje/kwv059Search in Google Scholar

[40] Freedman LS. Confidence intervals and statistical power of the “validation” ratio for surrogate or intermediate endpoints. J. Stat. Plann. Inference 2001;96:143–153.10.1016/S0378-3758(00)00330-XSearch in Google Scholar

[41] MacKinnon DP, Warsi G, Dwyer JH. A simulation study of mediated effect measures. Multivariate Behav. Res. 1995;30, 41–62.10.1207/s15327906mbr3001_3Search in Google Scholar PubMed PubMed Central

[42] Tein J-Y, MacKinnon DP. Estimating mediated effects with survival data. In: New developments in psychometrics. Springer, 2003:405–412.10.1007/978-4-431-66996-8_46Search in Google Scholar

[43] Belanger CF, Hennekens CH, Rosner B, Speizer FE. The nurses’ health study. Am. J. Nursing 1978;78:1039–1040.10.2307/3462013Search in Google Scholar

[44] Wolf AM, Hunter DJ, Colditz GA, Manson JE, Stampfer MJ, Corsano KA, Rosner B, Kriska A, Willett WC. Reproducibility and validity of a self-administered physical activity questionnaire. Int. J. Epidemiol. 1994;23:991–999.10.1093/ije/23.5.991Search in Google Scholar PubMed

[45] McCormack VA, dos Santos Silva I. Breast density and parenchymal patterns as markers of breast cancer risk: a meta-analysis. Cancer Epidemiol. Biomarkers Prev. 2006;15:1159–1169.10.1158/1055-9965.EPI-06-0034Search in Google Scholar PubMed

[46] Rice MS, Bertrand KA, VanderWeele TJ, Rosner BA, Liao X, Adami H-O, Tamimi RM. Mammographic density and breast cancer risk: a mediation analysis. Breast Cancer Res. 2016;18:94.10.1186/s13058-016-0750-0Search in Google Scholar PubMed PubMed Central

[47] Spiegelman D, Hertzmark E. Easy sas calculations for risk or prevalence ratios and differences. Am. J. Epidemiol. 2005;162:199–200.10.1093/aje/kwi188Search in Google Scholar PubMed

[48] Stone CJ. The dimensionality reduction principle for generalized additive models. Ann. Stat. 1986:590–606.10.1214/aos/1176349940Search in Google Scholar

[49] Steen J, Loeys T, Moerkerke B, Vansteelandt S. medflex: An r package for flexible mediation analysis using natural effect models. J. Stat. Softw. 2017;76:1–46.10.18637/jss.v076.i11Search in Google Scholar


Supplemental Material

The online version of this article offers supplementary material (https://doi.org/10.1515/ijb-2017-0006).


Published Online: 2017-9-20

© 2017 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 29.3.2024 from https://www.degruyter.com/document/doi/10.1515/ijb-2017-0006/html
Scroll to top button