Improved sensitivity bounds for mediation under unmeasured mediator–outcome confounding

1 Introduction

Mediation analysis is central in biostatistics, epidemiology, and related disciplines. The aim of such analysis is to decompose a total effect of an exposure on an outcome into an indirect effect, relayed through a particular mediator, and a direct (not through the mediator) effect. A common problem in mediation analysis is the presence of unmeasured confounding (common causes) of the mediator and the outcome, which makes direct and indirect effects non-identifiable, even if the exposure is not confounded with the mediator and the outcome [1,2].

Ding and Vanderweele [3], henceforth DV, proposed to handle this non-identifiability problem by computing bounds on the direct and indirect effects, i.e., ranges guaranteed to include the true values of the effects. In their proposed sensitivity analysis, the analyst specifies two sensitivity parameters, which measure the “maximal” conditional association between the confounders and the outcome and between the confounders and the exposure, given the mediator. Together with the observed data distribution, these sensitivity parameters translate into bounds on the direct and indirect effects according to simple analytic formulas.

We extend DV’s results in three directions. First, we derive the feasible space (i.e., the space of logically possible values) for the sensitivity parameters, given the observed data distribution. This space determines what values an analyst may, and may not, consider in the computation of the bounds. Second, we derive novel bounds on the direct and indirect effects. Even though DV referred to their bounds as “sharp,” they are not sharp in the absolute sense of being the tightest possible, given their assumptions. Our novel bounds are at least as sharp, and sometimes sharper, than DV’s bounds. We show by a numerical example that the improvement over DV’s bounds can be substantial. Furthermore, whereas DV only derived a lower bound for the direct effect and an upper bound for the indirect effect, we derive lower and upper bounds on both effects. Third, we derive Cornfield inequalities that the unmeasured confounders must satisfy to explain away an “observed” effect. These inequalities are at least as sharp, and sometimes sharper, than similar Cornfield inequalities for the direct effect derived by DV. Furthermore, whereas DV’s Cornfield inequalities only apply to direct effects of positive magnitude, our Cornfield inequalities apply to both direct and indirect effects, having either positive or negative magnitudes.

In line with previous literature, we ignore uncertainty due to sampling variability, and assume that the exposure and outcome are binary. However, whereas DV allowed for a mediator of arbitrary type, we restrict attention to binary mediators. Proofs of all theorems, and R code for all analyses, are provided in the Appendix.

2 Notation, definitions, and assumptions

Let A , M , and Y denote the binary exposure, mediator, and outcome, respectively. In addition to these, data may be available on a set of covariates C . However, since all DV’s results, and ours, are derived within strata of covariates, we keep the conditioning on C implicit in the notation. Therefore, the observed data distribution is given by p ( Y , M , A ) . Let M a and Y a be the counterfactual mediator and outcome, respectively, where A set to a , and let Y a M a ′ be the counterfactual outcome where A set to a and M set to M a ′ . Finally, define p a a ′ = p ( Y a M a ′ = 1 ) .

Like DV, we allow for unmeasured confounding of the mediator and the outcome, but we assume that the exposure is not confounded with the mediator and the outcome (given C ). Specifically, we assume that data were generated by a nonparametric structural equation model represented by the causal diagram in Figure 1. In the Appendix, we provide an SWIG illustrating the setting [4]. Here, U is the whole set of (unmeasured) common causes of M and Y . This model implies the following three assumptions. Consistency of counterfactuals:

Conditional independence between the set of mediator counterfactuals and the set of outcome counterfactuals, given U :

Independence between the set of counterfactuals and the exposure and between U and the exposure:

DV made essentially the same assumptions but phrased slightly differently.

Figure 1

Assumed causal diagram.

The target parameters are the natural direct effect (NDE) and natural indirect effect (NIE), which are defined as some contrasts

and

where g is a link function. The identity link, log link, and logit link give the effects as risk differences, log risk ratios, and log odds ratios, respectively. DV only considered the risk difference and risk ratio, but their results do in fact apply more generally.

