Sharp Nonparametric Bounds for Decomposition Effects with Two Binary Mediators

Abstract In randomized trials, once the total effect of the intervention has been estimated, it is often of interest to explore mechanistic effects through mediators along the causal pathway between the randomized treatment and the outcome. In the setting with two sequential mediators, there are a variety of decompositions of the total risk difference into mediation effects. We derive sharp and valid bounds for a number of mediation effects in the setting of two sequential mediators both with unmeasured confounding with the outcome. We provide five such bounds in the main text corresponding to two different decompositions of the total effect, as well as the controlled direct effect, with an additional 30 novel bounds provided in the supplementary materials corresponding to the terms of 24 four-way decompositions. We also show that, although it may seem that one can produce sharp bounds by adding or subtracting the limits of the sharp bounds for terms in a decomposition, this almost always produces valid, but not sharp bounds that can even be completely noninformative. We investigate the properties of the bounds by simulating random probability distributions under our causal model and illustrate how they are interpreted in a real data example. Supplementary materials for this article are available online.


Introduction
In randomized trials with full compliance, the observed association between the intervention and the outcome has a causal interpretation as the total intervention effect through all possible pathways.Once a total intervention effect has been established, there is often additional interest in specific pathways and mechanisms through which the intervention may affect the outcome.For settings with a single binary mediator, Robins and Greenland [1992] used counterfactual arguments to provide a formal framework for reasoning about direct effects and indirect (mediated) effects.Pearl [2001] proposed counterfactual definitions of direct and indirect effects for a single mediator.Specifically, Pearl [2001] distinguished between the controlled direct effect (CDE), which sets the mediator to a fixed value for each subject, and the natural direct effect (NDE), which sets the mediator to a counterfactual value that may differ across subjects.Although the natural direct effect is more difficult to conceptualize, it has the appealing property that it adds up, together with the corresponding natural indirect effect (NIE), to the total effect.
A problem with such effect decompositions is that the separate effect components are often not identified from data.Even when randomization rules out intervention-outcome confounding, there may still be unmeasured confounders for the mediator(s) and the outcome.If so, then any attempt to estimate the direct intervention effect by controlling for the mediator(s) will open back-door paths through the unmeasured confounders, thereby inducing an association between the intervention and the outcome [Robins and Greenland, 1992].
When causal effects are not identified, bounds for the target parameter of interest, i.e. a range that is guaranteed to include the true parameter value given the observed data, can be used to reduce the possible range of values that need to be considered.Balke and Pearl [1994] developed a method for deriving bounds for simple causal estimands based on linear programming techniques.For settings with a single binary mediator, Cai et al. [2007] used this linear programming technique to derive bounds on the CDE, assuming that the intervention is (or can be considered) randomized, but allowing for unmeasured confounding of the mediator and the outcome.Sjölander [2009] provided analogous bounds for the NDE.Kaufman et al. [2009] provided bounds for both the CDE and the NDE while relaxing the assumption about the exposure being randomized.Other authors have derived bounds for the CDE and NDE under certain monotonicity assumptions [VanderWeele, 2011, Chiba, 2010], and for settings where the mediator has more than two levels [Miles et al., 2017].
Recently, there has been a growing interest in effect decompositions with multiple mediators [e.g Avin et al., 2005, Albert and Nelson, 2011, VanderWeele and Vansteelandt, 2014, Daniel et al., 2015, Steen et al., 2017].However, this line of literature has focused on appropriate definitions and sufficient criteria for identification; to our knowledge, bounds for the nonidentified effect components in settings with more than one mediator have not been published.This is an important gap in the literature, since the introduction of multiple mediators poses identification problems that are not present in settings with a single mediator.Specifically, Daniel et al. [2015] showed that when a mediator has a direct effect on a subsequent mediator, the indirect effect through the former mediator is not generally nonparametrically identified, even in the complete absence of unmeasured confounders.
Hence, the importance of bounds is even stronger in setting with multiple mediators than in settings with a single mediator.
We derive bounds for the setting of a two-armed randomized trial with two causally ordered binary mediators that are confounded with the binary outcome of interest.Steen et al. [2017] proposed a two-way and a three-way decomposition of the total effect, for which we provide bounds for each component.The two-way decomposition is obtained by separating the total effect into the direct effect of the intervention on the outcome and total indirect effects through both of the mediators or the 'joint natural indirect effect', while the three-way decomposition further separates this joint natural indirect effect through the two mediators into two indirect effects.We also provide bounds for the controlled direct effects in this setting.In addition to providing sharp bounds for each term in the two-way and three-way decompositions, we demonstrate that in general the bounds for the separate terms of the decompositions cannot be combined to yield sharp bounds on their sum or difference.The exception is in the two-way decomposition when subtracting the limits of the joint indirect effect or the direct effect from the identified total effect, which we prove will always provide sharp bounds for the remaining effect.
Further decompositions have been derived, and as is shown in Daniel et al. [2015], one of the indirect effects in the three-way decomposition can be further decomposed into the effect through the second mediator due to the effect of intervention directly on that mediator and the effect of the intervention through the first mediator on the second, resulting in a four-way decomposition.In addition to the bounds we provide in the main text, we provide bounds for the exhaustive list of the thirty-two terms appearing in the twenty-four possible decompositions of the total intervention effect into the four-way decompositions of Daniel et al. [2015] in the supplementary materials.For each set of bounds we also provide easily downloadable R functions to facilitate their use.
In practice, there is a tendency towards point estimation, even in the absence of a discussed or well defined estimand or even if that estimand is only identifiable under strong untestable assumptions.Although sensitivity analyses is sometimes used to mitigate con-cerns about assumptions, these procedures rarely provide an assumption free range of the true causal effect.Due to the rising interest in the "estimands framework" for randomized clinical trials [Lipkovich et al., 2020], assumption free, i.e. nonparametric, bounds may become part of the standard for clinical trials analysis.In such settings, there are often multiple binary mediators, such as initiation of rescue mediation, withdrawal from treatment, or relapse or remission prior to death, making our bounds of practical relevance in these settings.Nonparametric bounds may sometimes be considered too wide to be useful in practice.In contrast, we believe that wide bounds are useful to present, since they highlight how little information the observed data contains about the target parameter, and how much a point estimate would have to rely on strong, potentially untestable assumptions.
The paper is organized as follows.In Section 2 we provide our notation and outline our estimands and settings of interest.In Section 3 we provide the bounds for each of the settings and estimands of interest.In Section 4 we conduct some numerical studies to gain insight into the bounds we provide.In Section 5 we illustrate our derived bounds in an illustrative data example from the mediation package in R. Code to reproduce all results in the simulations and real data example is available at [blinded-url] in addition to R functions to use the bounds in real data.Finally, in Section 6 we outline the limitations of our bounds and discussion future areas of research.

