Identification and Estimation of Joint Probabilities of Potential Outcomes in Observational Studies with Covariate Information
Year: 2022

Abstract:
The joint probabilities of potential outcomes are fundamental components of causal inference in the sense that (i) if they are identifiable, then the causal risk is also identifiable, but not vise versa (Pearl, 2009; Tian and Pearl, 2000) and (ii) they enable us to evaluate the probabilistic aspects of “necessity”, “sufficiency”, and “necessity and sufficiency”, which are important concepts of successful explanation (Watson, et al., 2020). However, because they are not identifiable without any assumptions, various assumptions have been utilized to evaluate the joint probabilities of potential outcomes, e.g., the assumption of monotonicity (Pearl, 2009; Tian and Pearl, 2000), the independence between potential outcomes (Robins and Richardson, 2011), the condition of gain equality (Li and Pearl, 2019), and the specific functional relationships between cause and effect (Pearl, 2009). Unlike existing identification conditions, in order to evaluate the joint probabilities of potential outcomes without such assumptions, this paper proposes two types of novel identification conditions using covariate information. In addition, when the joint probabilities of potential outcomes are identifiable through the proposed conditions, the estimation problem of the joint probabilities of potential outcomes reduces to that of singular models and thus they can not be evaluated by standard statistical estimation methods. To solve the problem, this paper proposes a new statistical estimation method based on the augmented Lagrangian method and shows the asymptotic normality of the proposed estimators. Given space constraints, the proofs, the details on the statistical estimation method, some numerical experiments, and the case study are provided in the supplementary material.

heading: 1 Introduction

heading: 1.1 Practical Background
text: In practical science, it is crucial to evaluate the likelihood of one event causing another event. For example, epidemiologists pay attention to determining the likelihood of a particular exposure being the cause of a particular disease. In order to assess such likelihood from statistical data, the probabilities of causation have been developed, which can be divided into “necessary causation”, “sufficient causation”, and “necessary-and-sufficient causation”. In addition, recently, in the field of explainable artificial intelligence (XAI), it has been pointed out that “necessity”, “sufficiency”, and “necessity and suffciency” are the important concepts of successful explanation and the probabilities of causation play an important role
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
in evaluating these concepts from the probabilistic aspects [10, 22, 35]. For the last decades, the probabilities of causation received increasing attention in the fields of epidemiology, risk analysis, legal reasoning, artificial intelligence, and policy analysis [e.g., 4–8, 10, 11, 22, 26, 29, 33, 35].

heading: 1.2 “Necessity”, “Sufficiency”, and “Necessity and Sufficiency”
text: Necessary cause is often used in the tort-law liability. When a 60-year-old man who worked in an asbestos textile plant develops lung cancer, we are interested in the probability that asbestos had caused it. This quantity to be estimated has been commonly referred to as “attributable fraction” or “assigned share” [4, 19]. It measures how necessary the exposure is for the development of the disease. Pearl [24] defined it as “Probability of Necessity (PN)”, which means “the probability that disease would not have occurred in the absence of exposure, given that disease and exposure did in fact occur”. The estimation of PN is important to provide a scientific basis for compensation in tort suits. Sufficient cause is often used in health risk assessment. In practice, it is important for public policy makers to predict how a health risk factor would endanger a healthy population [15, 16, 21]. The measure for it has been referred to as “susceptibility”, which measures how sufficient a risk factor is to make an individual develop a disease. Pearl [24] defined the measure of such prediction as “Probability of Sufficiency (PS)”, which means “the probability that a healthy unexposed individual would have developed the disease had he been exposed”. The assessment of PS helps public policy makers to make or revise regulations in order to minimize the occurrence of diseases. With the definitions of PN and PS, it would be a natural extension to inquire for the “Probability of Necessity and Sufficiency (PNS)”, which shows how likely an individual is to be affected by both ways [24]. PNS measures the probability that an exposure is a necessary and sufficient cause for a disease, that is, an exposure is the actual cause of a disease. Then, different from previous works where “probability of causation” is only referred to PN [e.g., 4, 7], Pearl [24] defined three types of “probability of causation”, PN, PS, and PNS. Pearl [26] and Tian and Pearl [33] developed formal semantics for the probabilities of causation based on the structural models of counterfactuals. They presented the formal definitions of the probability of necessity (PN), the probability of sufficiency (PS), and the probability of necessity and sufficiency (PNS). These probabilities of causation are formulated based on the joint probabilities of two potential outcomes. Since one cannot simultaneously observe the results of the same subjects exposed and unexposed in reality, even successful randomized experiments have failed to identify these quantities [26, pp.284–285]. Tian and Pearl [33] addressed the problem by showing how to bound these quantities using data obtained from experimental and observational studies. Their bounds are sharp under minimal assumptions concerning the data generating process. Kuroki and Cai [17] derived narrower bounds of the probabilities of causation than Tian-Pearl’s bounds using covariate information. Recently, Galhotra et al. [10] generalized Tian-Pearl’s bounds from binary to arbitrary variables in the context of explainable artificial intelligence (XAI). However, it has been often pointed out that these existing bounds are too wide to evaluate the probabilities of causation. To solve the problem, Tian and Pearl [33] and Galhotra et al. [10] also stated that the probabilities of causation are identifiable if the monotonicity can be assumed and the causal risks are identifiable and Pearl [26] showed that the specific functional relationships between cause and effect lead to the identification of the probabilities of causation.