DV showed that p 00 and p 11 are identified and equal p ( Y = 1 ∣ A = 0 ) and p ( Y = 1 ∣ A = 1 ) , respectively. Thus, identification of NDE and NIE crucially relies on identification of p 10 . Define

and the “observed” NDE and NIE

and

DV showed that, in the absence of unmeasured confounding U , p 10 is identified and equal to p 10 obs , so that NDE obs = NDE and NIE obs = NIE . However, in the presence of unmeasured confounding, p 10 , NDE and NIE are not identified without further assumptions.

3 DV’s lower bound on p 10

DV defined the sensitivity parameters RR U Y = max m RR U Y ∣ m and RR A U = max m RR A U ∣ m , where

and

They defined the bounding factor

which is a monotonically increasing function of RR U Y and RR A U . They showed that, given BF , p 10 is bounded from below by

Once the analyst has specified RR U Y and RR A U , these can be used to compute BF and the lower bound on p 10 . This bound can then be converted into a lower bound on NDE and an upper bound on NIE. For instance, a lower bound on NDE and an upper bound on NIE as risk ratios are given by l D V ∕ p ( Y = 1 ∣ A = 0 ) and p ( Y = 1 ∣ A = 1 ) ∕ l D V , respectively.

4 Feasible space of DV’s sensitivity parameters

To use DV’s bounds in practice, it is important to know the feasible space of their sensitivity parameters, given the observed data distribution. This space is established by the following theorem.

Part (ii) of the theorem implies that the feasible space of one of the sensitivity parameters is neither restricted by the value of the other sensitivity parameter, nor by the observed data distribution. This means that, regardless of what data distribution one may observe, one should consider all values ≥ 1 , as stated by part (i) of the theorem, as logically possible for both sensitivity parameters. This does not mean that all values ≥ 1 are plausible for any given application; hopefully, one may be able to use subject matter knowledge to significantly reduce this space when computing the bounds.

5 A sharper lower bound on p 10

It follows from results by Sjölander [5] that p 10 is bounded from below by

where

This bound is “assumption-free” in the sense that it does not require specification of BF (or any other sensitivity parameter). It follows from Theorem 1 that, for any given observed data distribution p ( Y , M , A ) , BF can be arbitrarily large. Furthermore, DV’s lower bound on p 10 bound approaches 0 as BF goes to infinity. Thus, it is possible to construct a scenario (e.g., a joint distribution p ( Y , M , A , U ) ) where DV’s lower bound is arbitrarily close to 0, whereas the assumption-free bound is substantially larger than 0, i.e., a scenario where the assumption-free bound is larger than DV’s bound.

A simple way to improve DV’s lower bound is to replace it with the assumption-free bound, whenever the latter is larger. However, the bound can be improved even further, as stated in the following theorem.

As an example, suppose that U is binary with p ( U = 1 ) = 0.54 . Suppose further that p ( M = 1 ∣ A , U ) and p ( Y = 1 ∣ A , M , U ) are as in Tables 1 and 2, respectively. For this example, the true values of p 10 and BF are 0.81 and 1.71, respectively. Suppose now that only the marginal (over U ) distribution p ( Y , M , A ) is given, but that we are able to correctly specify the true value of BF . From these, we obtain l D V = 0.44 , l A F = 0.59 , and l † = s 1 = 0.65 . Thus, for this example, the assumption-free lower bound offers substantial improvement over DV’s lower bound, and the lower bound in Theorem 2 offers a substantial improvement over the assumption-free bound.

6 An upper bound on p 10

To construct an upper bound on p 10 , we need to define a new sensitivity parameter relating the exposure to the confounders. In analogy with RR A U , we define the sensitivity parameter RR ˜ A U = max m RR ˜ A U ∣ m , where

and we define the corresponding bounding factor

The following two theorems establish the feasible space for { RR ˜ A U , RR U Y } , given the observed data distribution, and an upper bound on p 10 , given BF ˜ .