Notation and Setting
Let X and Y be the binary intervention and binary outcome of interest.Let M 1 and M 2 be two binary mediators on paths from X to Y .Let U be an unmeasured set of confounders between Y , M 1 and M 2 that are independent of X.The variables in this set have an unrestricted and unknown distribution, e.g. they may be a set of continuous and correlated unmeasured variables.We let Y (X = x) denote the potential outcome Y under the intervention which sets X to x, and Y (x, m 1 , m 2 ) be the equivalent under an intervention which sets X to x, M 1 to m 1 , and M 2 to m 2 .
Our setting of interest is as depicted in the causal model or directed acyclic graph (DAG) in Figure 1.We interpret the DAG in Figure 1 to represent a set of nonparametric structural equations such that: The set of 's represent 'errors terms' due to omitted factors, which are assumed independent of U and of each other.Given the values of the errors and the values of a variable's parents in the graph, the value of the variable is determined by its response function.The errors determine the manner in which the variable is determined from its parents.
In this setting, the total effect (TE) is identified and can be estimated.However, none of the mediation effects are identified from the observed data due to the unmeasured confounder(s) U , which open pathways between X and Y when conditioning on M 1 or M 2 .
Both Steen et al. [2017] and Daniel et al. [2015] provide rigorous proofs of this statement.Thus, in the following section we provide valid and sharp bounds for these estimands.
These bounds are functions of the observed data distribution p(Y, M 1 , M 2 |X).We say that the bounds are 'sharp' when each value within the bounds is a possible value for the estimand, given the true probabilities that are estimable from data.Similarly, we say that the bounds are 'valid' when no value outside the bounds is a possible value, given the true probabilities that are estimable from data.

