3.1 Counterfactual independence assumptions

The natural direct and indirect AFs are causal (counterfactual) parameters, and consistent estimation of these from observational data requires that we link the counterfactual variables to the observed variables, and that we make appropriate control for confounders. We thus make the standard ‘consistency assumption’ (Pearl 2009) that an intervention to set variable(s) V to value v does not affect the outcome, if V is already observed to be equal v:

We further assume that data have been collected on a set of variables C, which is sufficient for confounding control. Formally, we assume that that YA=0 is conditionally independent of A, given C:

and that YA=0,M is conditionally independent of A, given (M,C):

Assumption (5) requires that C is sufficient for control of A-Y confounding, and assumption (6) requires that C is sufficient for control of both A-Y confounding and M-Y confounding. In addition, assumption (6) requires that there are no confounders for M and Y that are affected by A. Assumptions (5) and (6) would thus, for instance, hold under the causal diagram in Figure 1. We refer to Pearl (2009) for general techniques to evaluate whether assumptions (5) and (6) hold under other causal diagrams. We note that the graphical identification criteria proposed by Shpitser and VanderWeele (2011) for natural direct and indirect effects also suffice for identification of our parameters of interest.

Figure 1:

A causal diagram under which assumptions (6) and (5) hold.

Under assumptions (4)-(6) we can rewrite the natural direct and indirect AFs, so that they are free from counterfactual probabilities. Towards this end we define ψ=E{p(Y=1|A=0,M,C)}, where the outer expectation is taken over the marginal distribution of (M,C). We can write the counterfactual probability p(YA=0,M=1) as

where the first equality follows from the law of iterated expectations, the second equality follows from assumption (6) and the third equality follows from assumption (4). Similarly, define ϕ=E{p(Y=1|A=0,C)}, where the outer expectation is taken over the marginal distribution of C. We can write the counterfactual probability p(YA=0=1) as

where the first equality follows from the law of iterated expectations, the second equality follows from assumption (5) and the third equality follows from assumption (4).

Define π=p(Y=1). It follows that we can write the natural direct and indirect AFs as

and

The expressions in (9) and (10) are both free from counterfactual variables, which enables us to estimate these from observational data, as described in the following two sections.

We end this section by emphasizing the importance of appropriate confounding control, when applying the estimation techniques developed in the following sections. Arguably, there is always unmeasured confounding in real observational studies. However, with a careful study design, the degree of unmeasured confounding may be reasonably low, in which case the crucial assumptions (5) and (6) may hold as an approximation.

3.2 Cross-sectional/cohort studies

In cross-sectional studies and cohort studies, data consist of an iid sample of size n from p(Y,A,M,C). Under this sampling scheme we may estimate π, ψ and ϕ separately, and plug the estimates into eqs. (9) and (10).

Trivially, we may use the sample prevalence πˆ as an estimate of the population prevalence π. A simple way to estimate ψ and ϕ would be to use regression models for p(Y|A,M,C) and p(Y|A,C), respectively. Since Y is binary it is natural to use logistic regression models, which can be fitted with standard ML methods. From the fitted models we may obtain predictions of p(Y=1|A=0,M,C) and p(Y=1|A=0,C) for each subject (i.e. for each observed value of (M,C) and C, respectively). By averaging these predictions we obtain ML estimates of ψ and ϕ, respectively.

A disadvantage of this ML method is that p(Y|A,M,C) and p(Y|A,C) are not variation independent, but related through the integral

This feature is somewhat awkward, since it implies that it may difficult to postulate models for p(Y|A,M,C) and p(Y|A,C) that are simultaneously correct. For instance, a logistic form for p(Y|A,M,C) generally implies a more complex functional form for p(Y|A,C), and vice versa.

One solution to this problem would be to replace the model for p(Y|A,C) with a model for p(M|A,C), then use the relation in eq. (11) to derive the implied model for p(Y|A,C). Since p(Y|A,M,C) and p(M|A,C) are variation independent, the implied model for p(Y|A,C) is guaranteed to be compatible with the model for p(Y|A,M,C). However, this method has important disadvantages for non-binary mediators. First, the integral in eq. (11) requires a model for the full conditional distribution of M (e.g. a mean model does not suffice), which may be difficult to well specify. Second, the integral in eq. (11) does not generally have an analytic solution, which makes the method computationally demanding.