In Theorem 4, we have used the notation u D V to indicate that this term is “analogous” to l D V in Theorem 2, in the sense that u D V is derived by using the same arguments as in the proof for the term l D V in Theorem 2. In the same sense, the terms v ˜ , s ˜ 0 , and s ˜ 1 are analogous to the terms v , s 0 , and s 1 , respectively, in Theorem 2. The term u A F is an assumption-free upper bound p 10 , which follows from results by Sjölander [5]. The reason why the corresponding term l A F does not appear in Theorem 2 is that this term is redundant, given { l D V , s 0 , s 1 } , since l D V ≥ 0 , s 0 ≥ q 0 , and s 1 ≥ q 1 ; hence, l A F ≤ max { l D V , s 0 , s 1 } .

Continuing on the example from Section 5, we have that u † = s ˜ 1 = 0.90 . Thus, for this example, we can tell that true value of p 10 is somewhere in the range ( 0.65 , 0.90 ) , provided that we are able to correctly specify the true value of BF .

7 Cornfield inequalities

DV used their lower bound l D V to derive Cornfield inequalities for NDE. These inequalities are lower bounds on the common value RR A U = RR U Y = δ required to “explain away” an observed NDE with a positive magnitude. Here, we say that an observed NDE has positive magnitude if p 10 obs > p ( Y = 1 ∣ A = 0 ) , and that the observed NDE is explained away by δ if the lower bound for the true NDE includes the null value, e.g., 0 for the risk difference and 1 for the risk ratio. For the risk ratio, DV’s Cornfield inequality has a simple analytic form, given by δ ≥ NDE obs + NDE obs ( NDE obs − 1 ) . This Cornfield inequality is valid in the sense that δ must indeed exceed this lower bound in order to explain away an observed NDE. However, since this Cornfield inequality is derived from the non-sharp bound l D V , it is also not sharp and may thus give an unnecessarily pessimistic assessment of the study at hand.

To illustrate this, we again return to the numeric example from Section 5. For this example, p 10 obs = 0.76 , p ( Y = 1 ∣ A = 0 ) = 0.54 , NDE obs = p 10 obs ∕ p ( Y = 1 ∣ A = 0 ) = 1.39 , and NDE obs + NDE obs ( NDE obs − 1 ) = 2.13 . Thus, in order to explain away the observed NDE, a common value of RR A U and RR U Y must be at least 2.13. This is a rather modest value, which may give the impression that the observed NDE is quite sensitive to unmeasured confounding. However, the assumption-free lower bound l A F is equal to 0.59, which is larger than p ( Y = 1 ∣ A = 0 ) . Hence, even if only the marginal (over U ) distribution p ( Y , M , A ) is known, and we are not able to specify BF , we can still in this example tell that there is a true non-null NDE. Put differently, no degree of unmeasured confounding could ever explain away the observed NDE in this example. This conclusion is much stronger than that obtained from DV’s Cornfield inequality.

DV’s Cornfield inequalities assume that the observed NDE has positive magnitude, and they are not applicable to NIE. Using the bounds ( l † , u † ) we can obtain Cornfield inequalities for both NDE and NIE, having either positive or negative magnitudes.

If the observed NDE has positive magnitude and l † = l D V , then the Cornfield inequality in Theorem 5 reduces to DV’s Cornfield inequality, but is otherwise sharper. The equations defining BF ( a ) and BF ˜ ( a ) in Theorem 5 have no solution if p ( Y = 1 ∣ A = a ) is outside the range ( l † , u † ) . This means that no degree of unmeasured confounding can explain away the observed NDE (for a = 0 ) or NIE (for a = 1 ), as in the numeric example above.

8 Simulation

We carried out a simulation to assess the numerical performance of the proposed bounds, and the potential for the relatively complex bounds ( l † , u † ) to improve on the simpler bounds ( l D V , u D V ) . We considered a binary confounder U and generated 1,000,000 distributions p ( Y , A , M , U ) from the model