Estimands
There are several other potential outcomes to define in the setting with two mediators when the mediators are controlled 'naturally'; what these estimands are depends on whether there is an effect of M 1 on M 2 .When there is an effect of M 1 on M 2 , as in our setting of interest, Figure 1, the potential outcome becomes relevant.This is the potential outcome of Y under the interventions that sets X to x, and M 1 to the value it would take on under the intervention that sets X to x 1 , and M 2 to the value it would take on under the intervention that sets X to x 2 and M 1 to the value it would take on if X was set to x 3 .If there is no effect of M 1 on M 2 , this potential outcome simplifies to In what follows, we focus on the setting of Figure 1 where we allow an effect of M 1 on M 2 .
The total effect of X on Y is defined as: In the randomized and perfect compliance setting, such as our setting of interest, the TE is identified, and equal to p(Y = 1|X = 1) − p(Y = 1|X = 0).We are also interested in the direct effect of X on Y holding the mediators to either a fixed level for all subjects, or the natural level they would have taken on had X been set to x. Holding a single mediator to a fixed level gives what is referred to as the controlled direct effect.This can directly be extended to multiple mediators.We define the controlled direct effects (CDE) for two mediators, which have four possible levels defined by the values to which the mediators are being held, by: This estimand fully describes the possible controlled direct effects in the setting of Figure 1 below.When you hold the mediators to the counterfactual value they would have taken had you intervened on X, this is referred to as the natural direct effect, for a single mediator.Similar to the multiple mediators in the controlled direct effect case, we define the natural direct effects (NDE), in the setting of Figure 1.These are the estimands discussed in Daniel et al. [2015], and here we follow their nomenclature directly.
The estimand NDE-000 was said to be the obvious extension of the pure natural direct effect in Daniel et al. [2015], and is the term in equation 7 in the decomposition of Steen et al. [2017].
We also consider the effect transmitted along either one or both mediators, the joint natural indirect effect: . JNIE 1 is equal to the term in equation 6 of the decomposition as given in Steen et al. [2017].
We can further decompose this joint effect into two indirect effects, where the first component is the effect through M 1 , directly and also through M 2 , in the notation of Daniel et al. [2015]: given in Steen et al. [2017].
Finally, we consider the indirect effect through M 2 , excluding the effect of M 1 through M 2 , as given in Daniel et al. [2015]: and is equivalent to We focus on the bounds for the terms of these two decompositions in the main text.
However, we also give the bounds on all terms in the twenty-four, four-term decompositions of Daniel et al. [2015] in Section S3 of the supplementary materials and the decompositions are in Section S2.For this we need to define two more estimands, where again we use the same notation as Daniel et al. [2015].
The indirect effect through M 1 , excluding the effect of and the indirect effect of M 1 through M 2 , excluding the direct effect of M 1 .
These terms add to the effect through M 1 and possibly also through M 2 , MS Define the short hand notation for the estimable probabilities as: For example,