heading: 1.3 Contribution
text: The aim of this paper is to discuss the identification and estimation problems of the joint probabilities of potential outcomes, which provides a general class including the probabilities of causation. The joint probabilities of potential outcomes are the building-blocks of causal inference, because the causal risk is also identifiable if they are identifiable, but not vise versa. In the context of the natural direct and indirect effects [25], under the assumption of no unmeasured confounding, Robins and Richardson [30] stated that the joint probabilities of potential outcomes are identifiable if (i) two potential outcomes are independent or (ii) one potential outcome is deterministically formulated as the function of the other potential
outcome. In addition, in the context of unit selection problems, Li and Pearl [20] showed that the linear combination of the joint probabilities of potential outcomes is identifiable under the gain equality. These existing researches show that the evaluation of the joint probabilities of potential outcomes play an important role in solving various problems of causal inference. Although considerable efforts have been devoted toward establishing the identifiability criteria and solving the confounding problem of causal risks, there has been limited discussion on how to identify the joint probabilities of potential outcomes when the present assumptions are violated. To solve the problem, this paper proposes two types of novel identification conditions for the joint probabilities of potential outcomes without existing assumptions — one uses both causal risks and one covariate (Theorem 1) and the other uses two covariates without causal risks (Theorem 2). The results of this paper also provide a counterexample to “sensitivity to the generative process” [26, pp. 284–285], in the sense that it is not required to specify the functional relationships that connect a cause and an effect. In addition, unlike the effect restoration proposed by Kuroki and Pearl [18], the assumption of the given number of categories of an unmeasured confounder is not required. Furthermore, from Theorem 2, it is remarkable that there are some situations where the causal risks are identifiable even if the existing identification conditions do not hold. When the joint probabilities of potential outcomes are identifiable through the proposed conditions, the estimation problem of the joint probabilities of potential outcomes reduces to that of singular models and thus they can not be evaluated by standard statistical estimation methods. To solve the problem, this paper proposes a new statistical method based on the augmented Lagrangian method [3] and shows the asymptotic normality of the proposed estimators. Given space constraints, the proofs, the details on the statistical estimation method, some numerical experiments, and the case study are provided in the supplementary material.