For each distribution p ( Y , A , M , U ) , we computed the probabilities p 00 , p 10 , and p 11 , the bounding factors BF and BF ˜ , and the bounds l D V , l † , u D V , and u † . When computing the bounds, we used the true bounding factors BF and BF ˜ multiplied with a constant k ≥ 1 . The value k = 1 represents the ideal scenario where the analyst uses the true bounding factors. This scenario is somewhat unrealistic, since the analyst would rarely know the true bounding factors in practice. Thus, we computed the bounds on both k = 1 and k = 1.3 . The latter is intended to represent the possibly more realistic scenario where the analyst provides a conservative guess of the bounding factors, in order not to invalidate the bounds.

We first assessed the distance between the bounds and the true value of p 10 , by constructing the standardized differences

and

Since ( l † , u † ) are at least as sharp as ( l D V , u D V ) , we have that 0 ≤ ( Δ l , Δ u ) ≤ 1 . The plots in Figure 2 show the empirical (over the generated distributions) complementary cumulative distribution functions for Δ l (left plot) and Δ u (right plot). We observe that even though the bounds ( l † , u † ) are usually not sharper than ( l D V , u D V ) , they sometimes are and the improvement can be substantial. For instance, with 4% and 6% probability, l † is more than 20% closer to p 10 than l D V , for k = 1 and k = 1.3 , respectively. The corresponding probabilities for u † are larger; 8% and 13% for k = 1 and k = 1.3 , respectively.

Figure 2

Complementary cumulative distribution functions for the standardized differences Δ l (left plot) and Δ u (right plot), for k = 1 (solid lines) and k = 1.3 (dashed lines).

Since l D V and u D V approach 0 and 1 as BF and BF ˜ approach infinity, one would intuitively expect the performance of ( l † , u † ) , relative to ( l D V , u D V ) , to increase with ( BF , BF ˜ ) . Figure 3 shows the same functions as in Figure 2, for k = 1 , but only based on those distributions p ( Y , A , M , U ) for which BF (for the lower bounds) and BF ˜ (for the upper bounds) exceed a threshold equal to 1, 2, 5, and 10. As expected, we observe that the complementary cumulative distribution functions increase as the bounding factors increase. This indicates that the stronger the unmeasured confounding, the higher the probability that the bounds ( l † , u † ) improves substantially on the bounds ( l D V , u D V ) .

Figure 3

Complementary cumulative distribution functions for the standardized differences Δ l (left plot) and Δ u (right plot), stratified by thresholds for BF and BF ˜ .

We next assessed the ability of the bounds to reject the null hypotheses that NDE = 0 and NIE = 0. The null hypothesis NDE = 0 is rejected by a lower (upper) bound on p 10 if the bound is larger (smaller) than p 00 . Similarly, the null hypothesis NIE = 0 is rejected by a lower (upper) bound on p 10 if the bound is larger (smaller) than p 11 . Tables 3 and 4 show the empirical probabilities of these rejections, for ( l D V , l † ) and ( u D V , u † ) , respectively. We observe that the rejection probabilities are very similar for l D V and l † and also very similar for u D V and u † , for both k = 1 and k = 1.3 , and for both NDE and NIE. This indicates that, even though the bounds ( l † , u † ) are potentially much sharper than the bounds ( l D V , u D V ) , their additional impact is modest on the power to reject the null hypotheses.

9 Application

We illustrate the proposed bounds by analyzing data from the Lipid Research Clinics Coronary Primary Prevention (LRC) trial [6]. These data were previously used by Cai et al. [7] and Sjölander [5] to illustrate the assumption-free bounds on controlled direct effects and NDEs, respectively. In the LRC trial, 1888 subjects were randomized to cholestyramine treatment and 1918 subjects to placebo. During a follow-up period of 1 year, each coronary heart disease (CHD) event was recorded. At the end of the follow-up, cholesterol levels were recorded for each subject. Both Cai et al. [7] and Sjölander [5] dichotomized cholesterol levels as “low” ( < 280 mg/dl) or “high” ( ≥ 280 mg/dl). The data are displayed in Table 5.