To reduce sensitivity to model misspecification, while avoiding having to model the mediator, we instead propose to use DR estimators for both ψ and ϕ. These estimators use one model for the outcome and one model for the exposure, and are consistent if either of these are correctly specified, not necessarily both.

Specifically, to estimate ψ we use one model for p(Y|A,M,C) and one model for p(A|M,C). After having fitted these models we estimate ψ as

where pˆ(Y=1|A=0,Mi,Ci) and pˆ(Ai|Mi,Ci) are model-based predictions for subject i. It follows from more general results in Robins (2000) that ψˆ is DR, in the sense that it is consistent if the model for p(Y|A,M,C) or the model for p(A|M,C) is correct, not necessarily both. Similarly, to estimate ϕ we use one model for p(Y|A,C) and one model for p(A|C). After having fitted these models we estimate ϕ as

where pˆ(Y=1|A=0,Ci) and pˆ(Ai|Ci) are model based predictions for subject i. Again, it follows from Robins (2000) that ϕˆ is DR, in the sense that it is consistent if the model for p(Y|A,C) or the model for p(A|C) is correct, not necessarily both.

We note that p(Y|A,M,C) and p(A|C) are variation independent. This means that the DR estimators give the researcher a ‘genuine chance’ of having a correct model when estimating ψ, while simultaneously having a correct model when estimating ϕ.

In Appendix A we use the theory for M-estimators (Stefanski and Boos 2002) to show that (πˆ,ψˆ,ϕˆ) has an asymptotic normal distribution, and we derive the asymptotic variance-covariance matrix for (πˆ,ψˆ,ϕˆ). We subsequently use the delta method to derive the asymptotic variances for NDAFˆ and NIAFˆ, obtained by plugging (πˆ,ψˆ,ϕˆ) into eqs. (9) and (10), respectively.

Table 1 summarizes what models are required to be correct, to guarantee consistency of the ML and DR estimators of NDAF, NIAF and AF. In this table we have explicitly spelled out NDAF, NIAF and AF as functions of the parameters being estimated by the models.

3.3 Case-control studies

In case-control studies, data consist of n1 iid cases sampled from p(A,M,C|Y=1), and n0 iid controls sampled from p(A,M,C|Y=0). Under this sampling scheme, it is not possible to obtain consistent estimates of π, ψ and ϕ, since this sampling scheme mutilates the outcome distribution. Instead, we use results in Bruzzi et al. (1985) to approximate NDAF and NIAF as functions of conditional odds ratios, which can be consistently estimated from case-control data (Prentice and Pyke 1979).

Let OR(M,C)=p(Y=1,A=1|M,C)p(Y=0,A=0|M,C)p(Y=0,A=1|M,C)p(Y=1,A=0|M,C) be the conditional odds ratio for A and Y, given (M,C), and define ψ˜=E{OR−A(M,C)|Y=1}, where the outer expectation is taken over the conditional distribution of (A,M,C), given Y=1. Similarly, let OR(C)=p(Y=1,A=1|C)p(Y=0,A=0|C)p(Y=0,A=1|C)p(Y=1,A=0|C) be the conditional odds ratio for A and Y, given C, and define ϕ˜=E{OR−A(C)|Y=1}, where the outer expectation is taken over the conditional distribution of (A,C), given Y=1. It follows from results in Bruzzi et al. (1985) that ψ/π≈ψ˜ and ϕ/π≈ϕ˜, provided that the outcome is rare. We thus have that

and

A simple way to estimate ψ˜ and ϕ˜ is to use logistic regression models for p(Y|A,M,C) and p(Y|A,C), respectively, which can be fitted with standard ML techniques. From the fitted models we may obtain predictions of OR−A(M,C) and OR−A(C) for each subject (i.e. for each observed value of (A,M,C) and (A,C), respectively). By averaging these predictions among the cases we obtain ML estimates of ψ˜ and ϕ˜, respectively.