heading: 2 Preliminary
text: This section introduces the potential outcomes used to discuss our problems. For simplicity, we consider the case where an exposure variable X ∈ {x1, x0} (x1 is exposed; x0 is unexposed) and an outcome variable Y ∈ {y1, y0} (y1 is disease; y0 is non-disease) are dichotomous. In addition, let p(X = x, Y = y) = p(x, y), p(Y = y |X = x) = p(y |x), and p(X = x) = p(x) be the joint probability of (X,Y ) = (x, y), the conditional probability of Y = y given X = x, and the marginal probability of X = x for x ∈ {x1, x0} and y ∈ {y1, y0}, respectively. Similar notation is used for other probabilities. Here, in this paper, we will assume that readers are familiar with the basic theory of causal inference [13, 26]. Especially, for the graph-theoretic terminology and the basic theory of structural causal models used in this paper, we refer the readers to Pearl [26]. In principle, for x ∈ {x1, x0}, the i-th of the N subjects has a potential outcome Yx(i) that would have resulted if X had been x, denoted as X(i) = x. Here, note that the subject ensures a deterministic relationship between two variables X and Y in the semantics of structural causal models [26]. In addition, the present paper assumes the stable unit treatment value assumption, which can be summarized as follows: (i) the exposure status of any subject does not affect the outcomes of the other subjects (no interference) and (ii) the exposures of all subjects are comparable (no variation in exposure). Thus, when the N subjects in the study are considered as random samples from the population of interest, X(i) and Yx(i) are referred to as the values of random variables X and Yx respectively, and thus the causal risk of X = x on Y = y is defined as p(Yx = y). Similar notation is used for other potential outcomes. For the i-th subject, the potential outcome Yx is observed only if X is x . This property is called the consistency [13, 26, 28], which is formulated as X(i) = x=⇒Yx(i) = Y (i). When a randomized experiment is conducted and compliance is perfect, since X is independent of (Yx0 , Yx1), the causal risk is identifiable and is given by p(Yx = y) = p(y |x). Here, “identifiable” means that the causal quantities, such as p(Yx = y), can be estimated consistently from a joint probability of observed variables. In contrast, when it is difficult to conduct
a randomized experiment and only observational data are available, we can evaluate the causal risk according to the conditionally-ignorable-treatment-assignment condition [31], or graphically, the back-door criterion [26]. In other words, for the exposure X, if there exists such a set S of observed covariates that X is conditionally independent of (Yx0 , Yx1) given S, we can say that treatment assignment is conditionally ignorable given S. Then, the causal effects are estimable using S as p(Yx = y) = Es[p(y |x,S)]. Here, Es[p(y |x,S)] is the expectation of p(y |x,S) regarding S. Although there are other identification conditions that can be used to solve our problem [e.g., 26, 34], this paper does not cover them due to space constraints. Following Pearl [26], we define three probabilities of causation, that is, the probability of necessity (PN), the probability of sufficiency (PS), and the probability of necessity and sufficiency (PNS). PN stands for the probability p(Yx0 = y0 |X = x1, Y = y1) that a disease would not have occurred in the absence of an exposure, given that the exposure and the disease occurred. PS stands for the probability p(Yx1 = y1 |X = x0, Y = y0) that a healthy unexposed subject would have contracted the disease upon exposure. PNS stands for the probability p(Yx1 = y1, Yx0 = y0) that an exposure is a necessary and sufficient cause for a disease. Since PN, PS, and PNS involve the joint probabilities of potential outcomes, even under successful randomized experiments, they cannot be estimated from the observed data without additional assumptions [26].

heading: 3 Identification