Results
As the variables of interest in the DAG in Figure 1 are all assumed binary, there exists a canonical partitioning of the unmeasured confounder U into finite states, as described by Balke and Pearl [1994].In this partitioning, the response function corresponding to each variable in the DAG is categorical with 2 2 k levels, where k is the number of parents of that variable in the DAG.Ignoring the response function variable for X since we condition on it, this leads to a total of 2 2 1 × 2 2 2 × 2 2 3 = 16, 384 probabilities associated with the joint response function variable distribution for (M 1 , M 2 , Y ), on which the thirty-two estimable probabilities form constraints that are used to bound the estimands of interest.The mediation effects of interest are the objectives that we maximize and minimize symbolically using vertex enumeration, resulting in bounds on the counterfactual probabilities in terms of the estimable data distribution.
Result 1: The bounds given in (1) are valid and sharp for CDE-m 1 m 2 under Figure 1.
With a small amount of algebra one can write the bounds in Result 1 in terms of the total effect, TE, with the lower bound as TE −B(m 1 , m 2 ), and the upper as TE . From this representation it is easy to see that when g(m 1 , m 2 ) = 0, this also implies B(m 1 , m 2 ) = 0, and the bounds collapse to the single point, TE.Additionally, it can be seen from this that if TE > B(m 1 , m 2 ) then CDE-m 1 m 2 > 0. It is also immediately evident that when p(M 2 = m 2 , M 1 = m 1 ) = 1, then all NDE collapse to the CDE for the same m 1 , m 2 , and the same result applies.However, in this setting this is less mediation and more a deterministic relationship.
Result 2: The bounds given in ( 2) and ( 3) are valid and sharp for NDE-000 under Figure 1.
and NDE-000 Result 3: The bounds given in (4) and ( 5) are valid and sharp for JNIE 1 under Figure 1. and Result 4: The bounds given in ( 6) and ( 7) are valid and sharp for MS 2 -NIE 1 -11 under Figure 1. and To compare the bounds presented above, it is useful to consider alternative valid bounds defined by 'subtraction procedures' or 'addition procedures', as follows.Consider an estimand of interest θ such that it can be decomposed into several terms θ = I i=1 γ i .Let (l i , u i ) be the valid and sharp bounds for each γ i , respectively.Given the decomposition of θ we have γ i = θ − j =i γ j so that alternative bounds for γ i are given by (l * i , u * i ), where l * i = θ i − j =i u j , and u * i = θ i − j =i l j .This is what we will call the subtraction procedure.Similarly, one can obtain valid bounds for θ by adding the bounds for γ i u * = i u i , and This is what we will call the addition procedure.These procedures will not always produce sharp bounds, as is demonstrated in the real data example, but there are at least some special cases where the subtraction procedure can be used to produce valid and sharp bounds.

Result 6:
For any two-way decomposition of the identified TE the subtraction of the limits of sharp bounds for one term in the decomposition from the TE will produce the sharp bounds for the other term.
Result 6 is easily proven by considering three estimands, θ, γ 1 and γ 2 , such that θ = γ 1 + γ 2 where θ is identified, but γ 1 and γ 2 are not.Suppose that (l 1 , u 1 ) and (l 2 , u 2 ) are valid and sharp bounds for γ 1 and γ 2 , respectively.We know that γ 2 = θ − γ 1 and that the validity of the bounds for γ 1 implies that θ − u 1 ≤ γ 2 ≤ θ − l 1 is valid.In addition, the sharpness of the bounds for γ 1 imply that any point in θ − u 1 ≤ γ 2 ≤ θ − l 1 is possible, which then implies that θ − u 1 = l 2 and θ − l 1 = u 2 , if the bounds for γ 2 are also valid and sharp.This is also easy to see by looking at the first panel of Figure 4.If the bounds are sharp, then the triangles must be of the same size and mirror images of themselves over the line defined by the TE.If one of the triangles is smaller on one side, then by projecting across the TE line you could produce narrower bounds for the other estimand.This contradicts the conditions of the result.Alternatively, if one of the triangles is too large, projecting from the smaller triangle towards the term you are subtracting results in the same contradiction, thus proving the result.
Result 6 does not generalize to a three-way decomposition.For example, θ = γ 1 +γ 2 +γ 3 , where again θ is identified, but γ 1 , γ 2 and γ 3 are not.We now also have sharp and valid bounds for γ 3 (l 3 , u 3 ).We can again bound γ 2 by θ − u 1 − u 3 ≤ γ 2 ≤ θ − l 1 − l 3 .By the validity of the bounds for γ 1 and γ 3 this bound is valid for γ 2 .However, they are only sharp if the limiting values l 1 and l 3 are simultaneously possible values of γ 1 and γ 3 , and u 1 and u 3 are also simultaneously possible.If the bounds are not constructed by considering the γ's jointly, there is no guarantee that γ 1 can take on the value l 1 , while γ 3 is equal to l 3 , or γ 3 taken on the value u 3 while γ 1 is equal to u 1 .Even in the case of a two-way decomposition it is clear that l 1 + l 2 and u 1 + u 2 will not always provide sharp bounds for the identified θ, unless l 1 = u 1 and l 2 = u 2 .Additionally, as is demonstrated in our real data example, if θ is also not identified with θ = γ 1 + γ 2 and has valid and sharp bounds (l 0 , u 0 ), the bounds derived for γ 2 using the subtraction procedure, l 0 − u 1 ≤ γ 2 ≤ u 0 − l 1 , have no guarantee of being sharp or even informative, i.e. not containing -1 and 1.This makes it clear that the bounds for JNIE 1 are generally narrower than those resulting from adding the limits of the bounds for each term in its decomposition; this is demonstrated in the lower panel of Figure 4.In this example, the bounds for NIE 2 -100 are wider than those for JNIE 1 , making clear that the addition procedure will not produce sharp bounds.Although in our real data example bounds for MS 2 -NIE 1 -11 are narrower than JNIE 1 , this will not always be the case.The bounds for MS 2 -NIE 1 -11 differ by a single term in each of the upper and lower bounds (4) and ( 5) for JNIE 1 .These terms, when active, tend to make the bounds narrower, although this does not often mean they exclude zero, as we demonstrate in simulations below.