However, to protect against model misspecification we may also use DR estimators. Specifically, to estimate ψ˜ we may use one logistic regression model for p(Y|A,M,C) and one logistic regression model for p(A|Y,M,C), which are constrained by having the same conditional odds ratio OR(M,C). We let θ be the parameter (vector) indexing OR(M,C) in these models. A standard choice of models would be logit{p(Y=1|A,M,C)}=θA+α0+α1C+α2M and logit{p(A=1|Y,M,C)}=θY+β0+β1C+β2, in which logit{OR(M,C;θ)}=θ is assumed to be constant across levels of M and C. Other examples are logit{p(Y=1|A,M,C)}=θ1A+θ2AM+α0+α1C+α2M and logit{p(A=1|Y,M,C)}=θ1Y+θ2YM+β0+β1C+β2M, in which logit{OR(M,C;θ)}=θ1+θ2M is allowed to depend linearly on M. The two logistic regression models can be fitted with ML techniques, and the fitted models may be combined, as described by Tchetgen Tchetgen, Robins, and Rotnitzky (2010), to produce an estimate of θ that is consistent if either model is correct, not necessarily both. We use the DR estimate of θ to obtain a prediction of OR−A(M,C;θ) for each subject. Finally, we obtain a DR estimate of ψ˜ by averaging these predictions among the cases. We obtain a DR estimate of ϕ˜ in a similar fashion, by using models for p(Y|A,C) and p(A|Y,C) instead of models for p(Y|A,M,C) and p(A|Y,M,C).

In Appendix A we use the theory for M-estimators (Stefanski and Boos 2002) to show that (ψ˜ˆ,ϕ˜ˆ) has an asymptotic normal distribution, and we derive the asymptotic variance-covariance matrix for (ψ˜ˆ,ϕ˜ˆ). We subsequently use the delta method to derive the asymptotic variances for NDAFˆ and NIAFˆ, obtained by plugging (ψ˜ˆ,ϕ˜ˆ) into eqs. (12) and (13), respectively.

Table 2 summarizes what models are required to be correct, to guarantee consistency of the ML and DR estimators of NDAF, NIAF and AF. In this table we have explicitly spelled out NDAF, NIAF and AF as functions of the parameters being estimated by the models.

We end this section with a technical remark. The validity of the OR-based estimators developed above depends on the approximate relations ψ/π≈ψ˜ and ϕ/π≈ϕ˜; the better approximations, the less (asymptotic) bias in the estimators. It can be shown (Bruzzi et al. 1985) that these approximations are exact, if the odds ratios OR(M,C) and OR(C) in ψ˜ and ϕ˜, respectively, are replaced by the corresponding risk ratios. Indirectly, the validity of the OR-based estimators is thus related to the concept of non-collapsibility, which is a well-known feature of the odds ratio (Greenland, Robins, and Pearl 1999). When the outcome is rare, the odds ratio approximates the risk ratio, which is collapsible. As the outcome becomes more common, the odds becomes increasingly more non-collapsible, and increasingly different from the risk ratio.

4.1 Cross-sectional/cohort studies

To investigate the performance of the proposed estimators for cross-sectional/cohort studies (Section 3.2) we generated 10,000 samples, of n=1000 iid observations each, from the model

For this model we have that p(Y=1)=0.02, NDAF=0.2, NIAF=0.05 and AF=0.25. Each sample was analyzed with the working models

When analyzing the simulated data, four scenarios were considered. In scenario I, no constraints were imposed on the working model parameters. In scenario II, the interaction terms α12 and γ12 were omitted from the working models. In scenario III, the interaction terms β12 and δ12 were omitted. In scenario IV, all interaction terms were omitted. For this analysis scheme, there are two types of model misspecification. First, the functional form for the probabilities p(A|M,C) and p(Y|A,C) is somewhat misspecified across all scenarios, since the logistic forms for p(Y|A,M,C) and p(A|C) in eq. (22), together with the normal distribution for M|(A,C), do not imply logistic forms for p(A|M,C) and p(Y|A,C). Second, the linear predictor for a particular working model is misspecified in a given scenario, if the interaction term between C1 and C2 is omitted from the working model in that scenario.

For each of the four scenarios scenario we used the working models to calculate the ML estimators and DR estimators of NDAF, NIAF and AF=NDAF+NIAF, as described in Section 3.2, together with their standard errors, as described in Appendix A. We finally calculated the mean (over the 10,000 samples) estimates, the mean standard errors and the standard deviation of the estimates. We note that, from the perspective of the ML estimators, scenarios I and III are identical and scenarios II and IV are identical, since the ML estimators only use the models for p(Y|A,M,C) and p(Y|A,C).