heading: 3.1 Case 1: Identification based on one covariate with causal risks
text: According to Pearl [26], letting U be the set of all discrete and continuous covariates that could affect X and Y , both observed and unobserved, we consider the problem of evaluating the joint probabilities of potential outcomes based on the directed graph shown in Figure 1a. Here, a covariate Z, in Figure 1a, is measured as a proxy variable of U . Note that Z can be a set of discrete and/or continuous variables. In Figure 1a, the directed edge from X to Y indicates that X may affect Y . In addition, the absence of a directed edge from Y to X indicates that Y cannot be a cause of X, and the directed path from U to Y through X indicates that some elements of U may effect Y mediated by X. Figure 1a also graphically represents the data generating process
Y = gy(X,U , y), X = gx(U , x), Z = gz(U , z), (1) where x, y, and y are independent random disturbances and they are also independent of U . When the structural equation models, such as equation (1), are used to represent the data generating process, the corresponding graph shown in Figure 1a is called a causal diagram. The situation such as Figure 1a is also discussed in Kuroki and Pearl [18] and Pearl [27]. However, in this paper, unlike Kuroki and Pearl [18] and Pearl [27], U can include uncertain number of all discrete and continuous covariates that influence the way a subject responds to exposures. Thus, in many situations, it is reasonable to assume the existence of a covariate Z that is independent of {X,Y } given U . Here, note that the independence assumptions between two observed variables would be affected by partitioning the states of U ∪{ x, y z}. In the situation shown in Figure 1a, irrespective of the complexity of U ∪ { x, y, z}, the impact of U on Y remains restricted to the modification of the functional relationships between X and Y . This yields four functions for two dichotomous variables X and Y , and thus the value taken by U ∪ { x, y, z} selects one of these four functions [26]. Considering these observations, according to Rothman et al. [32, p. 59], the states of U ∪ { x, y, z} are divided into the following four potential outcome types: u1 = (Yx1 = y1, Yx0 = y1) represents a “doomed” situation where exposure is irrelevant because disease occurs with or without exposure.
u2 = (Yx1 = y1, Yx0 = y0) represents a “causative” situation where disease occurs if and only if subjects are exposed.
u3 = (Yx1 = y0, Yx0 = y1) represents a “preventive” situation where disease occurs if and only if subjects are unexposed.
u4 = (Yx1 = y0, Yx0 = y0) represents an “immune” situation where exposure is irrelevant because disease does not occur, with or without exposure.
According to this partition of the states of U ∪ { x, y, z}, it is re-defined as a variable U taking a value u (u ∈ {u1, u2, u3, u4}). Then, similar to the setting of the non-compliance problem by Pearl [26, ch.8], throughout the paper, we assume that the recursive factorization of the positive joint probability of X, Y , and Z, p(x, y, z), is given as
p(x, y, z) = 4∑
i=1 p(y |x, z, ui)p(x |ui)p(z | ui)p(ui)
according to Figure 1b. It is also assumed that other joint probabilities of X, Y , and Z can be rewritten in a similar way. Then, this factorization implies that X is conditionally independent of Z, given U , denoted as X á Z | U , which is used in the context of the impact evaluation problem of the program [23]. In addition, from the consistency property, we have
p(x1, y1, z, ui) = p(x1, z, ui), p(x1, y0, z, ui) = 0, i = 1, 2, p(x1, y0, z, ui) = p(x1, z, ui), p(x1, y1, z, ui) = 0, i = 3, 4, p(x0, y1, z, ui) = p(x0, z, ui), p(x0, y0, z, ui) = 0, i = 1, 3, p(x0, y0, z, ui) = p(x0, z, ui), p(x0, y1, z, ui) = 0, i = 2, 4,
Thus, for example, p(x1, y1, z), p(z) and p(Yx1 = y1, z) can be rewritten as
p(x1, y1, z) = 2∑
i=1 p(x1 |ui)p(z |ui)p(ui), p(z) = 4∑ i=1 p(z |ui)p(ui)
p(Yx1 = y1, z) = p(Yx1 = y1 | z)p(z).
Additionally, the other joint and marginal probabilities of X, Y , and Z can be rewritten in a similar way. Then, when Z is a variable with the number of values k ≥ 4, say z1, . . . , z4, letting
P =  1 p(z1) p(z2) p(z3)p(x1, y1) p(x1, y1, z1) p(x1, y1, z2) p(x1, y1, z3)p(x1, y0) p(x1, y0, z1) p(x1, y0, z2) p(x1, y0, z3) p(x0, y1) p(x0, y1, z1) p(x0, y1, z2) p(x0, y1, z3)  , (2) Q =
( 1 p(z1) p(z2) p(z3) p(Yx1 = y1) p(Yx1 = y1, z1) p(Yx1 = y1, z2) p(Yx1 = y1, z3) p(Yx0 = y1) p(Yx0 = y1, z1) p(Yx0 = y1, z2) p(Yx0 = y1, z3) ) , (3)
we derive the following theorem:
Theorem 1 Letting Z be a variable with the number of values k ≥ 4, say z1, . . . , z4, p(Yx1 = y, Yx0 = y′) and p(Yx1 = y, Yx0 = y′ |x) are identifiable (y, y′ ∈ {y1, y0};x ∈ {x1, x0}) if the following conditions are satisfied:
Condition 1 p(Yx1 = y1 | z), p(Yx0 = y1 | z), p(x1, y1, z), p(x0, y1, z), p(x1, y0, z), and p(x0, y0, z) are available for z ∈ {z1, . . . , z4}.
Condition 2 {X,Y } á Z | U holds.
Condition 3 P is invertible and QP−1 consists of four distinct non-zero column vectors.
The proof is given in the supplementary material. Theorem 1 shows that PN, PS and PNS are identifiable. In addition, the causal risk is also identifiable from Theorem 1. Here, we would like to state some remarks. First, the violation of the assumptions can be testable: as seen from the supplementary material, if the simultaneous linear equation (A.12) in the supplementary material has no solutions, then at least one of the conditions does not hold. Second, although Condition 1 seems to require that some distributions are available, it is not the case. Theorem 1 is applicable with a single observed distribution, say, p(x, y, z, w) if both p(Yx1 = y | z) and p(Yx0 = y | z) are identifiable from the distribution. For example, regarding p(x, y, z, w), when {Z,W} satisfies the back-door criterion relative to (X,Y ), p(Yx = y | z) is identifiable and is given by
p(Yx = y | z) = ∑
w
p(y |x, z, w)p(w | z). (4)
Second, Condition 3 implies that all potential outcome types exist and thus the monotonicity assumption (e.g. no-prevention assumption; p(u3) = 0) does not hold. In addition, it implies that X is not independent of U , and thus X is not randomly assigned. Thus, in Theorem 1, the causal risks are assumed to be estimable by the identification conditions in observational studies (such as the back- and front-door criteria [26]) or by external studies (e.g., pilot studies or scientific judgement). If we have prior knowledge that a certain potential outcome type does not exist in the situation discussed this section, that is, the monotonicity assumption, then the PN, PS, and PNS would be identifiable without the covariate Z [33]. Theorem 1 demonstrates that, when we observe one covariate associated with potential outcomes and the causal risks are available, the joint probabilities of potential outcomes can be constructed from the covariate and causal risks; this can render the joint probabilities of potential outcomes identifiable without the monotonicity assumption.