Simulations
We randomly generated counterfactual probabilities and then use the constraints implied by the DAG to generate the true estimable probabilities, allowing us to ensure that we are generating distributions under the DAG.We generated the K = 16,384 dimensional counterfactual probability distribution vector by sampling from a Dirichlet distribution with the vector of parameters α = {α 1 , . . ., α K }, that has probability distribution function: where B(•) is the multivariate beta function, and where i q i = 1.
The generated counterfactuals describe a 16,384-dimensional space that is difficult to explore in any exhaustive manner.We consider two special cases, points at the vertices of the space where a single counterfactual probability is 1 and all others are zero, and the symmetric Dirichlet distribution, where all α i are equal.Additionally we consider only two characteristics of the bounds under each of these special cases, the width, with anything less than 2 providing information over the always valid (-1,1) bounds for any risk difference, and if the bounds cover zero, thus not providing evidence against the causal null hypothesis.
When considering the vertices of the 16,384-dimensional space, we found that the width of the bounds under these distributions is bimodal/trimodal, see Table 1.Specifically, for the NDE-000 and JNIE 1 , the widths of the bounds are 0 for 25% of the vertices and 1 for 75% of the vertices.For the MS 2 -NIE 1 -11, the widths of the bounds are 0 for half of the vertices and 1 for the other half.In all those cases, all of these bounds either exclude 0 or 0 is the lower or upper limit of the bounds.For the NIE 2 -100, the widths of the bounds are 0 for 25%, 1 for 25%, and 2 for 50% of the vertices.In the cases where the widths of the bounds are not 2, they either exclude 0 or 0 is the lower or upper limit of the bounds.
When the α i take a common value decreasing from 1 to 0, the counterfactual space has more and more of a bowl shape, such that the areas where all the counterfactual probabilities would be non-zero and equal has the lowest probability of being generated and the vertices where only one probability is non-zero has the highest.Figures 2 and 3 trace the path from all the α i being nearly zero (0.000001) to all being equal to 0.001, for the width and exclusion of 0, respectively.As can be seen in the figures, over the 1,000  For all estimands other than NIE 2 -100, 100 minus this % provides the bounds with width of exactly 1, as the bounds for these estimands never exceeds 1.For NIE 2 -100 it can be seen that the bounds width is often larger than 1 and in many cases is width 2.
simulations at each level of the α i only at very low values do all of the estimands have a large proportion of bounds with widths less than 1.However, at no value of the α i do the widths exceed 1 for the direct effect, the joint indirect effect, or the indirect effect through M 1 .The bounds for the effect through M 2 only are quite wide, often being greater than 1, with a large proportion having a width of 2 at all α i = 0.001.Additionally, we see that in most cases the bounds contain zero for the indirect effects, even when, in the case of the MS 2 -NIE 1 -111, there are a large proportion of the bounds that have width less than 1.The bounds for NDE-000, on the other hand, often exclude zero for low α i values.
This suggests that when many of the counterfactual probabilities are zero, the bounds are more informative.This is not surprising, as counterfactual probabilities being zero can be considered constraints, such as the well-known constraint in the instrumental variable setting called no defiers or monotonicity.In a 16,384-dimensional space, such constraints are more difficult to give intuitive names.When the α i are equal and held to be 1 this is the same as putting a flat uniform distribution on the probabilities.The trend in Figure 2 continues, with estimated probability of the width of the bounds being less than one approaching zero as the α i approaches one.