Table 3 displays the results. As expected, the ML estimator of NDAF is unbiased in scenarios I and III, where the model for p(Y|A,M,C) is correct. Interestingly, the ML estimators of NIAF and AF are unbiased in these scenarios as well, even though the functional form for p(Y|A,C) is misspecified for these (and all other) scenarios. The ML estimators of NDAF, NIAF and AF are biased in scenarios II and IV though, where the linear predictors for p(Y|A,M,C) and p(Y|A,C) are misspecified. We observe a similar pattern for the DR estimators. The DR estimator of NDAF is unbiased in scenarios I, II and III, where the linear predictor for either (or both) p(Y|A,M,C) and p(A|M,C) is (are) correct, even though the functional form for p(A|M,C) is misspecified. The DR estimators of NIAF and AF are also unbiased in scenarios I, II and III, where the linear predictor for either (or both) p(Y|A,C) and p(A|C) is (are) correct, even though the functional form for p(Y|A,C) is misspecified. In scenario IV, the DR estimators of both NDAF, NIAF and AF are biased. These results indicate that the DR estimators are indeed more robust to model misspecification than the ML estimators, and that both estimators have some robustness against misspecification of the functional form of the modeled probabilities. The mean standard errors agree well with the standard deviation of the estimates, for both the ML and DR estimators and across all scenarios. The mean standard error for the DR estimator is close to the mean standard error for the ML estimator across all scenarios, which indicates that the additional roubstness of the DR estimator comes at a relatively minor loss in precision.

We caution the reader that the observed robustness against functional form misspecification is not guaranteed by theory, and may depend heavily on the data generating model. In a setting similar to our, VanderWeele and Vansteelandt (2010) showed that functional form misspecification is minor when the outcome is rare and the mediator is continuous with constant variance, which is precisely the model we have simulated from above. Thus, one may suspect that bias due to functional form misspecification would be more pronounced if the outcome was more common and/or the mediator had a skewed distribution. We investigated this in an additional simulation (see Appendix B), but observed no more bias than in the simulation above. This gives some further support to the conclusions that, while misspecification of the linear predictor can give rise to serious bias, the impact of function form misspecification may be less severe.

Much of the literature on mediation has focused on the scenario where there are exposure-mediator interactions (see, for instance, VanderWeele (2015) and the references therein). We have added a simulation in Appendix B, in which exposure-mediator interaction is present. The conclusions from this simulation are again similar to those from the simulation above.

4.2 Case-control studies

To investigate the performance of the proposed estimators for case-control studies (Section 3.3) we generated a population of 10,000,000 subjects from the model in eq. (22). From this population we drew 10,000 samples of n1=500 cases and n0=500 controls. Each sample was analyzed with the working models

In these working models, OR(M,C)=exp(θ) and OR(C)=exp(η). We considered the same four scenarios as in the previous section. Thus, only the functional form for p(Y|A,M,C) is correct, since logistic forms for p(Y|A,M,C) and p(A|C), together with the normal distribution for M|(A,C), do not imply logistic forms for p(A|Y,M,C), p(Y|A,C) or p(A|Y,C). For each scenario we used the working models to calculate the ML estimators and DR estimators of NDAF, NIAF and AF=NDAF+NIAF, as described in Section 3.3, together with their standard errors, as described in Appendix A. We finally calculated the mean (over the 10,000 samples) estimates, the mean standard errors and the standard deviation of the estimates. Again we note that, from the perspective of the ML estimators, scenarios I and III are identical and scenarios II and IV are identical, since the ML estimators only use the models for p(Y|A,M,C) and p(Y|A,C).

Table 4 displays the results. We observe the same pattern as in the previous section. The ML estimators are unbiased in scenarios I and III, where the linear predictors for p(Y|A,M,C) and p(Y|A,C) are correct, but biased in scenarios II and IV where these linear predictors are misspecified. The DR estimator of NDAF is unbiased in scenarios I, II and III, where the linear predictor for either (or both) p(Y|A,M,C) and p(A|Y,M,C) is (are) correct, and the DR estimator of NIAF is also unbiased in scenarios I, II and III, where the linear predictor for either (or both) p(Y|A,C) and p(A|Y,C) is (are) correct. Misspecification of the functional form for the modeled probabilities does not seem to give bias for any of the ML or DR estimators.