heading: 3.2 Case 2: Identification based on two covariates without causal risks
text: We will tackle the more difficult problem of estimating the joint probabilities of potential outcomes without a prior knowledge of causal risks. We will show that, under certain conditions, the joint probabilities of potential outcomes can be evaluated using only covariates. To see this, consider a causal diagram shown in Figure 2, which is obtained by adding an observed variable W to Figure 1b. Then, when we consider the positive joint probability of X, Y , Z, and W , from the consistency property, we have
p(x1, y1, z, w, ui) = p(x1, z, w, ui), p(x1, y0, z, w, ui) = 0, i = 1, 2, p(x1, y0, z, w, ui) = p(x1, z, w, ui), p(x1, y1, z, w, ui) = 0, i = 3, 4, p(x0, y1, z, w, ui) = p(x0, z, w, ui), p(x0, y0, z, w, ui) = 0, i = 1, 3, p(x0, y0, z, w, ui) = p(x0, z, w, ui), p(x0, y1, z, w, ui) = 0, i = 2, 4,
for z and w taken by Z andW , respectively. Thus, for example, p(y1, z, w |x1) and p(z, w |x1) can be rewritten as
p(y1, z, w|x1) = 2∑
i=1 p(w|x1, ui)p(z|ui)p(ui|x1), p(z, w|x1) = 4∑ i=1 p(w|x1, ui)p(z|ui)p(ui|x1).
Additionally, the other joint and marginal probabilities of X, Y , Z, and W can be rewritten in a similar way. Then, when Z and W are variables with the number of values k ≥ 4, such as z1, . . . , z4, for Z, and w1, . . . , w4 for W , letting
Px1 =  1 p(z1 |x1) p(z2 |x1) p(z3 |x1)p(w1 |x1) p(z1, w1 |x1) p(z2, w1 |x1) p(z3, w1 |x1)p(w2 |x1) p(z1, w2 |x1) p(z2, w2 |x1) p(z3, w2 |x1) p(w3 |x1) p(z1, w3 |x1) p(z2, w3 |x1) p(z3, w3 |x1)  , (5)
Px0 =  1 p(z1 |x0) p(z2 |x0) p(z3 |x0)p(w1 |x0) p(z1, w1 |x0) p(z2, w1 |x0) p(z3, w1 |x0)p(w2 |x0) p(z1, w2 |x0) p(z2, w2 |x0) p(z3, w2 |x0) p(w3 |x0) p(z1, w3 |x0) p(z2, w3 |x0) p(z3, w3 |x0)  , (6) we derive the following theorem:
Theorem 2 Letting Z and W be variables by taking k (≥ 4) values, say z1, . . . , z4, for Z, and w1, . . . , w4 for W , p(Yx1 = y, Yx0 = y′) and p(Yx1 = y, Yx0 = y′ |x) are identifiable (y, y′ ∈ {y1, y0};x ∈ {x1, x0}) if the following conditions are satisfied:
Condition 4 p(x1, y1, z, w), p(x1, y0, z, w), p(x0, y1, z, w), and p(x0, y0, z, w), are available for z ∈ {z1, . . . , z4} and w ∈ {w1, . . . , w4}.
Condition 5 Both W á Z | {X,U} and X á Z | U hold.
Condition 6 Both Px1 and Px0 are invertible.
The proof is given in the supplementary material. Although the probabilities comprising the unobserved variables are not fully identified in Kuroki and Pearl [18], Theorem 2 shows that the joint probabilities of potential outcomes p(u1), . . . , p(u4) are identifiable and the causal risks are also identifiable and are given by, for example,
p(Yx1 = y1) = p(u1) + p(u2), p(Yx0 = y1) = p(u1) + p(u3). (7)
Clearly, Theorem 2 shows that PN, PS and PNS are also identifiable.