Real Data Example
We illustrate our bounds using the jobs dataset from the mediation package in R [Tingley et al., 2014].These data come from a randomized experiment designed to investigate the efficacy of a job training intervention on unemployed workers.In the experiment, 899 eligible unemployed workers completed a pre-screening questionnaire and were then randomly assigned to treatment, which consisted of participation in job skills workshops, or control, who received a booklet of job search tips.The randomization was 2:1 in favor of treatment with 600 assigned to treatment and 299 to control.The primary outcome is a binary variable representing whether the respondent had become employed, with 35% on the treatment arm and 29% on the control arm becoming employed at the end of the study.
The mediators of interest are a binary indicator of depressive symptoms after treatment, and a binary indicator of high job seeking self-efficacy.We believe that depression is the first mediator, or M 1 , which may have a causal effect on job seeking self-efficacy, M 2 .The bounds are displayed in Figure 4.
As can be seen in the first panel of Figure 4, the total effect was estimated to be barely above zero, 0.057, which by an exact test is not significantly different than 0, with a p-value of 0.10 and 95% CI (−0.01, 0.12).Also in this panel it can be seen that the bounds for the natural direct effect (−0.29, 0.71) and the joint indirect effect (−0.65, 0.34) are mirror images across the line of the total effect.This is as expected, given Result 6.Looking at the second panel in Figure 4 it can be seen that the bounds for the MS 2 -NIE 1 -111 (−0.50, 0.34) are narrower than both the natural indirect effect through high job seeking (−0.96,0.84) and the joint indirect effect.It can also be seen that the addition procedure used on the bounds for MS 2 -NIE 1 -111 and NIE 2 -100 is wider than the sharp bounds for JNIE 1 , as the bounds for NIE 2 -100 alone are wider than those for JNIE 1 .The subtraction procedure of the lower bound of NIE 2 -100 from the upper bound of the JNIE 1 , as well as the upper from the lower, produce much wider bounds for MS 2 -NIE 1 -111.This is, again, as expected as the bounds are sharp.In this case, the addition procedure and subtraction procedure both produce noninformative bounds, i.e. the bounds are outside the possible range of the causal effect (−1, 1).Additionally, although it was not the focus of this paper, Result 1 provides the bounds for the controlled direct effects and they are as follows in these data: CDE-00 ∈ (−0.71, 0.79), CDE-01 ∈ (−0.51, 0.54), CDE-10 ∈ (−0.84, 0.86), and CDE-11 ∈ (−0.87, 0.87).Each CDE bound covers zero and has a width greater than one, with CDE-01 having the narrowest width of 1.03.

Discussion
We present bounds for estimands of mediation effects in the two-mediator setting, including many of the estimands considered in Daniel et al. [2015], Steen et al. [2017].We show that in many cases the bounds will be narrower than one, and this occurs frequently when a large number of the counterfactual probabilities are zero.We also find that the addition or subtraction of sharp bounds for two estimands that are not identified does not produce, in general, sharp bounds for the estimand of their addition or subtraction.Additionally, the subtraction of sharp bounds does not produce sharp bounds for the remaining term(s) unless one of the estimands in the subtraction is identified.We prove that in a two-way decomposition of the identified TE, the subtraction of a set of sharp bounds will always produce sharp bounds for the other term.Daniel et al. [2015] showed that many of the mediation estimands are not identified even in the complete absence of unmeasured confounders, making bounds of particular relevance.
However, it should be noted that the bounds we provide are not generated under the assumption of no confounding.So, while the bounds will be valid under the assumptions of no confounding, one could possibly obtain narrower sharp bounds by including that assumption in the derivation.
Although we explore two specific scenarios in the simulations, exploring the full space, perhaps with some type of greedy algorithm, to determine if there is a pattern which corresponds to constraints is beyond the scope of this paper, but an area of future research for the authors.Additionally, in exploring the space in greater detail, we will likely find constraints to consider in the construction of bounds, which may result in narrower or different sets of bounds for these decomposition terms.
We provide a large number of bounds, and insight into their interpretation, but we do not discuss an estimation procedure.As the terms of the bounds are all linear combinations of conditional probabilities, estimation is straightforward.For inference, the same bootstrap procedure suggested in Gabriel et al. [2020] andHorowitz andManski [2000] can be used here.In both of those papers, extensive empirical results were provided showing nominal coverage of quantile bootstrap confidence intervals for similarly constructed bounds.Although we believe that the bootstrap is likely theoretically justified, there may be closed form variance estimates for many, if not all, of the bounds.We do not focus on accounting for the sampling variability here, as for any reasonable sample size, we expect the uncertainty in causal effects due to unmeasured confounding to be far greater than uncertainty due to sampling variability.
Although we only consider the setting where there are unmeasured and no measured confounders, in the simplest case where the measured confounders are discrete and bounded and do not restrict the impact of the unmeasured confounders on the mediators and the outcome, the bounds can be applied within levels of the measured confounders.More complex settings would need to be considered on a case-by-case basis and may result in narrower sharp bounds under particular assumptions.There are other scenarios where the same linear programming method may be used to determine if the bounds change, such as allowing the intervention to have a direct effect on the confounder.Defining the general conditions under which a DAG and target define a linear programming problem, and the automated checking of these conditions is an area of future research for the authors.Finally, our present results are limited to binary exposures, mediators, and outcomes so it would be worthwhile to extend these to multicategorical variables.