We coded all variables so that the “1” represents the “desirable” level, i.e., treatment ( A = 1 ), low cholesterol ( M = 1 ), and CHD absent ( Y = 1 ). We then computed each element of the min/max expressions for the bounds on p 10 given in Theorems 2 and 4. We finally constructed bounds on NDE and NIE as risk ratios, for each of these elements separately. For instance, the lower bound on NDE and the upper bound on NIE obtained from the element l D V are given by l D V ∕ p ( Y = 1 ∣ A = 0 ) and p ( Y = 1 ∣ A = 1 ) ∕ l D V , respectively.

Figure 4 shows the lower (left column) and upper (right column) bounds on NDE (top row) and NIE (bottom row), as functions of BF (for NDE) and BF ˜ (for NIE). In the absence of confounding ( BF = BF ˜ = 1 ), the lower and upper bounds collapse at the “observed effects” NDE = 1.02 and NIE = 1.01. These effects are explained away by a tiny amount of confounding; evaluating the Cornfield inequality in Theorem 5 gives that δ = 1.01 for both NDE and NIE. Further, the amount of assumed confounding has quite dramatic effects on some of the bounds. At BF = 1 , the bound l D V ∕ p ( Y = 1 ∣ A = 0 ) is largest, and thus most informative, among all lower bounds on NDE. However, already at BF = 1.6 , l D V ∕ p ( Y = 1 ∣ A = 0 ) becomes smaller than s 1 ∕ p ( Y = 1 ∣ A = 0 ) , and at BF = 2.6 , l D V ∕ p ( Y = 1 ∣ A = 0 ) becomes smaller than the assumption-free bound l A F ∕ p ( Y = 1 ∣ A = 0 ) . We observe a similar pattern for the upper bounds on NIE, where the bound p ( Y = 1 ∣ A = 1 ) ∕ l D V is smallest at BF = 1 , but becomes larger than p ( Y = 1 ∣ A = 1 ) ∕ s 1 at BF = 1.6 , and larger than the assumption-free bound p ( Y = 1 ∣ A = 1 ) ∕ l A F at BF = 2.6 . The upper bounds on NDE are less diverse; apart from u D V ∕ p ( Y = 1 ∣ A = 0 ) , all bounds are identical ( = 1 ∕ p ( Y = 1 ∣ A = 0 ) = 1.1 ) for all values of BF ˜ , and thus not be distinguished from each other in the top-right plot of Figure 4. The bound u D V ∕ p ( Y = 1 ∣ A = 0 ) is slightly smaller ( = 1.02 ) at BF ˜ = 1 , but becomes equal to all other bounds at BF ˜ = 1.1 . We observe a similar pattern for the lower bounds on NIE; apart from p ( Y = 1 ∣ A = 1 ) ∕ u D V , all bounds are identical ( = p ( Y = 1 ∣ A = 1 ) ∕ 1 = 0.93 ) for all values of BF ˜ . The bound p ( Y = 1 ∣ A = 1 ) ∕ u D V is slightly larger ( = 1.01 ) at BF ˜ = 1 , but becomes equal to all other bounds at BF ˜ = 1.1 .

Figure 4

Bounds on NDE and NIE, for the LRC data.

10 Discussion

The proposed bounds are applicable to any direct and indirect effects that can be written as contrasts of probabilities for binary outcomes. A useful extension would be to derive analogous bounds on contrasts of survival probabilities, with time-to-event outcomes.

We have shown that the proposed bounds are at least as sharp as those derived by DV. However, we have not been able to prove (or disprove) that our bounds are sharp in an absolute sense, i.e., that all values inside the bounds are logically possible, given the assumptions. We recognize this as an interesting topic for future research.