heading: 4 Estimation
text: When the joint probabilities of potential outcomes are identifiable through the proposed conditions, as seen from the proofs of Theorems (refer to the supplementary material A), the estimation problem of the joint probabilities of potential outcomes reduces to that of singular models and thus they can not be evaluated by standard statistical estimation methods. To solve the problem, we proposes a new statistical method based on the augmented Lagrangian method [3]. In this section, we focus on Case 2 in Section 3.2. For estimation in Case 1, see the supplementary material B.1.
Consider the matrices P̂x and Q̂x that are derived by replacing p(z |x), p(w |x), p(z, w |x), p(y |x), p(y, z |x), p(y, w |x), and p(y, z, w |x) of Px and Qx with sample probabilities p̂(z |x), p̂(w |x), p̂(z, w |x), p̂(y |x), p̂(y, z |x), p̂(y, w |x), and p̂(y, z, w |x), respectively (x ∈ {x1, x0}; y ∈ {y1, y0}; z ∈ {z1, . . . , z4}; w ∈ {w1, . . . , w4}). Then, for constant matrices Mx1 and Mx0 (see equation (A.23) of the supplementary material A.2), since we have
SP−1x Qx = MxS, x ∈ {x1, x0} (8)
from equation (A.26) of the supplementary material A.2, we estimate the matrix S (see equation (A.24) of the supplementary material A.2) as the solution Θ = (θij)1≤i,j≤4 of the optimization problem
minimize Θ 1 2‖ΘP̂ −1 x1 Q̂x1 −Mx1Θ‖ 2 F + 1 2‖ΘP̂ −1 x0 Q̂x0 −Mx0Θ‖ 2 F (9)
subject to 0 ≤ e>j (P̂xΘ−1)>e1 ≤ 1, j = 1, . . . , 4; 1>(P̂xΘ−1)>e1 = 1, x ∈ {x1, x0}, (10)
û = p̂(x1)(P̂x1 Ŝ−1)>e1 + p̂(x0)(P̂x0 Ŝ−1)>e1. where p̂(x) is the sample probability for p(x). For details on the asymptotic normality of the estimators and algorithms, see the supplementary material.