Figure 1 :
Figure 1: Causal diagram of the setting of interest.
100 is equal to the term in equation 9 of the decomposition given inSteen et al. [2017].Then the decomposition inSteen et al. [2017] is in our notation TE = NDE-000 + JNIE 1

Figure 2 :
Figure2: Width of bounds (points) with boxplots for 1,000 simulated distributions generated under different values of the parameter of a symmetric Dirichlet distribution.The orange numbers at the top of each figure indicate the % of bounds that have width less than 1.For all estimands other than NIE 2 -100, 100 minus this % provides the bounds with width of exactly 1, as the bounds for these estimands never exceeds 1.For NIE 2 -100 it can be seen that the bounds width is often larger than 1 and in many cases is width 2.

Figure 3 :
Figure 3: Proportion of bounds for the NDE-000 that exclude 0 out of 1,000 simulated distributions generated under different values of the parameter of a symmetric Dirichlet distribution.Grey ribbon indicates Clopper-Pearson 95% confidence limits for the proportion.

Figure 4 :
Figure 4: Bounds for the jobs dataset.The orange arrows show the addition procedure, summing the upper and lower bounds values of MSE 2 -NIE 1 -11 and NIE 2 -100 to obtain valid, but not sharp bounds for NIE 1 .The brown arrows demonstrate the subtraction procedure where an alternative lower bound for MSE 2 -NIE 1 -11 is derived by subtracting the upper bound value of NIE 2 -100 from the lower bound value of JNIE 1 , and an alternative upper bound for MSE 2 -NIE 1 -11 is derived by subtracting the lower bound value of NIE 2 -100 from the upper bound value of JNIE 1 .Again, this is valid but not sharp.

Table 1 :
Counts (proportion) out of the 16,384 vertices where the lower bound equals the label in the rows, and the upper bound equals the label in the columns.
Our simulations do not give a complete picture of what is occurring in the full 16,384dimensional space, as we have also found examples of sections of the space where two counterfactuals are both nonzero and the rest are equal but nonzero, where both the width of the bounds are less than 1 and the bounds exclude zero for the NDE-000, JNIE 1 and MS 2 -NIE 1 -11 effects.Given our limited exploration we have not uncovered a discernible pattern, other than the bounds tend to be more informative when a large number of the counterfactuals are zero.
2-NIE 1 -11 and NIE 2 -100 to obtain valid, but not sharp bounds for NIE 1 .The brown arrows demonstrate the subtraction procedure where an alternative lower bound for MSE 2 -NIE 1 -11 is derived by subtracting the upper bound value of NIE 2 -100 from the lower bound value of JNIE 1 , and an alternative upper bound for MSE 2 -NIE 1 -11 is derived by subtracting the lower bound value of NIE 2 -100 from the upper bound value of JNIE 1 .Again, this is valid but not sharp.