MissDAG: Causal Discovery in the Presence of Missing Data with Continuous Additive Noise Models
Authors: Erdun Gao
Ignavier Ng
Mingming Gong
Li Shen
Wei Huang
Tongliang Liu
Kun Zhang
Howard Bondell
Venue: NeurIPS 2022 Accept (2022)
Keywords: Causal discovery
Incomplete data
Additive noise models
EM method
Monte Carlo sampling

Abstract: State-of-the-art causal discovery methods usually assume that the observational data is complete. However, the missing data problem is pervasive in many practical scenarios such as clinical trials, economics, and biology. One straightforward way to address the missing data problem is first to impute the data using off-the-shelf imputation methods and then apply existing causal discovery methods. However, such a two-step method may suffer from suboptimality, as the imputation algorithm may introduce bias for modeling the underlying data distribution. In this paper, we develop a general method, which we call MissDAG, to perform causal discovery from data with incomplete observations. Focusing mainly on the assumptions of ignorable missingness and the identifiable additive noise models (ANMs), MissDAG maximizes the expected likelihood of the visible part of observations under the expectation-maximization (EM) framework. In the E-step, in cases where computing the posterior distributions of parameters in closed-form is not feasible, Monte Carlo EM is leveraged to approximate the likelihood. In the M-step, MissDAG leverages the density transformation to model the noise distributions with simpler and specific formulations by virtue of the ANMs and uses a likelihood-based causal discovery algorithm with directed acyclic graph constraint. We demonstrate the flexibility of MissDAG for incorporating various causal discovery algorithms and its efficacy through extensive simulations and real data experiments.

MissDAG: Causal Discovery in the Presence of Missing Data with Continuous Additive Noise Models
Erdun Gao
Ignavier Ng Mingming Gong'
Li Shen
Wei Huang'
Tongliang Liu Kun Zhang 2,5 Howard Bondell'
1 The University of Melbourne, Carnegie Mellon University,
3 JD Explore Academy, 'Tee University of Sydney
Mohamed bin Zayed University of Artificial Intelligence erdun.gaoos gaodstudent unimelb. edu. au, {ignavierng, kunz1) edu ggngg, gong, bei.huang, howard. bondell)@unnmeen edu, au chenlil00@jd. com, tongliang. liudsydney....... edu.au
{mingming.
Abstract
State-of-the-art causal discovery methods usually assume that the observational data is complete. However, the missing data problem is pervasive in many practical scenarios such as clinical trials, economics, and biology. One straightforward way address the missing data problem is first to impute the data using off-the-shelf imputation methods and then apply existing causal discovery methods. However, such two-step method may suffer from suboptimality, the imputation algorithm may introduce bias for modeling the underlying data distribution. In this paper, we develop general method, which call MissDAG, to perform causal discovery from data with incomplete observations. Focusing mainly on the assumptions of ig- norable missingness and the identifiable additive noise models (ANMs), MissDAG maximizes the expected likelihood of the visible part of observations under the (EM) framework. In the E-step, cases where comput- ing the posterior distributions of parameters in closed-form is not feasible. Monte Carlo EM is leveraged to approximate the likelihood. In the M-step. MissDAG leverages the density transformation to model the noise distributions with simpler and specific formulations by virtue of the ANMs and uses a likelihood-basee causal discovery algorithm with directed acyclic graph constraint. We demonstrate the flexibility of MissDAG for incorporating various causal discovery algorithms and its efficacy through extensive simulations and real data experiments.
appectationmmaximiaatie
Introduction
Discovering the underlying causal relations among variables of interest often occupies prominent position for supporting stable inference and rational decisions [42] in many applications such as medical diagnostics [48], recommendation systems [64] and economics [27]. To achieve this goal, conducting randomized controlled trials or using interventions is often acknowledged as the golden rule, which is effective but challenging in practice owing to high costs, ethical issues, difficulties in obtaining compliance [47]. To address this issue, causal discovery from purely observational data, realistic in specific settings, has drawn considerable attention in both academic
which may be and industrial fields [58, 24, 23, 21].
Existing causal discovery methods, such as constraint-based methods [57, 13], score-based methods [12, 46], and methods based on functional causal models [53, 25, 72, 44], typically focus on the
"qqual contribution.
i Work was done during an internship at JD Explore Academy.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

-xoll=
20
30
SHD
10-
20-
0.1
0.2
0.5
0.5
0.5
0.1
0.2
0.3
0.4
0.5
Missing rate (MCAR)
Missing rate (MCAR)
Complete data
MissDAG
Geussian-EN imputation
Figure 1: Example of linear Gaussian model with equal noise variance. With zero-mean noises, the recovered covariance matrix Σ (sufficient statistic) can be a criterion for distribution recovery.
settings which complete observations are available. However, practice, datasets often suffer from missing values caused by various factors such entry errors, deliberate non-responses, and sampling drops [31]. Following the definitions by Little and Rubin [31], the missing types can be categorized into three classes, namely missing completely at random (MCAR), missing random (MAR), and missing not at random (MNAR), according to different missing mechanisms. Many previous efforts have focused specifically on figuring more identifiable MNAR cases [6, 34, 37] and estimating the causal graphs from some specific MNAR cases [17, 62). while less attention has been paid to the M(C)AR cases extra assumption is required recover the ground-truth data distribution from incomplete observations [51, 34].
To perform causal discovery in the M(C)AR case, naive approach to handle the missing values is the listwise deletion method that simply drops the samples with missing value(s) in at least one of the variables. However, this may lead to unsatisfactory performance if the sample size is limited or the missing rate is high [59, 62] because of the decreased statistical power. Another straightforward approach is impute those missing values instead of listwise deleting them. However, these imputation methods may introduce bias for modeling the underlying data distribution [29]. Moreover, shown in Fig. (1), even though Gaussian-EM imputation method can consistently recover the data distribution as there is no model misspecification, it may lead to sub-optimal directed acyclic graph (DAG) estimation because it focuses solely on distribution recovery instead of structure learning. Therefore, a principled causal discovery approach that can handle the M(C)AR case is needed.
Contributions. In this work, develop practical and general EM-based framework, called MissDAG, to perform causal discovery in the presence of missing data, in which the underlying missing mechanism is independent from the observed information, which includes the M(C)AR case. Considering the data generating model, we focus on the identifiable additive noise models (ANMs) [44], including the linear non-Gaussian model [53], linear Gaussian model with equal noise variance [45], and nonlinear ANMs [25, 44]; as a byproduct, our framework also accommodates the typical non-identifiable case, namely the linear Gaussian model with non-equal noise variances that can only be identified to Markov equivalence class [57]. The resulting MissDAG framework flexibly accommodates different score based causal discovery algorithms [12, 74, 39, 71] developed for complete data and can be potentially extended to deal with more general cases (e.g., the log- determinant term is incorporated in the likelihood function). Moreover, we conduct extensive experiments on variety of settings, including synthetic and real data, against many baselines to verify the effectiveness of MissDAG.
2 Preliminaries
Additive noise models. We adopt the notion of structural causal model (SCM) [42, 43] characterize the causal relations among variables. Each SCM M = (2, F) consists of the exogenous variable Z (2,, Lz. .2.................... ]. the endogenous variable set X' {X1, X2..., Xa). and the function set JF = (f1, /2, faj. Here, each function fi computes the variable X, from its parents (or
causes) Pax, and an exogenous variable Z;, i.e., X; f( Zz). In this work, we focus on specific class of SCMs, called the ANMs [25, 44], given by
X; (Pax) Z: i == 1,2, d,
(1)
where Z,, interpreted as the additive noise variable, is assumed to be independent with variables in Pax, and mutually independent with variables in Z 24