heading: 5 Numerical Experiment
text: In this section, we present a numerical experiment to examine properties of the proposed estimation method through the PNS p(u2) and the causal risk difference p(u2)− p(u3). For simplicity, letting X, Y , Z, W , and U be discrete variables, we consider the causal diagrams shown in Figure 2, where the joint probabilities of (X,Y, Z,W,U) are given according to Table 1. Under the situation where (X,Y, Z,W ) can be observed but U can not, the properties of the proposed estimators p̂(u2) and p̂(u2)− p̂(u3) of p(u2) and p(u2)− p(u3), respectively, are verified in a numerical experiment using the setting with sample sizes n = 100, 200, 1000, and 5000. In this situation, since the p(u2) and p(u2)−p(u3) are 1/4 and zero, respectively, the sample means of p̂(u2) and p̂(u2)− p̂(u3) are expected to be close to 1/4 and zero, respectively. Table 2 and Figure 3 show the basic statistics and the box plots of p̂(u2) and p̂(u2)− p̂(u3) for 1000 replications with the given sample size n, respectively. The horizontal lines in Figure 3 show the true values of p(u2) and p(u2)− p(u3). From Table 2, the sample means of p̂(u2) and p̂(u2)− p̂(u3) are close to the true values and the sample standard deviations are smaller as the sample size is larger. Thus, it seems that the proposed estimation method provides the consistent estimators of p(u2) and p(u2)−p(u3). In addition, from Figure 3, the interquantile ranges for p̂(u2) and p̂(u2)− p̂(u3) are narrower and still include the true values even if the sample size is large. In contrast, it seems that p̂(u2) − p̂(u3) is symmetrically distributed and thus the asymptotic normality holds, but p̂(u2) may not. This is because the true value of p(u2) is relatively close to zero and Figure 3 is given as the zero truncated distribution of p̂(u2) for the finite sample size. Here, note that there are many outliers in each sample size. The outliers would occur when it is difficult to judge that Condition 6 holds from observed data. For further discussion on the simulation experiments, see the supplementary material.

heading: 6 Discussion
text: Joint probabilities of potential outcomes often appears in unit selection problems [20], the impact evaluation problem of social programs [12], the non-compliance problem of treatment effects [1, 2], the principal stratification [9], the identification problems of natural direct and indirect effects [25] and the prevented and preventable proportions [36], and the explainability problem of artificial intelligence [10, 22, 35]. In addition, if the joint probabilities of potential outcomes are identifiable, then the causal risks are also identifiable, but not vise versa. Therefore, the identification and estimation problems of these probabilities have been an important topic in causal inference. To solve the problems, the paper proposed (i) two types of novel identification conditions using covariate information and (ii) a new statistical estimation method based on the augmented Lagrangian method [3]. To the best of our knowledge, there has been previously much less discussion how useful covariate information is to identify these probabilities in nonparametric systems when existing identification conditions do not hold. In addition, the estimation problem of singular models often appears in the context of causal inference, but it seems that the solution has not been given in the literature. Thus, the results of this paper extend the range of solvable identification and estimation problems of causal inference under the nonparametric causal models. Finally, we assumed that both an exposure variable and an outcome variable is dichotomous. However, it is straightforward to extend our results from the case of dichotomous observed variables to the case of multivalued observed variables under certain assumptions. For example, multivalued or continuous outcome can be accommodated in the model using the event Y < y as an outcome variable. Then, for given y, the problem of the paper can be the identification and estimation of p(Yx > y, Yx′ ≤ y), . . . , p(Yx > y, Yx′ > y) and thus the result of this paper are applicable to the problems. For the related discussion, refer to Galhotra et al. [10] and Kada et al. [14]. In addition, when the exposure variable is continuous, according to Balke and Pearl [2], it is reasonable to assume that there exists an exposure interval around each x, within which a subject’s outcome is homogeneous. Under this assumption, it is possible to apply our results. However, in such cases, it may be difficult to obtain reliable statistics of the recovered probabilities due to data sparseness. It is left for future work.
Acknowledgments and Disclosure of Funding
We would like to acknowledge the helpful comments of the four anonymous reviewers. This research was partially funded by JFE Engineering Corporation and Japan Society for the Promotion of Science (JSPS), Grant Number 19K11856 and 21H03504.

Review summaries:
This paper studies the identification and estimation problem of the joint probabilities of potential outcomes. To solve this problem, this paper proposes two sets of identification conditions using covariate information. Furthermore, this paper proposes a novel statistical estimation method based on the augmented Lagrangian method to evaluate the joint probabilities of potential outcomes under the proposed identification conditions. 

This paper proposes two nonstandard identification strategies for the joint distribution of counterfactuals. One is based on a single proxy $Z$ for the unknown confounder with at least $4$ levels, such that $p(y|do(z))$ is identifiable. The other one is based two proxies $Z$ and $W$ with at least $4$ levels for each. Both identification strategies are fully nonparametric. They are given by solutions of a complex system which are nonlinear in the observed distributions. Beyond identification, the paper also proposes estimation procedures and proves the asymptotic normality of the estimators. The estimation procedures are essentially the plug-in version of the system involved in the identification strategy.