---

Causal graph. Each SCM M induces a causal graph, which we assume in this work to be DAG', The DAG GM (V.E) consists of a vertex set := {1, 21..................... d), in which each node corresponds to the variable X;, and an edge set IE S V where (i,j) E IE if and only if X; E Pax,- Let X = (2... Xa) be a random vector that includes the variables in X. and P(X) (with density p(=)) be the joint distribution of random vector X. We assume that there are no latent common causes of the observed variables (i.e., causal sufficiency), which, together with the acyclicity assumption, indicates that P(X) and induced DAG GM satisfy the causal Markov condition [42].
3 Problem definition
Notations. In this paper. we focus on the finite-sample setting with missing data. Consider a dataset D = (X, Y) that consists of N samples, where X E R Nd and Y E I xd Each row (X,, Y.). independently sampled from P(X, represents the i-th observation. Y is the indicator matrix that records the missing positions in X. i.c., Y O if X,y is missing and Yij presence of missing data, the fully observed data X is unavailable. For simplicity, let group all observed positions in X and M group all the missing positions. For each observation X,, let group the indexes of the observed part and m group the indexes of the missing part of X, Notice that o and m are different for different observations. Then, Xo = [X Yú 1, E [M], E [d]] includes all observed positions and XM = [X, Y,, = 0, i E j E [d] includes all missing positions X. Similarly, we have Xio [X,, Y, 1,j E [d]] and X,m [X: Yi = 0, E [d]].
otherwise. In the
Following [31], we define the full likelihood of the i-th sample (Xi.: Y,) as
(iin(Xion Y2,,,)) = / p( p(Xio, Ximi0) pl p(Y: v) dXim:
=P =P(Y. vh) / p(X.o. Xim 0) dXims
where the parameters vb govern the missing mechanisms and O include all the model parameters. The ignorable likelihood of X is defined as
Lignorable (Xioi 0) / p(Xio: Ximi 0) dX
Assumption 1 (Ignorable missingness [31]). The inference about parameters based on the ignorable likelihood evaluated only by X;o is the same as inference for 0 based on the full likelihood.
Ignorable missingness is important as it common assumption that is required by EM-style algorithms and also our method. It can be interpreted as belief that the available data is sufficient to "correct" the missing data, by assuming that the missingness and model parameters are distinct.
Task definition. Consider a distribution P(X, Y), where the marginal distribution P(X) is induced from a SCM M satisfying the assumption of ANM as defined in Eq. (1), and dataset D (X, Y). In practice, the observed distribution is only P(X., o, Y') satisfying Assumption I and the observational part of the dataset is D (Xo. Y). Our task is to learn the DAG G. from the dataset D.
4 MissDAG
In this section, we introduce our proposed method, called MissDAG, which leverages the penalized EM framework to iteratively identify the causal graph 9M and model parameters from the incomplete data D. In the M-step, MissDAG takes the log-likelihood applies penalty function as the score function to guide the search of model parameters. Instead of directly modeling the complex likelihood of the sample, MissDAG equivalently models the simpler noise distributions by the density transformation [355] since ANMs always limit the noise distributions to some specific distribution families. Moreover, the prior information of causal structure is also considered as an inductive bias reduce the variance of parameter estimation. In the E-step, the
the observational part of the sample and
p Indeed,, there may difference between DAG and a causal graph the directed edges of latter given a causal meaning that allows to answer interventional queries [28].
We use [N] = (1, 2, ... N represent the of all integers from to N.
3

---

log- seellikeihoodd function is integrated over the posterior of missing entries to obtain the expectation the closed-fore of the posterior available. Otherwise, Monte Carlo (MC) simulations [65] are adopted to numerically compute the expectation. The details are shown in the following subsections.
4.1 The overall framework of MissDAG
Leveraging the development of score based causal discovery methods, we follow the style that taking the log- likelihood of observations (only the observational part) and penalty function the score function. Then, the general form of the optimization problem can be written as
arg max S(O) = -log 0) PEN(!)), #
subject to Og E DAGs,
where # = (0g, fy), including graph learning part 0g and causal mechanisms learning part f denotes the parameters of an SCM M.. However, in some methods [74, 75, 39], fg can be absorbed into and induced from Zz. (EN(-) ()) is the penalty function and À is the penalty coefficient. With the i..dd assumption of each observation (X,, Y:), the score function S(0) can be written as
((0) =
Σ log p(X.... X.mi 0) dX.m - A (EN
Unfortunately, the closed-form solution of S(0) cannot be obtained. Since the score function in Eq. (4.1) gives rise to a penalized maximum log-likelihood estimation problem, we can take the iterative penalized EM method [11], which relates the parameters estimation of the SCM from logp(Xo: 0) to the same parameters estimation from the complete-data log-likeliiee lig-liieldhod log p(X; 0).
Different from the imputation methods that replace the missing entries XM by some specific values, EM based methods formulate a two-step iterative operation. We start with an initial value 6 and denote 8* as the estimate of 6 at the t-th iteration. Then, each iteration of the EM method can be represented as the following two steps:
E-step takes the estimated model parameters of the previous step and the observational part Xo to impute the missing entries by the distribution of XM. which, in other words, gets the expected log-likelihood of the complete-data follows.
Q 0*) f P(XM Xo: 0t) log p(opp(Xoo. XM: dXM À PEN(()
(3)
Exx Xu:0r (logp(Xo. XM; 0)} PEN(().
M-step calculates o++1 by maximizing the Q function follows.
a++1 arg max 2(0,0% 0*)
The E-step calculates the expected data log-likelihood @(0, 0 The following M-step, then, maxi- mizes Q00, 0²) in for the fixed Or with DAG constraint. The convergence analysis of MissDAG is provided in Appendix E.
4.1.1 Log-likelihood term of S())
The final problem comes to the exact formulation of p(X; 0), which is not straightforward to obtain in our problem since ANMs typically impose assumptions the noise distribution instead of the joint distribution P(X). Benefiting from the well-researched results of the change of variables rule of density transformation [35], we can equivalently formulate p(X; 0) by transforming p(Z; 8). Here, for simplicity, we take Z (Z,, Zz.... Za) and f (f1, f.....,ffu) Then, we have
px (X) (x f(x)) det (I Jn)1:
(4)
where I is the identity matrix and y represents the Jacobian of f evaluated on X. Proposition 1. 0g represents DAG, then det(I J 1.

---

The proof is included in Appendix B.1. With Proposition 1, the log- log-determinant term becomes zero if the candidate solution is acyclic and can be dropped to simplify Eq. (5). However, we would also like to point out that some recent ccore-based structure learning methods based on continuous optimization, e.g., GOLEM [39] and NOTEARS-ICA [73], have shown that including the log-determinant term (which corresponds to likelihood based on directed cyclic graphs) in its objective function may be desirable and lead to better performance. In this work, this term is ignored since only acyclic models considered. Then, with Eq. (4) and Proposition I, we can model the log-likelihood of the simpler mutually independent noises distributions. Then, we have
d 2(X;;0) = Σ
d Σ
logp2, (X,y - fz(X,)) + log det(I-JJp 1=1
log (Xi, - (X)). (5) j=l
In the E-step, with different data generation models, the exact formulation of the posterior may be unavailable. In the next section, therefore, we split up these cases named exact posterior and approximate posterior respectively for presentations.
5 Different posterior cases
5.1 Exact posterior
Firstly, we deal with the linear Gaussian models including linear Gaussian model with equal variance (LGM-EV) and linear Gaussian model with non-equal variance (LGM-NV), That is to say, in Eq. (I), each fi E F is a linear function and Zj ~ (1 of) with Σg = diag(o) or........ g2). Then, the model can be rewritten as
X= W X
(6)
where W [add is the weight matrix and Wij 0 means that X; is one of the causes of Xj- Based on the density of a linear transformation, we know P(X) belongs to a multivariate Gaussian distribution. For multivariate Gaussian, sufficient statistics consist of the mean vector u(X) (the first-order moment) and the covariance matrix cov(X) (the second-orde moment). With zero-mean assumption of Z, we have p (X) and estimate cov by T. In other words, we can equivalently replace p pX;; 0) by P(T; 0). Specifically, with full data, T can be directly calculated by T = X
In this case, we specify 0 (W, Σz) since the parameters govern the distribution P(X). According to Eq. (5), the complete log-likelihood L(X; W, Σz) can be sufficiently formulated by
C(:: W, E2) log WzEz)
1
Tr(log Ez)
log 2T.
(7)
2N
Since (7) is linear in T, the Q function can be formulated in a closed-form.
5.1.1 E-step: compute Q(W, Σz, W .12
As discussed before, the E-step calculates the expected log- likelihood Jogp(T: W, W2) with Xo and (W*, Σz)). From Eq. (6), we have (I WT)-1 Z. Then, the implicit parameter Mx can be estimated by Σxx (I-W)- (I W) Then, T E|T Xo; Ex] can be straightforwardly calculated by the well-known results on the conditional distributions of the multivariate Gaussian.
Fach entry of T is obtained by Tú = - E
&A Xo: Ex}, where
Exny
XAiXA if X missing, X1:Xes
&
otherwise.
X records the expectation of missing part of the k-th instance X, and also shares the same indexes with X X₂ is initialized as Xx. Then, X Et X' mo _ X
5

---

5.1.2 M-step: maximize Q(W W", E2)))
With T the M-step maximizes the score S(0) with L(T"; W Σ2) - log pT;; W, Eg). Here, also plug in the DAG constraint for optimization to ensure that the estimated graph is acyclic. Then, the overall optimization problem can be formulated as
W E++
arg max W, + PEN(WW, W.Z
subject to Iw E DAGs,
where 9w means the graph induced from W. Here, we do not restrict the use of any specific algorithm for solving this problem. For the LGM-EV, one can adopt the greedy search method by Peters and Bühlmann [45], GOLEM [39], NOTEARS [74] to estimate the DAG G, while for the LGM-NV, one can apply different search methods like GES [12], A [71], and GOLEM [39].
5.2 Approximate posterior
Unfortunately, for Non- Linear (NL)-ANMs and Linear Non-GAussian Model (LiNGAM) cases, the E-step is not available in a closed-form, which would make the likelihood inference for missing data more difficult. There are mainly two problems to compute the expectation: (1) computing the posterior distribution p( im ))); (2) computing the integral. Since the problem mainly comes from the noise modeling as we later show, we prefer to first introduce the M-step of our method.
5.2.1 M-step: maximize Q(0,
From Eq. (5), we know that two reasons lead to the non-closed-fome for integral. The first one is that f includes some complex non-linear functions. However, considering the likelihood issue, the non-linyarity problem can be well handled by taking neural networks to model f [75, 40, 30]. Therefore we leave f for both linear and non-linear models for brevity. The second problem is non-Gaussian noise, which needs be clarified before the M-step.
The problem of (Z). If noises are Gaussaas,, the exact formulation of each noise distribution unknown. To make the likelihood-basse based methods work, a fixed Super (Suss-Gaussian prior distribution can be set to model the noise distributions [26, 75]. The theoretical result that the maximum likelihood estimate is locally consistent even in the presence of small misspecification error is well-established [2]. Here, we take the
distribution as an example:
(Z.) cz exp log cosh (Z.)).
(8)
where is a constant. Notice that using (8) needs to assume the unit scale of noise, which may not hold for the real scenarios. Therefore, we prefer to take normalized likelihood for standardized noise variable Z,/o, Then, the standardized ligeeiold of the observation X, will be
X
((X.) Jj
Estandardized (x;;0) log PZ,
log oj
where o3 = Var(Z,) is the variance of noise Zj with empirical version of x [X.j fx))) Then, Eq. (9) serve the log-likelihood term and the overall optimization problem in the M-step can be written as
g++1=
arg max Lstandardized (X::0) PEN(O),
subject to fg DAGs,
where Lstandardized (X represents the expected log-likelihood function in the f-th iteration.
5.2.2 E-step: compute Q(0,0%)
The problem of the integral. For LiNGAM and NL-ANMs, there no closed-form for the integral to get the expectation of lig-likelihood in the E-step. Naturally, we take the Monte Carlo sampling
'The term "standardized log-likelln likelihood" is taken to follow the literature of independent component analysis.

---

method [50] to approximate the expectation. Then, Eq. (3) can be reformulated as
N N, log p( X mmmi 0);
2 0*) = -xx flogp(Xo. 0)) =
where a sampled from the posterior, represents the j-th value of the total N, sampling results for the missing part Xim of the observation X,
The problem of sampling from the posterior. Instead of directly sampling from p(X X 8*) to fill the missing part Xim. we use rejection sampling [8] sample from proposal distribution (Xmm) with the probability density function ()mm. from which we can readily draw samples. Then, constant eı is set to guarantee that (Xg((m) 2 p(Xim Xio: 0*) for v. For each sample vi from However, p(Xim X can
q(Xm) the accept rate would Paccept = p(r O')/:'))ee00
not be directly obtained. With Bayes Rule, we have
p(X,;0') p (Xi 8*)
p X xio 0')
while for each instance X,, po p(X,; is still not directly available. With Eq. (4) and Proposition 1, we find that may skip for obtaining the dosed-form of p(X;; but equivalently provide the value of pz (Z;
:f:00 0¹) is a constant marked as Co. However, the joint distribution
p p(X:;0') = [I Jg; Jgg))/ (Z,; 0²) = II pz, XX:j - 0f (X:)).
Then, just like in the M-step, we take a normalized distribution to model the noise distributions. Then, the probability that sample can be accepted is given by
I pz, (Xis 0j (X:))
Paccept
CoCEq(Tm))
Here, we take = CoCk as a re-normalized constant to activate the rejection sampling methods.
6 Experiments
We report the empirical results to verify the effectiveness of MissDAG on both synthetic and biological dataset.
Baselines. We mainly take imputation methods as baselines including Mean Imputation, MissForest Imputation [60], and Optimal Transport (OT) immuutttutt [36] to impute the incomplete data first,' and then apply the causal discovery methods including GOLEM [39], NOTEARS [74], the algorithm ('Ghoshal') ) by Ghoshal and Honorio [20], NOTEARS-MLE [75], and NOTEARS-ICA [73] to estimate the causal graph for different assumed models (see Appendix D for details.). For L.M,, also include the structural EM method for multivariate Gaussian distribution, which we called Gaussian-F imputation, to recover the complete data. For LGM-NV that aims to identify the Complete Partial DAG (CPDAG), we also include differrnnt searching methods such as A* [71] and GES [12] to solve the optimization problem of the M-step. Also, we include the Test-wise Deletion PC (TD-PC) a baseline for LGM-NV, For LiNGAM, we also use ICA-LiNGAM and Direct- LiNGAM as baseline methods and put the results in Appendix F.7.2. The detailed implementations of the imputation and structure learning methods, as well as the hyper-parameters of the proposed method, are presented in Appendix C.
Metrics. We report the widely used criterion named Structural Hamming Distance (SHD), which refers to the smallest number of edge additions, deletions, and reversals required to transform the recovered DAG into the true one, averaged over 10 random repetitions to assess how the edges differ between the estimated and ground-true -ruth DAG in the identifiable cases. For the non-identifialle case such LGM-NV, report the SHD-CPDAG measure the distances between different CPDAGs. Other criteria such as Fl and recall are included for the supplementary experiments in Appendix F.
'Results on more imputation methods including GAIN [68]. KNNImputer [3], MICE [66] are shown in Appendix F.I.

---

COLIM-N
NOCTARS
Clold
0312W-0
SHENDIAO
82 63
04
03
34
Midta (H
Goegan
lepictioo
OTiquss
Cooplns
(a) LGM-EV.
(b) LGM-NV.
Figure 2: Recovery of the true structure measured by SHD or SHD-CPDAG (.). (a) LGM-EV with d = 20. (b) LGM-NV with d = 15 since the searching time of A is too long. Rows: ERI, ER2, and SF2 graphs. Columns: different methods. Some results for Gaussian- EM imputation are truncated because its SHDs are too large in those cases.
Simulations. The synthetic data we consider here is generated according to the ANM in Eq. (1). As illustrated above, consider four cases, including LGM-EV, LGM-NV, LiNGAM, and NL-ANM. In each experiment, a ground ground-truth truth DAG & with d nodes and kd directed edges was first generated from one of the two graph models, Erdös-Rén Rénds (ER) or Scale-Free (SF). According to different edges, kd edges (k = 1, 2). the graph model is named ERA or SFk. We also simulate denser graphs in Appendix F.6. Then, for linear models, a weighted matrix W E Rdxd with coefficients sampled from Uniform 2.0.0 -0.5 [0.5,.00] 2.0]) with equal probability generated to assign values each edge in G. For the non-linear model, corresponding to each edge in G. a ʃa constructed from a fixed MLP with random coefficients. The non-Gaussian noise we take here follows Gumbel distribution. In line with the settings outlined in [73], we do not consider the scenario of non-equal scales, as the normalized likelihood makes the optimization problem hard be solved. Our framework, however, can be extended with future advanced methods to tackle this challenge. Experimental results on more different non-linear functions and different noise distributions are included in Appendix F.8 and Appendix F.7.1, respectively. In equal variance/scale sets, all independent noises belongs to their distributions with variance/scale as while non-equal variance/scale settings get the scale of each noise independently sampled from Uniform| 1.0, 2.0]. For each experiment, we sample 100 non-linear models. More results on different
observations for linear models and 200 observations numbers of samples with fixed number of nodes and different numbers of nodes with the fixed number of samples are shown Appendix F.3 and Appendix F.4, respectively. We also add the experiment with d nodes and 2d observations with run-time comparisons in Appendix F.5. The missing type in our experiments is MCAR while the results on MAR and MNAR are also provided in Appendix F.2.
Linear Gaussian case. In Fig. (2a), across all settings for LGM-EV, including different graphs and missing rates, MissDAG with GOLEM, NOTEARS and Ghoshal as baseline methods can show consistently the best performance performance comparable to the best performances. While all imputation methods are sensitive different baseline methods, all NOTEARS- -bsse methods show improvements compared to GOLEM-basee methods although GOLEM is the real full likelihood method for LGM. However, GOLEM solves the problem by soft-constraitt, which may suffer from the finite sample. For imputation methods, we can see that MissForest usually acquires the second- best place. Gaussian EM imputation can consistently recover the multivariate Gaussian distribution, but its performance varies a lot with different causal discovery methods. With NOTEARS, Gaussian- EM imputation can gain better, even the best results as compared to the others. The results of LGM-NV are shown in Fig. (2b). Three different searching strategies are considered. It is observed that Mis can still achieve the best performances across all settings. The capacities of A and GES may be limited by the finite observations while GES appears to perform the worst. Moreover, A am exact search method, may suffer from high computing complexity. The comparisons on running time are provided in Appendix G.
Visualization of the learned DAG of MissDAG. We take an example of the MissDAG optimization process LGM-EV and plot the change in estimated parameters in Fig. (3), which shows that the

---

12
11
(a) Ground truth.
Optimization process; ((()) refers to estimated graph n-th EM iteration.
1N Web
10
(floratoc)
(Chaatioo
#(Deniis)
(c) LGM-EV.
Figure 3: Visualizations of the optimization process.
learned causal graph asymptotically approximates the ground-truth DAG g, including the existence of edges and their weights. The data distribution can also be well recovered.
310
40-
38
20
10-
14 Reatame MODR
MastaTer men MOCAD
MuDkO
HADIO
(a) LiNGAM.
(a) LGM.
20
Miskggnss (SCAR)
03
6.1
0.5
0.
0.3
013
0.5
05
Coopax
Muttan
(b) Non Linear ANMs.
(b) LINGAM.
Figure 4: MissDAG with approximate posterior.
Figure Dream4 results.
LiNGAM and NL-ANMs. Fig. (4a) shows the performances of MissDAG and different baseline methods with NOTEARS-ICA as the causal discovery algorithm on LiNGAM while Fig, (4b) provides the results on NL-ANM of 20 variables. One can see that MissDAG always occupies the best one of the best methods across all settings. MissForest [60] always shows the best or the second-best performance. Especially for NL-ANM, MissForest can get comparable results to MissDAG.
Biological dataset. We take biological dataset named Dream4 and provided [22], which simulates gene expression measurements from five sub-networks of transcriptional regulatory networks of E. coli and S. cerevisiae. Here we consider the 10-node networks, which however include feedback loops. We
that MissDAG occupies the first place for most of the settings and LGM is more suitable to learn from this dataset. The results in Fig. (5) show that there are misspecifications between our models and real data, probably due the cycles in the real data.
Related works
Causal discovery from complete data. Two lines methods prevail in causal discovery research, namely constraint-based methods, such PC and fast causal inference (FCT) [58], and score based methods like GES [12]. The first branch reads the (conditional) independencies information encoded in the data distribution, which can also be viewed as equality constraint [55], to decide the existence

---

of edges and directions of some edges. However, the constraint-baser method can only reach the Markov Equivalence Class (MEC) of the ground-truth DAG since DAGs in the same MEC share the totally same conditional independencies information. The second branch searches the model parameters in the DAG space by maximizing the penalized likelihood (score) on the observational data. For a long time, these methods suffer from the high searching complexities of combinatorial optimization. Recently, NOTEARS [74] recasts this problem as a continuous optimization by introducing an algebraic characterization of DAG. Then, NOTEARS has been extended to handle non- linear cases [69, 40, 76, 30, 63], time-series data [41], unmeasured confounder [7] and interventional data [9]. By imposing further assumption on the data generating model, NL-ANMs [25, 44], Post- Nonlinear Model (PNL) [72), LGM-EV [45]. and LiNGAM [53], etc, are proposed to learn the ground-truth DAG with identifiability guarantees.
Causality with incomplete data. MissDeeppCasal [33] leverages the deep latent model to estimate the causal effects of treatment, intervention or policy from incomplete data. GINA [32] system- atically analyzes the identifiability of generative models under MNAR case and designs a practical deep generative model which can provide identifiability guarantees for certain MNAR mechanisms. DECI [18] proposes a general deep latent model to perform both causal discovery and inference. Moreover, the theoretical results can guarantee that this model can identify the causal graph under standard causal discovery assumptions.
BN learning with incomplete data. Previous methods [23, 56, 15, 16, 49] mainly inherit the Expectation-Maxxiiizatoonn (EM) method framework [51], which conducts likelihood inference in an iterative optimization way. Friedman [15], Singh [56] iteratively refine the conditional distributions and sampling the missing values from these distributions. [61] proposed data augmentation method by stochastic simulation-based method that draws the filled-in value from predictive distribution. The augmentation method was accelerated by recasting the problem into two phases, parent set identification by an exact search and structure optimization by an approximate algorithm [1]. However, existing BN learning methods from incomplete data focus on identifying the Markov equivalence classes (i.e., discrete cases) under suitable assumptions and usually formulate the structure learning problem as a discrete optimization program, while our work focuses on continuous identifiable ANMs of which the structure is fully identifiable and includes recent structure learning approaches based on continuous optimization. More discussions can be found in Appendix A.
Conclusion and Future Work
In this paper, we propose new approach named MissDAG learn the underlying causal relations from incomplete data. MissDAG, leveraging the EM-based paradigm, iteratively maximizes the likelihood of the observational part of data with the inductive bias of DAG structure. Existing score- based causal discovery methods can be directly integrated into our framework for graph and model parameter learning. Moreover, MCEM is introduced to address the challenge that the closed-form posterior of missing entries is unavailable. The experiments show that our method works well across various of settings. However, our method inherits the time inefficiency issue EM algorithm. Future works include (1) improving the sampling efficiency with more efficient sampling or variational inference techniques to approximate the posterior to scale up to larger problems,
(2) incorporating other advanced causal discovery methods into the MissDAG framework, and (3) allowing unobserved confounders and cycles, c.g., using the methods by Bhattacharya et al. [7], Ghassami et al. [19].
Acknowledgements
LS is supported by the Major Science and Technology Innovation 2030 "Brain Science and Brain-like Research" key project (No. 2021ZD0201405), IN and KZ were partially supported by the National Institutes of Health (NIH) under Contract R01HL159805, by the NSF- Coveeennee Accelerator Track-D award #21 #2134901, by grant from Apple Inc., and by grant from KDDI Research Inc.. EG is supported by an Australian Government Research Training Program (RTP) Scholarship. This research was undertaken using the LIEF HPC-GPGPU Facility hosted at the University of Melbourne. This Facility was established with the assistance of LIEF Grant LEI 70100200. MG was supported by ARC DE210101624. TL was partially supported by Australian Research Council Projects DP180103424, 19010177,,,,, IC- 1901000,,,,,, 2P-220102121, and FT-22010031 18. HB was supported by ARC FT FT190100374.
10

---

References
[1] Tameem Adel and Cassio P De Campos. Learning bayesian networks with incomplete data by augmentation. In Association for the Advancement of Artificial Intelligence, 2017.
[2] Shun-ichi Amari, Tian-Ping Chen, and Andrzej Cichocki. Stability analysis of learning algo- rithms for blind source separation. Neural Networks, 10(8):1345- 1351, 1997.
[3] Rebecca R Andridge and Roderick JA Little. A review of hot deck imputation for survey non-response. International Statistical Review, 78(1):40- 64, 2010.
[4] Dimitri P Bertsekas. Constrained Optimization and Lagrange Multiplier Methods. Academic Press, 1982.
[5] Dimitri P Bertsekas. Nonlinear Programming. Athena Scientific, 2nd edition, 1999.
[6] Rohit Bhattacharya, Razieh Nabi, Ilya Shpitser, and James M Robins. Identification in missing data models represented by directed acyclic graphs. In Uncertainty in Artificial Intelligence, pages 1149-1158. PMLR, 2020.
[7] Rohit Bhattacharya, Tushar Nagarajan, Daniel Malinsky, and Ilya Shpitser. Differentiable causal discovery under unmeasured confounding. In International Conference on Artificial Intelligence and Statistics, 2021.
[8] Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
[9] Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin. Differentiable causal discovery from interventional data. In Advances in Neural Information Processing Systems, 2020.
[10] Tony Cai, Weidong Liu, and Xi Luo. A constrained ll minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106(494):594 607, 2011.
[11] Lin Chen, Ross L Prentice, and Pei Wang. A penalized em algorithm incorporating missing data mechanism for gaussian parameter estimation. Biometrics, 70(2):312-322, 2014.
[12] David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507-554, 2002.
[13] Diego Colombo, Marloes Maathuis, Markus Kalisch, and Thomas Richardson. Learning high- dimensional directed acyclic graphs with latent and selection variables. The Annals of Statistics, 10:294-321, 2011.
[14] Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical Lasso. Biostatistics, 9:432-41, 2008
[15] Nir Friedman. Learning belief networks in the presence of missing values and hidden variables. In International Conference on Machine Learning, 1997.
[16] Nir Friedman. The bayesian structural em algorithm. In Conference on Uncertainty in Artificial Intelligence, pages 129-138. 1998.
[17] Alexander Gain and Ilya Shpitser. Structure learning under missing data. In International Conference on Probabilistic Graphical Models, pages 121-132. PMLR, 2018.
[18] Tomas Geffnnr, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, et al. Deep end-to-end causal inference. arXiv preprint arXiv: 2222002295. 2022.
[19] AmirEmad Ghassami, Alan Yang, Negar Kiyavash, and Kun Zhang. Characterizing distribution equivalence and structure learning for cyclic and acyelic directed graphs. International Conference Machine Learning, 2020.

---

[20] Asish Ghoshal and Jean Honorio. Learning linear structural equation models in polynomial time and sample complexity. In International Conference on Artificial Intelligence and Statisties, pages 1466-1475. PMLR, 2018.
[21] Clark Glymour, Kun Zhang, and Peter Spirtes, Review of causal discovery methods based on graphical models. Frontiers Geneties, 10:524, 2019.
[22] Alex Greenfield, Aviv Madar, Harry Ostrer, and Richard Bonneau. Dream4: Combining genetic and dynamic information to identify biological networks and dynamical models. PloS one, 5 (10)ee3337,, 2010.
[23] David Heckerman. A tutorial on learning with bayesian networks. Innovations in Bayesian networks, pages 33-82, 2008.
[24] David Heckerman, Christopher Meek, and Gregory Cooper. A bayesian approach to causal discovery. In Innovations in Machine Leurning, pages 1-28. Springer, 2006.
[25] Patrik Hoyer, Dominik Janzing, Joris M Mooij. Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, 2008.
[26] Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer. Estimation of a structural vector autoregression model using non- non--aussiaaitt. Journal of Machine Learning Research, 11 (5), 2010.
[27] Guido W Imbens and Donald Rubin. Causal inference statisties, social, and biomedical sciences. Cambridge University Press, 2015.
[28] Daphne Koller and Nir Friedman. Probabilistic graphical models: principles and techniques. MIT press, 2009.
[29] Trent Kyono, Yao Zhang, Alexis Bellot, and Mihaela van der Schaar. Miracle: Causally- aware imputation via learning missing data mechanisms. In Advances in Neural Information Processing Systems, volume 34, 2021.
[30] Sébastien Lachapelle, Philippe Brouillard, Tristan Deleu, and Simon Lacoste-Julien Gradient- based neural dag learning. In International Conference on Learning Representations, 2020.
[31] Roderick JA Little and Donald B Rubin. Statistical analysis with missing data, volume 793. John Wiley & Sons, 2019.
[32] Chao Ma and Cheng Zhang. Identifiable generative models for missing not at random data imputation. Advances in Neural Information Processing Systems, 34:27645-27658, 2021.
[33] Imke Mayer, Julie Josse, Félix Raimundo, and Jean-Philippe Vert. Missdeepcausal: Causal infer- ence from incomplete data using deep latent variable models. arXiv preprint arXiv:2002. 2020.
[34] Karthika Mohan and Judea Pearl. Graphical models for processing missing data. Journal of the American Statistical Association, pages 1-16, 2021.
[35] Joris M Mooij, Dominik Janzing, Tom Heskes, and Bernhard Schölkopf. On causal discovery with cyclic additive noise model. In Advances in Neural Information Processing Systems, 201
[36] Boris Muzellec, Julie Josse, Claire Boyer, and Marco Cuturi. Missing data imputation using optimal transport. In International Conference on Machine Learning, pages 7130-7140. PMLR, 2020.
[37] Razieh Nabi, Rohit Bhattacharya, and Ilya Shpitser. Full law identification in graphical models of missing data: Completeness results. In International Conference on Machine Learning, pages 7153-7163. PMLR, 2020.
[38] Ronald C Neath. On convergence properties of the monte carlo em algorithm. Advances in Modern Statistical Theory and Applications: a Festschrift in Honor of Morris Eaton, pages 43-62, 2013.
12

---

[39] Ignavier Ng, AmirEmad Ghassami, and Kun Zhang. On the role of sparsity and dag constraints for learning linear dags. In Advances in Neural Information Processing Systems, volume 33, 2020.
[40] Ignavier Ng. Shengyu Zhu, Zhuangyan Fang. Haoyang Li, Zhitang Chen, and Jun Wang. Masked gradient-based causal structure learning. In SIAM International Conference on Data Mining, pages 424-432, 2022.
[41] Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series data. In International Conference on Artificial Intelligence and Statistics, 2020.
[42] Judea Pearl. Causality: Models, reasoning, and inference. Cambridge University Press, 2009.
[43]
Judea Pearl, Madelyn Glymour, and Nicholas P Jewell. Causal inference in statistics: A primer. John Wiley & Sons, 2016.
[44] J Peters, JM Mooij, D Janzing, and Schölkopf. Causal discovery with continuous additive noise models. Journal of Machine Learning Research, 15(1):2009-2053 2014.
[45] Jonas Peters and Peter Bühlmann. Identifiability of gaussian structural equation models with equal eπror variances. Biometrika, 101(1):219-228 2014.
[46] Garvesh Raskutti and Caroline Uhler. Learning directed acyclic graph models based on sparsest permutations. Stat, 7(1) 2018.
[47] David B Resnik. Randomized controlled trials in environmental health research: ethical issues. Journal of Environmental Health, 70(6):28. 2008
[48] Jonathan G Richens, Ciarán M Lee, and Saurabh Johri. Improving the accuracy of medical diagnosis with causal machine learning. Nature communications, 11(1):1-9, 2020.
[49] Carsten Riggelsen, Learning bayesian networks from incomplete data: An efficient method for generating approximate predictive distributions. In SIAM International Conference on Data Mining, pages 130-140. SIAM, 2006.
[50] Christian Robert, George Casella, and George Casella. Monte Carlo statistical methods, volume 2. Springer, 2004.
[51] Donald B Rubin. Inference and missing data. Biometrika, 63(3):581-592, 1976.
[52] Gideon Schwarz, Estimating the dimension of a model. The Annals of Statistics, pages 461-464, 1978.
[53] Shohei Shimizu, Patrik O Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
[54] Shohei Shimizu, Takanori Inazumi, Yasuhiro Sogawa, Aapo Hyvärinen, Yoshinobu Kawahara, Takashi Washio, Patrik O Hoyer, and Kenneth Bollen. DirectLiNGAM: A direct method for learning linear -an-aassiann structural equation model. Journal of Machine Learning Research, 12(Ap)): 2011.
[55] Ilya Shpitser, Robin Evans, Thomas Richardson, and James M Robins. Introduction to nested markov models. Behaviormetrika, 41(1):3-39, 2014.
[56] Moninder Singh. Learning bayesian networks from incomplete data. In Association for the Advancement of Artificial Intelligence, pages 534-539, 1997.
[57] Peter Spirtes and Clark Glymour. An algorithm for fast recovery of sparse causal graphs. Social Science Computer Review, 9:62-72, 1991.
[58] Peter Spirtes, Clark Glymour, Richard Scheines, al. Causation, Prediction, and Search, volume 1. The MIT Press, 2001.
13

---

[59] Nicolas Städler and Peter Bühlmann. Missing values: sparse inverse covariance estimation and an extension sparse regression. Statistics and Computing. 22(1):219-235, 2012.
[60] Daniel Stekhoven and Peter Bühlmann. Missforest
Missforest-non-paraamee -on--noonnnrtt missing value imputation for mixed-type data. Bioinformatics, 28(1): 2012.
[61] Martin A Tanner and Wing Hung Wong. The calculation of posterior distributions by data augmentation. Journal of the American statistical Association, 82(398):528-540, 1987.
[62] Ruibo Tu, Cheng Zhang, Paul Ackermann, Karthika Mohan, Hedvig Kjellström. and Kun Zhang. Causal discovery in the presence of missing data. In International Conference on Artificial Intelligence and Statistics, pages 1762-1770. PMLR, 2019.
[63] Xiaoqiang Wang, Yali Du, Shengyu Zhu, Liangjun Ke, Zhitang Chen, Jianye Hao, and Jun Wang. Ordering-basse causal discovery with reinforcement learning. In International Joint Conference on Artificial Intelligence, 2021.
[64] Yixin Wang, Dawen Liang, Laurent Charlin, and David M Blei. Causal inference for recom- mender systems. In ACM Conference on Recommender Systemx, pages 426-431, 2020.
[65] Greg CG Wei and Martin A Tanner. A monte carlo implementation of the algorithm and the poor man's 's data augmentation algorithms. Journal of the American statistical Association, 85 (411):99-7704, 1990.
[66] lan R White, Patrick Royston, and Angela M Wood. Multiple imputation using chained equations: issues and guidance for practice. Statistics in medicine, 30(4):377-399, 2011.
[67] CF Jeff Wu. On the convergence properties of the em algorithm. The Annals of Statistics, pages 95-103, 1983.
[68] Jinsung Yoon, James Jordon, and Mihaela Schaar. Gain: generative adversarial nets. In International Conference on Machine Learning, pages 5689- 5698. PMLR, 2018
Missing data imputation using
[69] Yue Yu, Jie Chen, Tian Gao, and Mo Yu. DAG-GNN: DAG structure learning with graph neural networks. In International Conference Machine Learning, 2019.
[70] Changhe Yuan and Brandon Malone. Learning optimal Bayesian networks: A shortest path perspective. Journal of Artificial Intelligence Research, 48(1):23 65, 2013.
[71] Changhe Yuan, Brandon Malone, and Xiaojian Wu. Learning optimal bayesian networks using a* search. In International Joint Conference on Artificial Intelligence, 2011.
[72] Kun Zhang and Aapo Hyvarinen. On the identifiability of the post-nonlinear causal model. In Conference on Uncertainty in Artificial Intelligence, 2009.
[73] Xun Zheng. Learning DAGs with Continuous Optimization. PhD thesis, Carnegie Mellon University, 2020.
[74] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Erie P Xing. DAGs with NO TEARS: Continuous Optimization for Structure Learning. In Advances in Neural Information Processing Systems, volume 31, 2018.
[75] Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric P. Xing. Learning sparse nonparametric DAGs. In International Conference Artificial Intelligence and Statistics, 2020.
[76] Shengyu Zhu, Ignavier Ng, and Zhitang Chen. Causal discovery with reinforcement learning. In International Conference on Learning Representations, 2020.
14

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full set of assumptions all theoretical results? [Yes]
(b) Did you include complete proof' of all theoretical results? [Yes]
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any new assets either in the supplemental material or URL? [No]
(d) Did you discuss whether and how consent obtained from people whose data you're using/curating? [No]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you deseribe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
15


Reviews:
[{"review_id": "m6zsO0RiN5Q", "reviewer": "Reviewer_woaA", "summary": "The paper proposes a method called MissDAG for learning causal structures in the presence of missing data. It assumes the data is missing at random (MAR) and the underlying causal model is an additive noise model (ANM). The paper derives expectation-maximization (EM)-based algorithms for maximizing a penalized log-likelihood score under different classes of ANMs including the linear non-Gaussian models, linear Gaussian models with equal noise variances, and nonlinear ANMs. The proposed MissDAG method is evaluated through synthetic and real data experiments.\n", "questions": "Why are MissDAG methods not empirically compared with existing work on causal discovery or Bayesian network (BN) structure learning from missing data?\n\nWhat are the sizes of the models (the number of variables $d$) in LiNGAM and NL-ANM experiments shown in Fig. (3)?\n\nWhat is the running time of MissDAG vs. baselines in the experiments?\n\nHow well does MissDAG scale to larger models?\n", "limitations": "The paper did not discuss the limitations of the proposed approach. One limitation I see is that the proposed approach may be too slow to learn large models which I believe is a well-known issue in the EM-type of algorithms for learning BN structures from missing data. The largest models tested in the experiments contain 20 variables. How well does MissDAG scale to larger models?\n", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper addresses an important problem of causal discovery in the presence of missing data. The existing work mostly focuses on structure discovery from discrete data while this work studies learning causal structures from incomplete data under ANMs. The proposed technical approaches are mostly standard techniques – its main contributions are deriving EM steps under different classes of ANMs. \n\nThe experiments only compared with imputation methods as baselines, but not with existing work on causal discovery or Bayesian network (BN) structure learning from missing data, e.g. Structural EM [16], MVPC [60].  Why are MissDAG methods not empirically compared with existing works?\n\nThe writing is overall clear, but there are grammar issues here or there.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "_wo92Vd4ot-", "reviewer": "Reviewer_vkvi", "summary": "This paper studies the problem of causal discovery when parts of the observation matrix are missing. The proposed method called “MissDAG” is an EM-like procedure: the E-step estimates the expected data log-likelihood as a function of the unknown causal graph and parameters by marginalizing out/sampling the missing values given the graph and parameters estimated in the previous step; the M-step then maximizes this expected data log-likelihood function with respect to the unknown causal graph and parameters. In an empirical evaluation, the work shows that MissDAG results in better graph recovery in terms of SHD/F1 scores than using other imputation methods.", "questions": "Method:\n\nA main question arising in Section 4 concerns the likelihood term $p(\\mathbf{X} | \\theta)$ (Section 4.1.1). Throughout the paper, we assume that the SCMs considered are *acyclic* (line 74 and 79). Thus, irrespective of the functional form and the noise distribution, the likelihood factorizes into a product of local conditionals $p_j(x_j | pa_j)$ (the Markov kernels). As a result, the likelihood is always trivial to evaluate conditional on the causal graph and parameters, especially for the additive noise models considered in this work.\n\nHence, I do not understand why this work introduces the complicated generalization in Eqs. (7)-(8) involving determinants and the Jacobian. This general form was used by [32] to handle *cyclic* SCMs, in which the likelihood is not trivially computable in the above, recursive form. However, the cyclic setting is not studied in this work. \n\nLater in the manuscript (Proposition 1, line 216), the authors even note that for a DAG, this additional term is 1 and can thus be dropped from Eqs. (7)-(8), which reduces the likelihood computation to the simple factorized computation always used in BN structure learning (This is even explicitly spelled out in lines 217-218). Unless I missed something essential about MissDAG and the incomplete data setting, this unnecessary complication in the likelihood somewhat indicates of a lack of familiarity with related work and the task of causal structure learning more generally. The whole of Section 4 could be significantly simplified and thus made more accessible when this is dropped.\n\nRelated to the above, I do not see why Eq. (10) needs to be so complicated. The likelihood in the linear Gaussian is trivial to evaluate and simply a product of Gaussian distributions with densities that are tractable to compute, without needing to introduce $\\mathbf{T}$. It would also be straightforward to add a bias term to the linear Gaussian model instead of using the structural equation model notation in Eq (9), where each causal variable has zero mean.\n\nRelated work:\n\nPlease clarify the precise relationship to and contribution over the prior works [15] and [54]. It appears that their algorithms conceptually do the same EM procedure as MissDAG, except that MissDAG substitutes in more modern likelihood model assumptions and structure learning methods into the subcomponents of the algorithm. For instance, [54] states: “[Our algorithm] uses the current estimate of the structure and the incomplete data to refine the conditional probabilities, then imputes new values for missing data points by sampling from the new estimate of the conditional probabilities, and then refines the network structure from the new estimate of the data using standard algorithms for learning Bayesian networks from complete data.”\n", "limitations": "Overall, the paper presents an intuitive idea based on EM for the causal discovery problem from incomplete data, which is validated convincingly in experiments on simulated data but is not new in the literature [15, 54]. In addition, given the points brought up above, the presentation and derivation have room for improvement both in terms of writing and in terms of clarity and motivation of the mathematical derivations, which sometimes indicates a lack of familiarity with existing approaches to the problem (in particular, the point on computing the likelihood). This makes an otherwise clear and intuitive method, which appears to be effective empirically, very hard to parse and internalize. Since the idea and results are convincing overall, I am open to revising my scores if the issues raised have been addressed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The motivation of the paper and the introduction are convincing. It makes sense to use EM to account for missing values in this context, and it has been done by previous work (e.g. [15, 54]). The array of imputation baselines and considered settings are likewise comprehensive and show the efficacy of MissDAG in the context of causal discovery. Since we are interested in the _causal_ graph, it would make sense to also show metrics like the structural intervention distance (SID), which is closer to our interests in causal inference.\n\nWhile the overall presentation structure of the paper is fine, the presentation of the various subparts and subsettings in the methods section (Section 4) is quite convoluted, making it hard to parse the key challenges and their proposed solutions. I think the exposition could be improved by first introducing the E-step and M-step and their different levels of difficulty (exact vs. approximate posteriors) without reference to a specific generative model yet (such as LGM, LiNGAM, NL-ANM). Then, to separate the MissDAG algorithm from the specific generative model assumptions, specific instantiations of these steps could be presented in an additional Section 5 that gives certain closed forms, as for instance in the linear Gaussian case.\n\nFurther points/questions about the presentation:\n- Figure 3: “$\\Sigma$ (sufficient statistic) can be a criterion for distribution recovery.” What is your justification for that? In this paper, the goal is causal structure learning and not (observational) distribution recovery, so the work should give some motivation to show this figure and metric. \n- Line 46: “the consistent approach”. What is the consistent approach that is meant here? This is related to the previous point about Figure 1, which is referenced in this sentence.\n- Figure 6(b): What is a “history”? This is never explained anywhere.\n- In line 56, you mention the concept of “ignorability” when stating your main contribution, even though this concept was never introduced at this point of the manuscript. It would be good if a sentence earlier in the introduction explains this notion and why it is needed.\n- I do not understand the point of introducing Eq. (2) and/or $\\phi$, since it is never used in the later sections of the paper. Can you provide some motivation? Also, why is “ignorable missingness” (Assumption 1, line 95) important? This is not motivated explicitly and never picked up again later in the manuscript, and unless the reader is very familiar with this literature already, it may be difficult to relate to the other parts of the paper.\n- Line 124: What is “$\\theta_\\mathcal{M}$”? This notation is not used anywhere else.\n- Line 128: What do you mean by imputing “the missing entries by the functions of $X_M$”? Do you mean the distributions of $X_M$?\n- Line 168-169: “we also plug in the DAG constraint as an inductive bias for optimization.” The DAG constraint is hardly an inductive bias, more a hard constraint (line 74), which usually makes the structure learning problem harder to solve and not easier. Please consider rephrasing.\n- Line 189-190: “… the exact formulation of each noise distribution is agnostic.” What is meant by this? This is most likely a wrong usage of the word “agnostic”.\n- Line 205-206: This is imprecise. You should indicate that $x^j_{im}$ are samples from the posterior $p(\\mathbf{X}_{m} | X_o, \\theta^t)$.\n- Lines 299-300: “Different from the BN learning methods, our method is specifically designed for causal discovery, where identifiability of the underlying causal graph is critical.”\nI do not agree with this statement. There is no inherent difference between causal discovery algorithms from observational data using a score- or constraint-based approach and algorithms for learning a BN. For certain classes of BNs (such as the ANM ones studied in this work), the true causal structure achieves highest likelihood, which is why these models are deemed “identifiable” from observational data. In fact, NOTEARS itself is a method for BN structure learning and used pervasively in this work as a causal discovery algorithm.\n\nTypos and minor suggestions:\n- Abstract (lines 6-7): “… a two-step method may suffer from suboptimality, as the imputation algorithm is unaware of the causal discovery step.”\nThis sentence is vague, especially since MissDAG is not “aware” of the causal discovery step either: MissDAG uses an estimated causal model to approximate the expected data log-likelihood. What does it mean to be “aware” of the causal discovery step? Formulations about efficiency and/or bias as in lines 42-45 would be more convincing.\n- Line 21: “concerned variables” -> consider rephrasing as, e.g. “variables of interest”\n- Line 83: typo, the second $\\mathbf{X}$ should be a $\\mathbf{Y}$\n- Line 87: typo, I believe “observed portions” should be “observed positions”\n- Line 118: “a graph learning part $\\theta_\\mathcal{G}$”. It appears as though this is the causal graph. If so, it may be clearer to use $\\mathcal{G}$ as the notation (similar to line 75 or 171), since the causal graph is usually considered to be separate from the parameters or functions parameterizing an SCM.\n- Line 119: “While in …, … .“ This sentence is grammatically incorrect.\n- Line 125: “which relates the estimation … to the same parameters estimation”. This sentence is grammatically incorrect. \n- In my opinion, the related work section would work better as Section 3. Otherwise, it is not transparent from the beginning what the gap in the state of the art is that MissDAG sets out to fill, especially since EM-like approaches have previously been used for BN structure learning.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "N8udch3wbd", "reviewer": "Reviewer_GKqy", "summary": "An EM-based method is described for doing causal discovery when some data entries are missing (MAR or MCAR). Different versions are proposed, which can deal with different variants of additive noise models. Experimental results show it tends to improve over imputation methods.", "questions": "* What can you say about the running time of your method compared to other methods for the same task?\n\n* The conclusion of appendix E (Convergence analysis) was unclear to me. Under what conditions is convergence guaranteed? Does this suffice to guarantee overall correct behaviour of the method (e.g. in the sense of finding the same graph as would be found without missingness, given enough data)?\n\nMinor points:\n* line 77: \"include\" -> \"includes\"; I recommend adding a comma before \"and\"\n* line 78: \"is no\" -> \"are no\"\n* line 80: \"indicate\" -> \"indicates\"\n* line 83: 2nd formula about X should be about Y, and in {0,1}$^{N \\times d}$\n* line 106: \"applys\" -> \"applies\"\n* line 115: \"follow the fashion that\": pick another word in place of \"fashion\" (then make the sentence grammatical)\n* line 124: $\\theta_{\\mathcal{M}}$: is this a new notation?\n* line 178: NL-ANM and LiNGAM are used without having been introduced, please mention the abbreviation when the work is cited earlier.\n* line 231: \"Completely Partial DAG\" -> \"Complete Partial DAG\"\n* line 278: I guess the \"is\" should be removed", "limitations": "-", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The problem is practically relevant and the theoretical analysis looks sound, though I did not check it closely. I think the clarity of the paper suffers somewhat from the large amount of material compared to the page limit.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "2_rLwW3R9Zg", "reviewer": "Reviewer_sgce", "summary": "The method proposed copes with the missing data problem in causal discovery by maximising the expected likelihood of the visible part of observations using the expectation-maximisation framework. ", "questions": "Questions. In Section 3, in notations,  X \\in R^{N \\times d} is repeated twice; the second should be Y? \n\nI found it not easy to follow Section 4.1.1. It is stated that \"benefiting from the well-researched results of density transformation...\": it might be well-researched but not extremely well known, and it is not easy to understand where equation 7 comes from. \n\nI am not sure that the notation in equation 8, where X comes twice (on the left and on the right) is successful. \n\nThe explanation what \"SHD\" stands for comes a little bit late in the text.  \n", "limitations": "I guess that the method is not scalable (since it relies on the EM).\nFrom the causal viewpoint: it is assumed that there is not any latent confounders and there are no cycles in the graph. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths. The paper is well written. The technical details are provided, and the results of numerical experiments are convincing. \n\nWeaknesses. The EM algorithm to treat missing data is quite a standard approach. Its extension to the current problem is rather straightforward. Difficult cases, such as latent common causes or cyclic graph are not considered.\n\nA DAG is not necessarily a causal graph, as far as I know. Section 4 is focused on an extension of EM to DAGs. It seems that there is nothing specific for causal inference in this contribution. \n\nIn Related Work: the papers cited are quite old. There is much recent work on missing data processing (but not necessarily for DAGs). \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept