Gene Regulatory Network Inference in the Presence of Dropouts: a Causal View
Authors: Haoyue Dai
Ignavier Ng
Gongxu Luo
Peter Spirtes
Petar Stojanov
Kun Zhang
Venue: ICLR 2024 oral (2024)
Keywords: Gene regulatory network
Single-cell RNA-sequencing
Dropout
Zero-inflated data
Causal model
Causal discovery
Nonparametric
Primary area: causal reasoning
TL;DR: The conditional independence (CI) relations in the data with dropouts, after deleting the samples with zero values for conditioned variables, are identical to the CI relations in the original data without dropout.

Abstract: Gene regulatory network inference (GRNI) is a challenging problem, particularly owing to the presence of zeros in single-cell RNA sequencing data: some are biological zeros representing no gene expression, while some others are technical zeros arising from the sequencing procedure (aka dropouts), which may bias GRNI by distorting the joint distribution of the measured gene expressions. Existing approaches typically handle dropout error via imputation, which may introduce spurious relations as the true joint distribution is generally unidentifiable. To tackle this issue, we introduce a causal graphical model to characterize the dropout mechanism, namely, Causal Dropout Model. We provide a simple yet effective theoretical result: interestingly, the conditional independence (CI) relations in the data with dropouts, after deleting the samples with zero values (regardless if technical or not) for the conditioned variables, are asymptotically identical to the CI relations in the original data without dropouts. This particular test-wise deletion procedure, in which we perform CI tests on the samples without zeros for the conditioned variables, can be seamlessly integrated with existing structure learning approaches including constraint-based and greedy score-based methods, thus giving rise to a principled framework for GRNI in the presence of dropouts. We further show that the causal dropout model can be validated from data, and many existing statistical models to handle dropouts fit into our model as specific parametric instances. Empirical evaluation on synthetic, curated, and real-world experimental transcriptomic data comprehensively demonstrate the efficacy of our method.

Published as a conference paper at ICLR 2024
GENE REGULATORY NETWORK INFERENCE IN THE PRESENCE OF DROPOUTS: A CAUSAL VIEW
Haoyue Dai' Ignavier Ng' Gongxu Luo2 Peter Spirtes'
Petar Stojanov
Kun Zhang
Department of Philosophy, Carnegie Mellon University
2 Machine Learning Department, Mohamed bin Zayed University of Artificial Intelligence Cancer Program, Eric and Wendy Schmidt Center, Broad Institute of MIT and Harvard
ABSTRACT
Gene regulatory network inference (GRNI) is challenging problem, particularly owing to the presence of zeros in single-cell RNA sequencing data: some are biological zeros representing no gene expression, while some others are technical zeros arising from the sequencing procedure (aka dropouts), which may bias GRNI by distorting the joint distribution of the measured gene expressions. Existing approaches typically handle dropout error via imputation, which may introduce spurious relations as the true joint distribution is generally unidentifiable. To tackle this issue, we introduce causal graphical model to characterize the dropout mechanism, namely, Causal Dropout Model. We provide a simple yet effective theoretical result: interestingly. the conditional independence (CI) relations in the data with dropouts, after deleting the samples with zero values (regardless if technical or not) for the conditioned variables, are asymptotically identical to the CI relations in the original data without dropouts. This particular test-wise deletion procedure, in which we perform CI tests on the samples without zeros for the conditioned variables, can be seamlessly integrated with existing structure learning approaches including constraint-baset and greedy score-based methods, thus giving rise to principled framework for GRNI in the presence of dropouts. We further show that the causal dropout model can be validated from data, and many existing statistical models to handle dropouts fit into our model as specific parametric instances. Empirical evaluation on synthetic, curated, and real-world experimental transcriptomic data comprehensively demonstrate the efficacy of our method.
INTRODUCTION
Gene regulatory networks (GRNs) represent the causal relationships governing gene activities in cells (Levine & Davidson 2005), essential for understanding biological processes and diseases like cancer (Ito et al., 2021; Huang et al. 2020 Parrish et al.) 202T). Typically, for p genes, the representing gene expressions, and directed edges representing cross-gene regulations. Traditional lab-based GRNI involves gene knockout experiments, but conducting all combinatorial interventions is challenging. In contrast, observational expression data is abundant by ANN-sequunc. In the past decade, single-cell RNA-sequencing (scRNA-seg) has become prevalent, enabling comprehensive studies in cancer Neftel et al. 2009);; Boiarsky et al. cancer[eetel
GRN is graph consisting of nodes Z V.
(202)) and genomic atlases Regev et al.) (2017); Consortiim et al. (2018) 2022) at the individual cell Causal discovery techniques for GRNI have also been steadily developed leverage these advances (Wang et al. 2017 Zhang et al.
2017 Belyaeva et al. 2021 Zhang et al. 2021 2023).
Despite the advantage of scRNA-seq, fundamental challenge known as the dropout issue arises. scRNA-seq data is known to exhibit an abundance of zeros. While some of these zeros correspond to
genuine biological absence of gene expression (Zappia et al. 2017 Alberts al. 2002), some others are technical zeros arising from the sequencing procedure, commonly referred to as dropouts (Jiang ctt al. 2022 Ding et al. 2020). Various factors are commonly acknowledged to contribute to the occurrence of dropouts, including low RNA capture efficiency (Silverman et al. 2020 Jiang et al.) 2022, Kim et al. 2020; Hicks et al. 2018). intermittent degradation of mRNA molecules (Pierson & Yau 2015 Kharchenko et al. 2014), and PCR amplification biases (Fu et al., 2018 Tung et al. 2017). The dropout issue has been shown to introduce biases and pose threat to various downstream
tasks, including gene regulatory network inference (Jiang al. 2022 Van Dijk al. 2018).

---

Published as a conference paper at ICLR 2024
Dealing with dropouts in scRNA-seq data has been approached through two main strategies. One
approach involves using probabilistic models such as zero- inflated models (Pierson & Yau [2015]
Kharchenko et al. 2014 Saced et al. 2020 Yu et al. [023] Min & Agresti 2005 Li & [i| Tang
et [ or hurdle models (Finak et al., 2015 Qiao et al. (023) to characterize the distribution
of gene expressions with dropouts. However, these methods may have limited flexibility due the
restrictive parametric assumptions (Svensson et al., 2018: Kim et al. (0220) Another approach is
imputation (Van Dijk et al. 2018 Li & Li 20188 Huang al. 2018 Amodio et al. 2019 Lopez et al. 2018), where all zeros are treated as missing values and imputed to estimate the underlying distribution of genes without dropouts. However, imputation methods often lack a theoretical guarantee due to the inherent unidentifiability of the underlying distribution. Empirical studies also demonstrate mixed or no improvement when using imputation on various downstream tasks (Hou et al., 2020) Jiang et al. 2022). Overall, despite various attempts, there is currently still not principled and systematic approach to effectively address the dropout issue in scRNA-seq data.
While GRNs represent the causal regulations among genes in the expression process, can also ex- tend this understanding to the causal mechanisms of dropouts in the sequencing process? With this mo- tivation, abstract the common understanding of dropout mechanisms by proposing causal graph- ical model, called causal dropout model ($2.2) We then recognize that the observed zeros in scRNA- seq data resulting from dropouts are non-ignorable, implying that the distribution of the original data is irecoverable without further assumptions ($2.4). However, surprisingly, we show theoretically that, given such qualitative understanding of dropout mechanism, we could simply ignore the data points in which the conditioned variables have zero values, leading to consistent estimation of conditional independence (CI) relations with those in the original data (2.5). This insight then readily bridges the gap between dropout-tainted tainted measurements and GRNI with an asymptotic correctness guarantee ($ We further provide a systematic way to verify such dropout mechanisms from observations ($7).
Contributions. 1) Theoretically. the proposed causal dropout model is, to the best of our knowledge, the first theoretical treatment of the dropout issue in a fully non-parametric setting. Most existing parametric models fit into our framework as specific instances. 2) Empirically, extensive experimental results on synthetic, curated, and real-world experimental transcriptomic data with both traditional causal discovery algorithms and GRNI-specific algorithms demonstrate the efficacy of our method.
2 CAUSAL MODEL FOR DROPOUTS AND MOTIVATION OF OUR APPROACH
In this section, we develop a principled framework to model the dropout mechanisms in scRNA-seq data, and describe the motivation of our proposed GNRI approach that will be formally introduced in $3 We first introduce causal graphical model to characterize the dropout mechanism in 2. and demonstrate in $2.3 how existing statistical models commonly used to handle dropouts fit into our framework as specific parametric instances. We then discuss in $2.4 the potential limitations of imputation methods, demonstrating that the underlying distribution without dropouts is generally unidentifiable. Lastly, we demonstrate in 22.5 that, despite the unidentifiability of the distribution, the entailed conditional (in)dependencies can be accurately estimated, which serves as the key motivation of our proposed GNRI approach in (grounded on the causal graphical model presented in 22.)
2.1 WARMUP: CAUSAL GRAPHICAL MODELS AND CAUSAL DISCOVERY
Causal graphical models use directed acyclic graphs (DAGs) to represent causal relationships (di- rected edges) between random variables (vertices) (Pearl 2009 Spirtes al. 2000). Causal discovery aims to identify the causal graph from observational data. In this paper, by "observational" we assume the data samples are drawn i.i.d., without interventions or heterogeneity. In the context of GRNI, we focus on cells produced in same scRNA-seq run without different perturbations (Dixit et al. 2016).
typical kind of causal discovery methods is called constraint- methods, e.g., the PC algorithm (Spirtes et al.; (000). These methods associate the conditional (in)dependencies (X T YJS) in data with the (un)connectedness in graph, characterized by the d-separation patterns (X LS)(S), and accordingly recover the underlying causal graph structure, up to a set of equivalent graphs known as Markov equivalence classes (MECs), represented by complete partial DAGs (CPDAGs). Note that CI tests do not require any specific assumptions about the underlying distributions of the variables (Rosenbaum) 2200 Zhang et al. (01T), constraint- based methods can in theory work nonparametrically for arbitrary distributions consistent with the graph.
2

---

Published as a conference paper at ICLR 2024
2.2 CAUSAL DROPOUT MODEL: A CAUSAL GRAPHICAL MODEL FOR DROPOUTS
Let G be a DAG with vertices partitioned into four sets, ZUXUDUR, Denote by Z = { the true underlying expression variables of the p genes. The subgraph of G on Z is the gene regulatory network (GRN) that we aim to recover from data. However, we do not have full access to true expressions Z. but only the sequenced observations X {X,) contaminated by dropout error. We use boolean variables D = P to indicate for each gene whether technical dropout happens in sequenced individual cells: D; = 1 indicates dropout and otherwise. Hence X are generated by:
X; = (1 - D,) Zz.
Similar to Z. we also do not have full access to D in practice (shown as the gray color in Figure I]: whether a zero is biological or technical is still unknown to us. To help address such issue, we further introduce zero observational indicators R [R: defined as R; D; OR 1 (Z; = 0), i.e., zeros are either technical or biological. By definition Equation @. R2 can be fully read from X, by R₂ 1 (X, 0). For any subset S [p] [1. ...... pp. we denote the corresponding random vector by e.g., Xs = {X, : E S}. Specifically, for boolean indicator variables, we use e.g., Rg- 1 to denote eR, 1, and Rs = 0 for ANDres R. = 0. We call such a proposed causal graphical model for modelling dropouts Causal Dropout Model (CDM) throughout this paper.
7A
D2
R
X,
Figure 1: Causal graph for dropouts. Gray nodes are un- derlying partially observed variables and white nodes are observed ones.
Note that in the causal graph depicted in Figure []] the node R. is represented by square, while the others are represented by circles. This distinction is made because R2 serves as an auxiliary indicator instead of a variable with an atomic physical interpretation. The directed edges into X, and Rg from Z, and D, are deterministic, i.e., the causal model can be sufficiently represented by Zr and any one of D;, R,, or X,, while the redundant inclusion of all four variables merely for the convenience of derivation. The key component is the edge Z; > Di, representing the self-masking dropout mechanism, i.e., whether a gene is undetected in a particular cell is determined by the gene's true expression level in the cell. This aligns with the several commonly acknowledged reasons of dropouts (Silverman et al. 2020 Jiing et al. 2022, Kim et al. 2020 Hicks et al.) 2018). But note that such self- masking assumption can also be relaxed, allowing for edges Zj Rs, discussed in $4.2
2.3 EXISTING MODELS AS SPECIFIC INSTANCES OF THE CAUSAL DROPOUT MODEL
In addition to the common biological understandings on the dropout mechanisms, many existing statistical models to handle dropouts also fit into our introduced framework specific parametric instances. Below, We give a few illustrative examples (detailed analysis is available in Appendix [B.T):
Example (Dropout with the fixed rates). Dj ~ Bernoulli(p)). i..., the gene Z, gets dropped out with a fixed probability p; across all individual cells. The representative models in this category
are the zero-inflated models (Pierson & Yau 2015 Kharchenko et al. 2014 Saeed et al., 2020
Yu et al. 2023 Min & Agresti| (005) and the hurdle models (Finak et al. 2015 Qiao et al. 2023)
Biologically, this can be explained by the random sampling of transcripts during library preparation,
regardless of the true expressions of genes. Graphically, in this case the edge Z. - D; is absent.
Example 2 (Truncating low expressions to zero). D, = 1Z, < c.). i.e., the gene Z, gets dropped out in cells whenever its expression is lower than a threshold c,. A typical kind of such truncation models
is, for simple statistical properties, the truncated Gaussian copula (Fan al. 2017 Yoon et al. 2020
Chung et al. (2022). Biologically, such truncation thresholds e; (quantile masking UJang et al.] 2022)
can be explained by limited sequeneing depths. Graphically, in this case the edge Z. D. exists.
Example 3 (Dropout probabilistically determined by expressions). D ~ Bernoulli(F.
a.)).
where B < O and F; is monotonically increasing (typically as CDF of probit or logistic (Cragg 1971: Liu 2004 Miao al. 2016)), i.e., gene Z; may be detected (or not) in every cell, while the higher it is expressed in a cell, the less likely it gets dropped out. Biologically, this can be explained by inefficient amplification. Graphically, the edge Z. + D; also exists, and non-deterministic.
As shown above, many existing statistical models to handle dropouts can fit into our proposed causal dropout graph as parametric instances. Furthermore, it is important to note that our causal model itself, as well as the corresponding causal discovery method proposed in ] are nonparametric. Hence, our model is more flexible and robust in practice. As is verified in $5.1 even on data simulated according to these parametric models, our method outperforms the methods specifically designed for them.
3

---

Published as a conference paper at ICLR 2024
2.4 POTENTIAL THEORETICAL ISSUES WITH THE IMPUTATION METHODS
While many parametric models exist ($2.) imputation methods have become more prevalent to handle dropouts in practice (Hou et al. 2020 Van Dijk et al. 2018 8 Li & Li 2018; 2 et al. 2018) Lopez et al. 2018). Imputation methods treat the excessive zeros as "missing holes" in the expression matrix and aim to fill them using nearby non-zero expressions. However, fundamental question arises: it theoretically possible to fill in these "missing holes"? In the context of the causal dropout model in Equation (T) we demonstrate that, in general, the answer is negative, as the underlying true distribution p () is unidentifiable from observations p(X) (detailed analysis in Appendix B.2):
Firstly, according the missing data literature, the underlying joint distribution (Z) is irrecoverable due to the self- dropout mechanism (Enders) 2022; Mohan et al.) 2013 Shpitser 2016):
Example Consider the following two models consistent with Example Miao et al. 2016):
1. Z; Exp(2), D. ~ (sigooiiliiiimiiid/ 2 Z,)). 2. Z ~ Exp(1), D; ~ Bernoulli ( ZZ log 2)).
With even the same parametric logistic dropout mechanism and different distributions of Zi. the resulting observations X; share exactly the same distribution, as 2e 22= T-clogz- e "4c"
Secondly, different from missing data, the "missing" entries here cannot be precisely located, i.e. D variables are latent. Therefore, even under the simpler fixed-dropout-rate scheme in Example T (dropout happens completely-at-rundee (CAR) (Little & Rubin 2019). ). (Z) remains irrecoverabl:
Example 5. Consider the following two models consistent with Example !
1. Z; (ernoulli(0.5), D; Bernoulli
2. Z, ~ Bernoulli D; ~ vernoulli(0.75)
With the CAR dropout mechanism and different p(Zz), the observed X; follow same Bernoulli0022)
2.5 MOTIVATION: PRESERVATION OF CONDITIONAL (IN)DEPENDENCIES IN NON-ZERO PARTS
We have demonstrated the general impossibility of recovering the true underlying distribution p(Z) from
observations X. However, for the purpose of GRNI, structure learning task, it is not necessary accurately recover p (Z). Alternatively, the graph structure among Z can be sufficiently recovered by the CI relations entailed by p(Z). Hence, the key question becomes whether the Cls of p( (Z) can survive the inestimable distortion introduced by dropouts and leave trace in p (X).
Example 6. Consider the linear structural equation model (SEM) consistent with a chain structure Zi + Zz Zg:
7, 7:
Zi
0.69
0,10
Zi = Eg
Zz Z1 + Eg
[8]
(b)
Z₃ = 3Zz + Eg
Z
F1, E2, Fg (0)
XaIX
X5. K:
and the logistic dropout scheme:
p- 0.28
0.00
D, Bernoulli (oaiiliiioiiiiiii-i1i )
for = 1,2.3. The scatter- plots (a) and (d) in Figure 2 demonstrate Z X Zg and X1 1 X3 marginally. Now condition on the middle variable.
(4)
The d-separation Z, 1 a Zz2 implies that Zi 1L Z3 Zz, as supported by (b) and (c) with
1: 1.1
PUXiA PiX L, - POGIN 1.N.-1
Desiity
Ze = 0, respectively. How- ever, in contrast, one observes X1 L X2) as supported by the scatters in (e): with the concentration of points along
Figure 2: Above: scatterplots of Z1: Z3 and X1; X3 under differ- ent conditions. Below: density plot of Xa under different condi- tions (vertical slices of scatters in (f)) to show X3 1 X = 1. Kernel width is set 0.02 to prevent from oversmoothing.
4

---

Published as a conference paper at ICLR 2024
axes (representing excessive zeros) and distorted, non-elliptical distribution away from the axes, X1 and Xg are still dependent when X2 = 0, because the zero entries of X2 mix different values of Zz. Nonetheless, this observation leads to an important insight: while zero entry of X2 may be noisy (i.e., may be technical), a non-zero entry of X2 must be accurate, i.e., biological. Although the non-zero parts undergo a distribution distortion (which is why imputation fails), this distortion does not impact the estimation of conditional independence estimation, which only cares about each specific value. That is, conditioning on X2 with any non-zero value æ2 + 0 is equivalent to conditioning on Z2 with the same 22 value, thereby eliminating the spurious dependence between X1 and X3. (f) shows the case where X2 = 1. Note that even though there are still dropout points of X1 and X3 on axes in (f), these respective dropouts are also independent given Zz. The identical density illustrates the conditional independence X1 1 X3|X2 =
plots (vertical slices of (f)) in Figure 2
CAUSAL DISCOVERY IN THE PRESENCE OF DROPOUT
Building upon the intuition of "conditioning on non-zero entries" in Example 6 we now formally introduce our deletion- CI test and the corresponding causal discovery methods.
Assumptions 1. We first list all the assumptions that may be used throughout this paper:
(AI.) Common assumptions including causal sufficiency, Markov condition, and faithfulness over Z U X UU R. acyclicity of G. and a pointwise consistent CI testing method.
(A2.) j E [p]. there is edge D; - Zz. i.e., dropout, technical artifact in the sequencing step, does not affect the genes' expressing process.
(A3.) Vi E [p]. D, has at most one parent, Z, i.c., a gene's 's dropout can only be directly affected by the expression value of itself, if not completely at random.
(A4.) For any variables A, E XUZ, C X Z. S C [p). A LL B|C, Rs = A f. BJC, Rs
O, i.e., dependencies conditioned on R are preserved in corresponding non-zero values. (A5.) For any conditioning set S C [p] involved in the required CI tests in the algorithm, p(Rs
0) > O, i.e., asymptotically, the remaining sample size is large enough for test power.
We will elaborate more on these assumptions' plausibility in $ Generally speaking, these assump- tions rather mild and are either theoretically testable empirically supported by real data.
3.1 CONDITIONAL INDEPENDENCE ESTIMATION WITH ZEROS DELETED
llsstrrtes how the conditional (in) dene)ddeeddcee relations may differ between true underlying expressions Z and dropout-contanne contomttoott X. Now let us first investigate the scenario where the dropout issue is disregarded, and CI tests are directly conducted on the complete data points of X:
Proposition (Bias to a denser graph). Assume (AI) (A2) Vi E [p].j E jp]. S c [p], we have Ze 1. 2, 1s X, 1 X, Xs. The reverse direction does not hold in general, and holds only when S Ø.
By Proposition D if one directly applies existing structure learning algorithms on observations X with- out dropout correction, the d-separation patterns entailed in the true GRN are generally undetectable. Consequently, the inferred graph tends to be much denser than the true one, posing significant bias, given that gene regulatory relations are usually sparse (Dixit et al.) 2016 Levine & Davidson 2005).
To address this bias, we introduce the following clean and essential finding: ignoring the data points with zero values for the conditioned variables does not change the CI relations. While the remaining samples exhibit different distribution from the true underlying distribution (as demonstrated in $2.4.. they maintain identical conditional independence relations. Formally, we state the following theorem:
Theorem 1 (Correct CI estimations). Assume
(AI)
(A2)
(A3) 0.
and (A4) For every
[p].j
[p]. S C [p], we have Z, 1 23 + X; XX, Zs.
Rs
Note that in Theorem [ the right hand side expresses the conditional independence relations as X; 1 X, Zs. Rs 0. rather than Z; 1 Z, Zs, Rs = 0. This is because we only delete the samples with zeros for conditioned variables, while the zeros in estimands X, and X, are retained, and the underlying Z; and Z, remain unobservable. Generally speaking, for any Z; explicitly involved in a CI estimation, samples with the corresponding zeros must be deleted, i.e., R; = 0 must be conditioned on. Alternatively, one may wonder what if the samples with zeros in all variables are deleted, not just the conditioning set. Interestingly, under assumption [A3] both approaches are asymptotically equivalent, meaning that the CI relations can still be accurately recovered: Z; 1
5

---

Published as a conference paper at ICLR 2024
Zj|Zg Z; 1 Zz R.jus = O, though empirically, the latter may delete more samples and potentially reduce the statistical power of the test. It is also worth noting that in a more general setting where iA dropped. the two approaches not equivalent, and the latter provides tighter bound for estimating conditional independencies. Further details will be discussed in Theorem 3
3.2 CAUSAL DISCOVERY METHODS WITH DROPOUT CORRECTION
Now, we integrate the above consistent CI estimation into established causal discovery methods:
Definition 1 (The general procedure for causal discovery with dropout correction). Perform any consistent causal discovery algorithm (c.g., PC (Spirtes et al. (0000) based on CI relations estimated by test-wise delction as in Theorem Infer Z; رزبة Zg if and only if X; 1 X, Zg, Rg = 0, and use these d-separation patterns to inter the graph structure among Z.
Constraint- methods. The test-wise deletion can be seamlessly incorporated into any existing constraint- based methods for structure learning, as shown in Definition [ Assumption [A5) ensures that asymptotically, for each CI relation to tested, the remaining sample size after lest-wise deletion remains infinite, and the consistency of the method is still guaranteed.
Greedy score-based methods. While the proposed GRNI on dropout data with integrated test-wise deletion in constraint-based methods is asymptotically consistent, empirical reliability may be limited due to order-dependency and error propagation of constraint- 2014). In contrast, ccore-based methods generally search over the graph space and is immune to order-dependency, and thus may provide more empirically accurate results on large scale GRNs.
methods (Colombo & Maathuis
One typical bosre-based method is the Greedy Equivalence Search (GES (Chickering 2002 algo- rithm. GES uses a scoring function to assign score to each directed acyclic graph (DAG) given data, and finds the optimal score by traversing over the space of CPDAGs. Consisting of forward phase and a backward phase, GES iteratively performs edge additions and deletions in two phases respectively to optimize the total score until further improvement is not possible.
Score functions are typically assumed to exhibit three attributes: global consistency, local consistency, and decomposability. However, a closer examination on Chickering (2002)'s proof reveals that for GES's s consistency, only a locally consistent score is needed. As another way of rendering faithfulness, local consistency enforces adding or deleting edges based on CI relations. This insight, as is also affirmed from the connection between BIC score and Fisher-Z test statistics (Lemma 5.1 of (Nandy) al. 2018)) and further echoed in (Shen et al. 2022), suggests that GES can be completely reframed as constraint-based method using CI tests without defined score. Thus, our deletion-basee CI test can be easily incorporated into GES, ensuring asymptotic consistency and order-independence.
ON TESTABILITY OF THE CAUSAL DROPOUT MODEL
In 3 we propose principled approach for GRNI scRNA-seq data with dropouts, i.e., deleting samples with zeros for the conditioned variables in each CI test. The asymptotic correctness this
approach depends assumptions listed in Assumptions I However, one may wonder whether these assumptions on the dropout mechanisms hold in practice. Therefore, in this section, we delve deeper into these assumptions and show that they can be either theoretically empirically verified from data.
4.1 OVERVIEW OF THE ASSUMPTIONS
200
To enhance clarity, we first provide a detailed explanation of each assumption. The common assumptions for causal discovery in (A7) are thoroughly discussed in (Spirtes et al. 2000). and thus we focus more on the remaining four. Specifically. [A2) [A2] and J structural assumptions on the dropout mechanism, indicating that whether a gene dropped out or not in the sequencing procedure (temporally later) does not affect the genes' expressions and interactions (temporally earlier), and that each gene's dropout can only be directly affected by the expression value of itself, not by other genes' dropouts or expression values. [A4] ] is called faithful observability in previous work on missing data (Strobl et al. 2018 TTu et al.) 2019): while conditional independence means independence everywhere (on every conditioned value), conditional dependence only requires
Sample size with Re-
6030
0:30
2000
2000
2ic 400 aid Abn
ude 220 ibe
The number of tess PC
Figure 3: The remaining sam- ple size after deleting zero con- ditioning samples, for each CI test during a run of PC real data (Dixit et al.) 2016).. with 15 genes and 9843 cells.

---

Published as a conference paper at ICLR 2024
dependence at some conditioned value. Thus (A4) is assuming that any dependence conditioned on R are preserved in corresponding non-zero values. This is reasonable as biological regulation is usually performed by genes' expression instead of non-expression. and is unlikely that all the samples that capture the dependence happen to be dropped out. Lastly, (A5] assumes that for each required CI test during the algorithm, asymptotically there exists sufficient number of samples for the test power. This assumption is empirically validated in real data, as shown in Figure 3 and also validated by a synthetic experiment with varying dropout rates (Figure 12 in Appendix D.I
4.2 VALIDATING AND DISCOVERING THE STRUCTURAL MECHANISM ON DROPOUT
Among the assumptions discussed above, the structural assumption [A3) concerning the dropout
mechanism is of particular interest: each gene's 's dropout solely affected by the gene itself and not by any other genes. While this assumption is in line with most of the existing parametric models ($2.)
and plausible in the context of scRNA-seq data, where mRNA molecules are everse-transeribe
independently across individual genes, it is still desirable to empirically validate this assumption
based on the data. Therefore, here we propose a principled approach to validate this assumption. In other words, in addition to discovering the relations among genes, we also aim to discover the inherent mechanism to explain the dropout of each gene from the data.
Surprisingly, have discovered that even without assumption [A3] it is still possible to identify the regulations among Z and the dropout mechanisms for R. The approach is remarkably simple: we perform causal structure search among Z using the procedure outlined in Definition with only one modification: instead of deleting samples with zeros only for the conditioned variables, we delete samples with zeros for all variables involved in the CI tests. By doing so, the causal structure among Z, representing the GRN, as well as the causal relationships from Z to R. representing the dropout mechanisms, can be identified up to their identifiability upper bound. Formally, we have Definition 2 (Generalized GRN and dropout mechanisms discovery). Perform the procedure outlined in Definition D except for inferring Z; L. [Zg if and only if 7, 1 Z Rsu(is) = 0.
Theorem 2 (Identification of GRN and dropout mechanisms). Assume (AT) In the CPDAG among Z output by Definition2 2 if Zi and Zj are non-adjacent, it implies that they are indeed non-adjacent in the underlying GRN, and for the respective dropout mechanisms, Zi does not cause R, and Zj does not cause R₃. On the other hand, if Zí and Zj are adjacent, then in the underlying GRN Z and Z, are non non-adjace adjacenj only in one particular case, and for dropout mechanisms, the existence of the causal relationships Z; + R3 and Z, - R² must be naturally unidentifiable.
(A2) (A4) and [A5)
Due to space limit, here we only give the above conclusion. For illustrative examples, proofs, and detailed elaboration on the identifiability in general cases without A3] please refer to Appendix C
5 EXPERIMENTAL RESULTS
In this section, we conduct extensive experiments to validate our proposed method in $ demon- strating that it is not only theoretically sound, but also leads to superior performance in practice. We provide the experimental results linear SEM simulated data, more 'rcalistic' synthetic and curated
scRNA-seq data, and real-world experimental data in $5.1 $5. and $5.3 respectively!
5.1 LINEAR SEM SIMULATED DATA
Methods and simulation setup. We assess our method (testwise delction) by comparing it to other dropout-handling approaches on simulated data, including MAGIC Van Dijk al. (2018) for imputa- tion, mixedCCA Yoon et al. (2020) as parametric model, and direct application of algorithms on full samples (without deleting or processing dropouts, corresponding to Proposition]] I.. We also report the results on the true underlying Z, denoted as Oracle*, to examine the performance gap between these methods and the best case (but rather impossible) without dropouts. PC PC[Stisses et al. (2000) and GES (Chickering] 2002) with FisherZ test are used as the base causal discovery algorithms. We randomly generate ground truth causal structures with p E {10, 20, 30} nodes and degree of to 6. Accordingly. Z is simulated following random linear SEMs in two cases: jointly Gaussian, and Lognormal distri- butions to better model count data (Bengtsson et al. 2003). We apply three different types of dropout mechanisms on Z to obtain X: (1) dropout with the fixed rates, (2) truncating low expressions to zero, and (3) dropout probabilistically determined by expression, which are described in Examples ] to 3 respectively. We report the structural Hamming distance (SHD) wrt. the true and estimated
CPDAGs calculated from 5 random simulations. More details are available in Appendix D.I
Codes available at https: / htt//gttuuuu wtwpth..ccco/aaaaaaaaaaaaaaaaaaa Caaaa--aaaaaaa -orrg

---

Published as a conference paper at ICLR 2024
Ornele*
Testwise deletion
Full sumples
MAGIC
mixedCCA
Fleed dioppsst salo
200
100
SID ofCPDNG (lover io bette))
Truncated
200
200
100
100
0
Degree
Degree
(a) |2] = 30; Gaussian.
(b) |Z| = 30; Lognormal.
Figure 4: Experimental results (SHDs of CPDAGs) of 30 variables on simulated data, where three dropout mechanisms are considered. The variables Z follow Gaussian Lognormal distribution.
Systhefic
Curated
F00s 50 8050 2005 30 3555 Gopd
GB8
SINCERITIES
15
TPCOR
MIC
GENBOOST22
SCODE
Figure 5: Experimental results (Fl-scores of estimated skeleton edges) BoolODE simulated data and BEELINE benchmark framework (Pratapa et al. 2020). The 9 rows correspond to PC, GES, and 7 other GRNI-specific SOTA algorithms benchmarked. The 10 colored column blocks correspond to all synthetic and 4 curated datasets (Pratapa et al. 2020). The 5 column strips in each block correspond different dropout-handling strategies. Cell colors indicate the corresponding values (brighter higher, i.e., better). The maximum and minimum of each strip are annotated.
Experimental results. The results 30 nodes are provided in Figure Our method (1) consistently outperforms other competitors across most settings, particularly when using PC the base algorithm, and (2) leads to performance close to the oracle. Directly applying algorithms on full samples often performs poorly due to potential false discoveries as suggested by Proposition]] [ Notably, mixedCCA performs relatively well in the Gaussian and truncated dropout setting (as parametric model with no model misspecification) but still outperformed by our method. This further validates our method's effectiveness and flexibility. For more experimental results, including results in different scales, on different metrics, and with varying dropout rates, please kindly refer Appendix D.I
5.2 REALISTIC BOOLODE SYNTHETIC AND CURATED DATA
While synthetic data is preferred for assessment thanks ground-truth graphs, to truly gauge real- world applicability, we aim to conduct simulations that resemble scRNA-seq data more closely than just linear SEM. Furthermore, integrated into algorithms specifically designed for scRNA-seg data, rather than just general causal discovery methods like PC and GES. To accomplish this, conduct array of experiments based
seek investigate the efficacy of our zero-deletion approach when
8

---

Published as a conference paper at ICLR 2024
on the widely recognized BEELINE framework (Pratapa et al., 2020), which offers scRNA-seg data simulator called BoolODE, and benchmarks various SOTA ERNI-specifee algorithms.
Following the BEELINE paradigm, for data simulation, simulates gene expressions with pseudotime indices from Boolcan regulatory models. We simulate all the 6 synthetic and 4 literature-curated datasets in (Pratapa et al. 2020), each with 5,000 cells and a 50% dropout rate, following all the default
the BoolODE simulator that basically
For algorithms, in addition to PC and GES, we examined all SOTA algorithms (if executable; with default nypper-prameters) benchmarked in (Pratapa et al., 2020). For each algorithm, five different tropout-handling strategies are assessed, namely, oracle*, testwise deletion, full samples, imputed, and binarization (Qiul 2020). Method details, e.g., how is zero deletion incorporated into various algorithms, can be found in Appendix D.2
The Fl-scores of the estimated skeleton edges are shown in Figure [ We observe that: 1) Dropouts do harm to GRNI. Among the 90 dataset-algorithm pairs (9 algorithms x 10 datasets), on 65 of them (72%), 'full samples' (with dropouts) performs worse than 'oracle'. 2) Existing propout-handling strategies (imputation and binarization) don't work well. Among the 90 pairs, 'imputed' and 'binary' even worse than 'full samples' on 45 (50%) and 58 (64%) of them, respectively, i.e., they may even be counterproductive. As discussed in 2.4 such strategies may indeed introduce additional bias. 3) The proposed zero-deletion is effective in dealing with dropouts, with consistent benefits across different integrated algorithms and different datasets. Among the 90 pairs, 'zero-deletion' is better than 'ull samples' and than 'imputed' on 64 (71%) and 61 (68%) of them, respectively, i.e., it effectively helps GRNI with dropouts, and is more effective than imputation in dropout-handling
5.3 FEAL-WORLDD EXPERIMENTAL DATA
To examine the efficacy of our approach in real-world setting, we applied our method on gene expression data collected via single-cell Perturb-seq (Dixit et al. 2016).
As in prior work Saeed et al. 20 Dixit et al. 2016), we focus on the bone marrow dendritic cells (BMDC)
in the unperturbed setting (9843 cells in total), and on 21 transcription factors believed to be driving differentiation
in this tissue type Dixit et al. (0166): Garber et al. (2012).
We compared the predicted causal edges with known relations derived as union of statistically significant
predictions from interventional data Dixit et al.) (2016) and prior knowledge Han et al. (2018). For this amalysis, we ran a baseline of the PC algorithm with all the samples (PC-full), compared to the PC algorithm with testwise deletion (PC-test-del). From the figure, one can appreciate that PC-full introduces numerous edges that are not backed by prior knowledge (labeled red in Figure 6). For example, the E2f1 - Cebpbp causal link is opposite of what has been reported in the literature Gutsch al. (201T) T) Furthermore, Stat3 and E2f1 are known to cooperate in downstream regulation, rather than E2FT being a cause of Stat3 Hutchins et al. (2013). On the contrary, majority of the edges that testwise deletion provides are previously known (labeled in blue), with a few that are not in the list of previously seen interactions that we used (labeled in black). For experimental details and results on another two real-world experimental datasets, please kindly refer to Appendix D.3
Figure 6: Experimental results on 21 key genes from Pertub-seq data (Dixit et al. 2016). Red edges are returned by PC-full but not by PC-test-del. Blue and black edges are returned by PC-test-del. Specif- ically, blue edges are the regulatory inter- actions priorly known (Dixit et al. 2016). while black edges are not priorly known.
6 CONCLUSION AND DISCUSSIONS
Building upon common understanding of dropout mechanisms, we develop the causal dropout model to characterize these mechanisms. Despite the non-ignorable observed zeros resulting from dropouts, we develop testwise-deletie procedure to reliably perform CI test, which can be seamlessly inte- grated into existing causal discovery methods to handle dropouts and is asymptotically correct under mild assumptions. Furthermore, our causal dropout model serves as systematic framework to verify the qualitative mechanism of dropout studied in the literature is valid, and learn such a mechanism from observations. Extensive experiments method leads to improved performance in practice. A possible limitation is the decreasing sample size after testwise deletion, and future work includes developing practical method to resolve it.
simulated and real-world datasets demonstrate that our

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGEMENTS
This material is based upon work supported by the AT Research Institutes Program funded by the National Science Foundation under AI Institute for Societal Decision Making (AI-SDM), Award No. 2229881. The project is also partially supported by the National Institutes of Health (NIH) under Contract ROIHLI59805, and grants from Apple Inc., KDDI Research Inc., Quris AI, and Infinite Brain Technology. Petar Stojanov was supported in part by the National Cancer Institute (NCI) grant number: K99CA277583-00 and funding from the Eric and Wendy Schmidt Center at the Broad Institute of MIT and Harvard.
REFERENCES
Britt Adamson, Thomas M Norman, Marco Jost, Min Y Cho, James K Nuñez, Yuwen Chen, Jacqueline E Villalta, Luke A Gilhert, Max A Horlbeck, Marco Y Hein, al. multiplexed single-cell CRISPR screening platform enables systematic dissection the unfolded protein response. Cell, 167(7):1867-188, 2016.
Bruce Alberts, Alexander Johnson, Julian Lewis, Martin Raff, Keith Roberts, and Peter Walter. Molecular biology of the cell. Garland Science, 2002.
Matthew Amodio, David Van Dijk, Krishman Srinivasan, William S Chen, Hussein Mohsen, Kevin R Moon, Allison Campbell. Yujiao Zhao, Xiaomei Wang, Manjunatha Venkataswamy, et Exploring single-cell data with deep multitasking neural networks. Nature Methods, 16(11);1139) 11:11199-1145 2019.
Tallulah Andrews and Martin Hemberg. False signals induced by single-cell imputation. F1000Research, 2018.
Tallulah Andrews and Martin Hemberg. M3Drop: dropout-based feature selection for scRNA-seq. Bioinfor- matics, 35(16)22865 2867, 2019.
Anastasiya Belyaeva, Chandler Squires, and Caroline Uhler. DCT: learning causal differences between gene regulatory networks. Bioinformatics, 37(18):3067-3069 2021.
Martin Bengtsson, Anders Stählherg, Patrik Rorsman, and Mikael Kubista. Gene expression profiling in single cells from the pancreatic islets Langerhans reveals lognormal distribution of mRNA levels. Genome Research, 15(10):1388-133 2005.
Rebecca Boiarsky, Nicholas Haradhvala, Sann-Baptiste Alberge, Romanos Sklavenitis- Pistofidis, Tarek Mouhieddine, Oksana Zavidij, Ming-Chieh Shih, Danielle Firer, Mendy Miller, Habib El-Khoury, et al. Single cell characterization of myeloma and its precursor conditions reveals transcriptional signatures of early tumorigenesis. Nature Communications, 13(1):7040, 2022.
David Maxwell Chickering. Optimal structure identification with greedy search. Journal of Machine Learning Research, 3(Nov):507-554 2002.
Junsouk Choi, Robert Chapkin, and Yang Ni. Bayesian causal structural learning with zero-inflated Poisson Bayesian networks. Advances in Neural Information Processing Systems, 33:5887-5897,, 2020.
Li-Fang Chu, Ning Leng, Jue Zhang, Zhonggang Hou, Daniel Mamott, David T Vereide, Jeea Choi, Christina Kendziorski, Ron Stewart, and James A Thomson. Single-cell RNA-seq reveals novel regulators of human embryonic stem cell differentiation definitive endoderm. Genome Biology, 17:1-20, 2016.
Hee Cheol Chung, Yang Ni, and Irina Gaynanova. Sparse semiparametrie discriminant analysis for high- dimensional zero-inflated data, arXiv preprint arXiv:2208.0372 2022.
Diego Colombo and Marloes H Maathuis. Order-inde constraint- based causal structure learning. Journal of Machine Learning Research, 15(1):3741-3782, 2014.
The Tabula Muris Consortium, Overall coordination, and Logistical coordination. et al. Single-cell transcrip- tomies of 20 mouse organs creates Tabula Muris. Nature 562(7727):367 372, 2018.
The Tabula Sapiens Consortium, Robert C Jones, Jim Karkanias, Mark A Krasnow, Angela Oliveira Pisco, Stephen R Quake, Julia Salzman, Nir Yosef, Bryan Bulthaup. Phillip Brown, et al. The Tabula Sapiens: multiple-organ, single-cell transcriptomic atlas of humans. Science, 3abb6594),.bb48866 2022.
John G Cragg. Some statistical models for limited dependent variables with application to the demand for durable goods. Econometrica: Journal of the Econometric Society, pp. 229,844, 1971.
10

---

Published as a conference paper at ICLR 2024
Haoyue Dai, Peter Spirtes, and Kun Zhang. Independence testing-based approach causal discovery under measurement error and linear non-gaussian models. Advances Neural Information Processing Systems, 35: 77524-27536 2022.
Jiarui Ding, Xian Adiconis, Sean Simmons, Monika Kowalczyk, Cynthia Hession, Nemanja Marjanovie, Travis K Hughes, Mare H Wadsworth, Tyler Burks, Lan T Nguyen, et al. Systematic comparison of single-cell and single-nucleus sequeencgn methods. Nature Biotechnology, 38(6):737-746, 2020.
Atray Dixit, Oren Parnas, Biyu Li, Jenny Chen, Charles P Fulco, Livnat Jerby-Aron, Nemanja D Marjanovie, Danielle Dionne, Tyler Burks, Raktima Raychowdhury, et al. Perturb-Sec Seq: dissecting molecular circuits with scalable single-cell RNA profiling of pooled genetic screens, cell, 167(7):1853-1866 2016.
Craig K Enders. Applied missing data analysis. Guilford Publications, 2022.
Paul Erdös, Alfred Rényi, al. On the evolution of random graphs. Publication of the Mathematical Institute of the Hungarian Academy of Sciences, 5(1):17-60, 1960.
Jianqing Fan, Han Liu, Yang Ning, and Hui Zou. High dimensional semiparametric latent graphical model for mixed data. Journal of the Royal Statistical Society. Series (Statisticul Methodology), pp. 405-421, 2017.
Greg Finak, Andrew McDavid, Masanao Yajima, Jingyuan Deng, Vivian Gersuk, Alex K Shalek, Chloe Slichter, Hannah W Miller, M Juliana McElrath, Martin Prlic, et MAST: flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Riology, 16(1 16,3, 2015.
Yu Fu, Pei-Hsuan Wu. Timothy Beane, Phillip D Zamore, and Zhiping Weng. Elimination of per duplicates in RNA-seq and small RNA-seq using unique molecular identifiers. BMC Genomics, 19:1-14, 2018.
Wayne A Fuller. Measurement error models. John Wiley & Sons, 2009.
Erdun Gao. Ignavier Ng, Mingming Gong, Li Shen, Wei Huang, Tongliang Liu, Kun Zhang, and Howard Bondell. MissDAG: Causal discovery in the presence of missing data with continuous additive noise models. Advances Neural Information Processing Systems, 35:5024 5038, 2022.
Manuel Garber, Nir Yosef, Alon Goren, Raktima Raychowdhury, Anne Thielke, Mitchell Guttman, James Robin- son, Brian Minic, Nicolas Chevrier, Zohar Itzhaki, et al. high-throughput chromatin immunoprecipitationn approach reveals principles of dynamic gene regulation in mammals. Molecular Cell, 47(5):8008822, 2012.
Romina Gutsch, Judith D Kandemir, Daniel Pietsch, Christian Cappello, Johann Meyer, Kathrin Simanowski, René Huber, and Korbinian Brand. CCAAT/enhancer binding protein 3 inhibits proliferation in monocytic cells by affecting the retinoblastoma protein/E2F/eyclin E pathway but is not directly required for macrophage morphology. Journal of Biological Chemistry, 286(26) 201
Heonnong Han, Jae-Wor Won Cho. Sangyoung Lee, Ayoung Yun, Hyojin Kim, Dasom Bae, Sunmo Yang. Chan Yeong Kim, Muyoung Lee, Eunbeen Kim, et al. TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Research, 46(D1):D380-DD DD36, 2018.
Stephanie C Hicks, F William Townes, Mingxiang Teng, and Rafael A Trizarry. Missing data and technical variability in single-cell RNA-sequencing experiments. Biostaristics, 19(4):562-578, 2018.
Wenpin Hou, Zhicheng Ji, Hongkai Ji, and Stephanie C Hicks. A systematic evaluation of single-cell RNA- sequencing imputation methods. Genome Biology, 21:1-30, 2020.
Alan Huang, Levi A Garraway, Alan Ashworth, and Barbara Weher. Synthetic lethality as an engine for cancer drug target discovery. Nature Reviews Drug Discovery, 19(1):3338, 2020.
Mo Huang, Jingshu Wang, Eduardo Torre, Hannah Dueck, Sydney Shaffer, Roberto Bonasio, John I Murray, Arjun Raj. Mingyao Li. and Nancy R Zhang. SAVER: gene expression recovery for single-cell RNA sequencing. Nature Methods, 15(7)5339 542, 2018.
Andrew Paul Hutchins, Diego Diez, Yoshiko Takahashi, Shandar Ahmad, Ralf Jauch, Michel Lucien Tremblay, and Diego Hiranda-Saaveddra Distinct transcriptional regulatory modules underlie STAT3's cell type- independent and cell type-specific functions. Nucleic Acids Research, 41(4):2155-217) 2013.
Takahiro Ito, Michael J Young. Ruitong Li, Sidharth Jain, Andreas Wernitznig, John M Krill-Burger, Burger, Christo- pher T Lemke, Davide Monducci, Diego J Rodriguez. Liang Chang, et al. Paralog knockout profiling identifies DUSP4 and DUSP6 a digenic dependence in MAPK pathway-driven cancers. Nature Genetics, 53(12): 1664-1672, 2021.
11

---

Published as a conference paper at ICLR 2024
Ruochen Jiang, Tianyi Sun, Dongyuan Song. and Jingyi Jessica Li. Statistics biology: the zero-inflation controversy about scRNA-seq data. Cenome Biology, 23(1):1-24 2022.
Peter Kharchenko, Lev Silherstein, and David T Scadden. Bayesian approach to single-cell differential expression analysis. Nature Methods, 11(7):740-742, 2014.
Hyun Kim, Xiang Zhou, and Mengjie Chen. Demystifying "drop-outs" in single-cell UMI data. Genome Biology, 21(1):196, 2020.
Manabu Kuroki Judea Pearl. Measurement bias and effect restoration causal inference. Biometrika, (2):423-437, 2014.
Michael Levine and Eric Davidson. Gene regulatory networks for development. Proceedings of the National Academy of Sciences, 102(14):4936 4942, 2005.
Wei Vivian Li and Jingyi Jessica Li. An accurate and robust imputation method scimpute for single-cell RNA-seq data. Nature Communicationss 9 91):9977, 2018.
Roderick JA Little and Donald B Rubin. Statistical analysis with missing data, volume 793. John Wiley & Sons, 2019.
Chuanhai Liu. Robit regression: a simple robust alternative to logistic and probit regression. Applied Bayesian Modeling and Casual Inference from Incomplete- Damapeeee-DDta Perspectives, pp. 227-238, 2004.
Romain Lopez, Jeffrey Regier, Michael B Cole, Michael Jordan, and Nir Yosef. Deep generative modeling for single-cell transcriptomies. Nature Methods, 15(12):0053-1058, 2018.
Wang Miao, Peng Ding, and Zhi Geng, Identifiability of normal and normal mixture models with nonignorable missing data. Journal of the American Statistical Association, (11(516):6733 1683, 2016.
Yongyi Min and Alan Agresti. Random effect models for repeated measures of zero-inflated count data. Statistical Modelling, 5(1):1-19, 2005.
Karthika Mohan. On handling self-masking and other hard missing data problems. In AAAI Symposium, 2018.
Karthika Mohan, Judca Pearl, and Jin Tian. Graphical models for inference with missing data. Advances in Neural Information Processing Systems, 26, 2013.
Preetam Nandy, Alain Hauser, and Marloes H Maathuis. High-dimentn sich-dimmnsionall consistency in score- score-based and hybrid structure learning. The Annals of Statistics, 46(6A):3151 3183, 2018.
Cyril Neftel, Julie Laffy, Mariella G Filbin, Toshiro Hara, Marni E Shore, Gilbert Rahme, Alyssa R Richman, Dana Silverbush, McKenzie L Shaw, Christine M Hebert, al. An integrative model of cellular states, plasticity, and genetics glioblastoma. Cell, 178(4):835-849, 2019.
Phoebe CR Parrish, James D Thomas, Austin M Gahel, Shriya Kamlapurkar, Robert K Bradley, and Alice H Berger. Discovery synthetic lethal tumor suppressor paralog pairs the human genome. Cell Reports, 109597, 2021.
Judea Pearl. Causality: models, reasoning and inference. Cambridge University Press, 2009.
Judea Pearl. On measurement bias causal inference. arXiv preprint urXiv:1203.3504 2012.
Emma Pierson and Christopher Yau. ZIFA: Dimensionality reduction for zero-inflated single-cell gene expression analysis. Genome Biology, 16(1):1-10, 2015.
Aditya Pratapa, Amogh P Jalihal, Jeffrey N Law. Aditya Bharadwaj. and TM Murali. Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data. Nature Methods, 17(2):147-154, 2020.
Zhili Qiao, Elle Barnes, Susannah Tringe, Daniel P Schachtman, and Peng Liu. Poisson hurdle model-based method for clustering microbiome features. Bioinformatics, 39ta) 2023.
Peng Qiu. Embracing the dropouts in single-cell RNA-seq analysis. Nature Communications, 11(1):1 166, 2020.
Joseph Ramsey, Madelyn Glymour, Ruben Sanchez-Romero. and Clark Glymour. A million variables and more: the fast greedy equivalence search algorithm for learning high sigh--imensionll graphical causal models, with an application to functional magnetic resonance images. International Journal of Data Science and Analytics, 3: 121-129, 2017.
12

---

Published as a conference paper at ICLR 2024
Aviv Regev, Sarah A Teichmann, Eric S Lander, Ido Amit, Christophe Benoist, Ewan Birney, Bernd Bodenmiller, Peter Campbell, Piero Carninci, Menna Clatworthy, The human cell atlas. elife, 6:c27041, 2017.
Paul R Rosenbaum. Overt hias in observational studies. Springer, 2002.
Basil Saeed, Anastasiya Belyaeva, Yuhao Wang, and Caroline Uhler. Anchored causal inference in the presence of measurement error. In Proceedings of the Thirty-Sixth Conference on Uncertainty in Artificial Intelligence, pp. 619-628. AUAI Press, 2020.
Richard Scheines and Joseph Ramsey. Mcasurement error and causal discovery, In CEUR procedding volume 1792, pp. 1. NIH Public Access, 2016.
Xinwei Shen, Shengyu Zhu, Jiji Zhang, Shoubo Hu, and Zhitang Chen, Reframed GES with neural conditional dependence measure. In Uncertainty in Artificial Intelligence, pp. 1782-1791. PMLR, 2022.
Ilya Shpitser, Consistent estimation of functions of data missing non mon-ooootonicall and not at random. Advances in Neural Information Processing Systems, 29, 2016.
Justin D Silverman, Kimberly Roche, Sayan Mukherjee, and Lawrence A David. Naught all zeros in sequenee count data are the same. Computational and Structural Biotechnology Journal, 18:2789-2798, 2020.
Peter Spirtes. Clark N Glymour. Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Eric V Strobl, Shyam Visweswaran, and Peter L. Spirtes. Fast causal inference with non-random missingness by test-wise deletion. International Journal of Data Science and Analytics, 6:47-62, 2018.
Valentine Svensson, Roser Vento-Tormo, and Sarah A Teichmann. Exponential scaling of single-cell RNA-seq in the past decade. Nature Protocols, 13(4):599-604 2018.
Wenhao Tang, François Bertaux, Philipp Thomas, Claire Stefanelli, Malika Saint, Samuel Marguerat, and Vahid Shahrezaci. bayNorm: Bayesian gene expression recovery. imputation and normalization for single-cell 2020.
RNN-ssqqeniing data. Bioinformatiic,
Ruibo Tu, Cheng Zhang, Paul Ackermann, Karthika Mohan, Hedvig Kjellström, and Kun Zhang. Causal discovery in the presence missing data, In The 22nd International Conference on Artificial Intelligecce and Statistics, pp. 1762-1770. PMLR, 2019.
Po- Yuan Tung, John Blischak, Chiaowen Joyce Hsiao, David Knowles, Jonathan E Burnett, Jonathan Pritchard, and Yoav Gilad. Batch effects and the effective design single-cell gene expression studies. Scientific Reports, 7(1):39921,001 2017.
David Van Dijk, Roshan Sharma, Juozas Nainys, Kristina Yim, Pooja Kathail, Ambrose J Carr, Cassandra Burdziak, Kevin R Moon, Christine L. Chaffer, Diwakar Pattabiraman, et al. Recovering gene interactions from single-cell data using data diffusion. Cell, 174(3):716 729, 2018
Yuhao Wang, Liam Solus, Karren Yang, and Caroline Uhler. Permutation-basee causal inference algorithms with interventions. Advances in Neural Information Processing Systems, 30, 2017.
Mekuanint Simeneh Workie and Ahehaw Gedef Azene. Bayesian negative binomial logit hurdle and zero-inflated model for characterizing smoking intensity. Journal of Big Data, 8:1-20, 2021.
Feiyi Xiao, Junjie Tang, Huaying Fang, and Ruibin Xi. Estimating graphical models for count data with applications to single-cell gene network. Advances in Neural Information Processing Systems, 35:29038- 29050, 2022.
Grace Yoon, Raymond Carroll, and Irina Gaynanova. Sparse semiparametric canonical correlation analysis for data of mixed types. Biometrika, 107(3):009-625, 2020.
Shiqing Yu, Mathias Drton, and Ali Shojaie. Directed graphical models and causal discovery for zero-inflated data. Conference on Causal Learning and Reasoning, pp. 27-67. PMIR, 2023.
Luke Zappia, Belinda Phipson, and Alicia Oshlack. Splatter: simulation of single-cell RNA sequencing data. Genome Biology, 18(1):174, 2017.
Jiaqi Zhang, Chandler Squires, and Caroline Uhler. Matching desired causal state via shift interventions. Advances Neural Information Processing Systems, 34, 2021.
Jiagi Zhang, Louis Cammarata, Chandler Squires, Themistoklis P Sapsis, and Caroline Uhler. Active learning for optimal intervention design in causal models. Nature Machine Intelligence, pp. 1-10, 2023.
13

---

Published as a conference paper at ICLR 2024
Kun Zhang. Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kemel-based conditional independence and application in causal discovery. Proceedings of the Twenty-Seventh Conference Uncertainty Artificial Intelligence, pp. 804-81 AUAI Press, 2011.
Kun Zhang, Biwei Huang. Jiji Zhang, Clark Glymour, Bernhard Schölkopf. Causal discovery from honstationaryyheeerggeennnu data: Skeleton estimation and orientation determination. In IJCA:: Proceedings of the Conference, volume 2017, pp. 1347. NIH Public Access, 2017.
Yujia Zheng, Biwei Huang. Wei Chen, Joseph Ramsey, Mingming Gong, Ruichu Cai, Shohei Shimizu, Peter Spirtes, and Kun Zhang. Causal-learn Causal discovery in python, arXiv preprint 0tt:::2377.110000 2023.
14

---

Published as a conference paper at ICLR 2024
Appendix
Table of Contents
A Proofs of Main Results
15
A.1 Proof of Proposition
15
A.2 Proof of Theoremil
15
A.3 Proofs of Results in
16
Discussions
Existing Parametric Models Instances of the Causal Dropout Model
B.2 Relation to Imputation Methods and Missing Data Problem
a General Causal Dropout Model with Relaxed Assumptions
Elaboration on Theoorm 3 with an Illustrative Example
19
Motivation on the General Causal Dropout Model
20
Proofs of the Results in the General Causal Dropout Model in $4.2
21
D
Supplementary Experimental Details and Results
22
D
Linear SEM Simulated Data
22
D.2 Realistic Synthetic and Curated Data under the BEELINI INE Framework
27
D.3 Real-World Experimental Data
28
A PROOFS OF MAIN RESULTS
A.1 PROOF OF PROPOSITION
!
Proposition (Bias to denser graph). Assume (AI) Z, Zs = X; f Xj Xs. The reverse direction does not hold in general, and holds only when S Ø
(A2)
Vi [p.je [p].S [p], we have Z4 1
Proof of Proposition ! The = direction: By Markov condition in (A7) we have Z; 1 Z, Zg = Zz, by (A2) all non-colliders not Xs
Zi LaZz Zzs. Consider the open path Zi and all colliders have descendants in Xs, and thus X; <- Z.
Zj + X2 is still open, i.e.,
X, XXx, By faithfulness in (A7) X, XxX, X, I X
The < direction doesn't 't hold: Because the non-colliders in Zg cannot be conditioned children Xs. X, Xxs only when every path connecting Z; and Z, has collider, i.e., Z, 1 Zj.
A.2 PROOF OF THEOREM
Theorem 1 (Correct CI estimations). Assume (A]
AA2)
(A3)
and (A4) For every [pl.j
[p]. S C [p), we have Z, Zs + X; X, Zs. Rs 0.
Proof of Theorem 7 The + direction: By faithfulness in (A/) Z, 2, = 7 L By (A3) the blocked paths remain blocked for the two children, i.c., X, LZs, Since conditioning on Rs does not introduce any collider, we further have X; 1,X, /Zs, Rs. By Markov condition in
(A)) X; X2 Zs, Rs X, 1 X,|Zs, R.s, which implies X; 1 X, Zg, Rs 0.
The + direction: Show by its contrapositive. By Markov condition in (AI) Z, Z = Z; Lz IZs. By (A2), Rg does not contain any non-collider the open path connecting Z; and Zj. and thus still Z, L.Z, Rs. By (A3) the open paths remain for the two children, i...
X, L Ls, Rs. By faithfulness in (AI)
X, Zz, Rs. By (A4) X, X, Zg, Rş 0.
15

---

Published as a conference paper at ICLR 2024
A.3 PROOFS OF RESULTS IN $42
For results in $4.2 (including Proposition and Theorems 2 and (
regarding the discovery and
validation of the dropout mechanisms within a general causal dropout model without assumption [A3) please refer to Appendix C where we will first provide the necessary background and introduce the general framework to support our findings.
B DISCUSSIONS
B.1 EXISTING PARAMETRIC MODELS AS INSTANCES OF THE CAUSAL DROPOUT MODEL
In $2. we discussed several existing parametric models as specific instances of our proposed causal dropout model. Here we give some detailed analysis on them:
1. Dropout with the fixed rates (Example D.. D; ~ Bernoullipt).. i.e., the gene Zj gets dropped out with a fixed probability p: across all individual cells. The representative models in this category are the zero-inflated models (Pierson & Yau 2013; Kharchenko et al.
2014 [aeed et al. 2020; Yu et 2023 Min Agresti (2005) (or with slight difference, the hurdle models (Finak et al. 2015 Qiao et al. (023)), where sequenced data X; are assumed to follow a mixture distribution with two components: one point mass at zero for dropouts, and one common distribution, e.g., Gaussian, Poisson negative binomial, for gene counts. Biologically, dropouts with the fixed rates can be explained by the random sampling of transcripts during library preparation regardless of the true expressions of genes. Graphically, in this case the edge Z; + D; is absent. is worth noting that in some (e.g., the Michaelis-Menter (Andrews & Hemberg 2019) models, the dropout rate is determined by the average expression of the gene, i.e., Pi = f(E[Z))). In this case however, the edge Z; D; is still absent, as apparently Zi D; (this p, is fixed across all cells). 2. Truncating low expressions to zero (Example 2) D; 1(2 < <c). i.e., the gene Z; gets dropped out in cells whenever its expression is lower than threshold q. A typical kind of such truncation models is, for simple statistical properties, the truncated Gaussian copula (Fan et al.] 2017 Yoon et al.) 2020 Chung et al.] 2022). where Z are assumed to be joint Gaussian (usually with additional assumptions e.g., standardized), and the covariance among is estimated from X. Biologically, such truncation thresholds e; (quantile masking (Jiang et al. 2022)) can be explained by limited sequencing depths. Graphically, in this case the edge Z, + D; exists.
3. Dropout probabilistically determined by expressions (Example B
D: Bernoulli (F, + as)). where usually F; are strictly monotonically increasing functions in range [0. 1], and 3,, a, R are parameters with B; < 0, i.e., gene Z, may be detected (or not) in every cell, while the higher it is expressed in a cell, the less likely it will get dropped out. F, is typically chosen as CDF of some common distributions, e.g., probit or logistic (Cragg 1971 Liu 2004 Miao et al. 2016)). Biologically, the mechanism can be explained by inefficient amplification. Graphically, the edge Zi + D, also exists, and is non-deterministic.
Specifically regarding the several existing zero-inflated models (Pierson & Yau 2015 Kharchenko Yu et al. 2023 Min & Agresti 2005) (or related hurdle models (Finak et al. 2015 Qiao et al. 2023 that can be seen as parametric instances within our proposed causal dropout model, they align with Example I where each gene Z, experiences dropout with a fixed probability p; across all cells. However, it is worth noting that there are also some other zero-in inflatnde models and variations that fall outside the scope of our causal dropout model. We categorize these deviations based on the following three reasons:
al. 2014 Saeed al.
2020
1. Dropout is influenced by other covariates, rather than solely by the gene itself. In certain models, the excessive zero rate not fixed across all cells, meaning that Z, J. D,. Further- more, this dependency cannot be explained solely by Z;, as Z Z, 1. D An example of such a model one that models the zero rate using a logit link, as in (Workie & Azene 22021:: D; Bernoullip..) with p; = 3, Z. where 3, is the parameter vector in IR. Here, the non-zero entries of 33 are not limited to the i-th entry, and all other entries are also considered as "parents" of Z, in GRN (Choi al. 2020). This aligns with our findings in $4.2 and Appendix[C.. [C.2 where the causal effects within the GRN or as dropout mechanisms are unidentifiable.
16

---

Published as a conference paper at ICLR 2024
2. Measurement errors are involved. In our model, X, =1 (1- D,) Z₁, indicating that the latent variable Z is partially observed. with the non-zero observations being the true values. However, in representative models such the post Poisson model (Xiao et al. 2022: Saced et al. 2020), the data generating process is formulated as follows:
X (Possson(Zy (Z)
w.p. p:
(B.1)
w.p. 1 - P;
In this case, even if p; is fixed across all cells, it does not fit into our model as the non-zero parts of the observation X; do not reflect the true values of Z, due to the post-Poisson noise.
This scenario should be addressed separately as the problem of causal discovery in the
presence of measurement ervor (Fuller 2009
Pearl 2012 Kuroki & Pearl 2014 Scheines
& Ramsey 2016) Dai et al. (2022).
3. Excessive zeros are incorporated into the latent generation process, rather than being modeled as post-sequencing procedure. In our model, we have two sets of variables: Z, corresponding to the true underlying expressions, and X, corresponding to the observations. E.g., in (Xiao et al. 2022, Saeed al. (020),, each X, conditioned on Z; is assumed to follow a zero-inflated Poisson distribution. However, in certain models, there is no such differentiation. The observations X are assumed to be exactly the true expressions, meaning there are no technical dropouts. For example, in (Choi et al., 2020). the true expression of each gene, conditioned on its parents in the GRN, directly follows zero-inflated Poisson distribution, i.e., these zeros are inflated in the underlying genes interactions, not in post sequencing procedure. It is important address the dropout issue but rather to provide a more accurate characterization (beyond just Poisson) of the excessive zeros in the true gene expressions.
note that this type of model is not designed to
B.2 RELATION TO IMPUTATION METHODS AND THE MISSING DATA PROBLEM
In 22.4 we examined the potential theoretical issues associated with imputation methods and demonstrated the general unidentifiability of the true underlying joint distribution p(Z) from the observational distribution p(X). Now, we delve deeper into the derivation of this result and provide more comprehensive discussion of the broader missing data problem.
Firstly, according to the missing data litera- ture, the underlying joint distribution pZ)) is irrecoverable due to the self-masking
Row
Inguttd [a]
Imputed (b)
dropout mechanism (Enders] 2022 Little &
Rubin 2019 Mohan et al.) 2013] Shpitser 2016 Mohan 2018). To "fill in the holes" on any subset of variables S from non- zero observations, i.e.. to recover the joint distribution p (Zy) from
E214
E214
Ez14
Pans
imp ted
Impated [b]
0) (which is p( /Zs Rg = 0)), it generally re- quires (Zs) to be factorizable into compo- nents where the missingness is ignorable, ie.. p(-)= = p = O). unbiasedly es- timable from non-zero parts. By induc- tion, at least one gene's 's marginal distribu- tion p(Z,) should be recoverable, which is however impossible: since Z; directly af- fects the dropout by itself, Z, and D;, R& are dependent conditioning on any other variables, owing to the presence of the di- rected edges between Z; and D,, R& in Fig- ure I This irrecoverability persists even if the dropout mechanism is assumed to be parametrically fixed, as shown by the counterexample Example
45
Junb
6.0
Junb
Junb
Figure 7: Examples illustrating the introduction of spu- rious relationships through imputation. The two rows correspond to two gene pairs, namely (E2f4, Maff) and (Junb, Spil), obtained from unperturbed cells in the data by (Dixit et al. 2016). The three columns repre- sent the raw counts after normalization and the imputed data using MAGIC (Van Dijk et al. 2018) with (b) in- corporation of all genes or (c) considering only the 24 important genes discussed in $5.3 It can be observed that imputation introduces spurious relationships, such the presence of fork-shaped and highly nonlinear as- sociations. Furthermore, the disparity between (b) and (c) demonstrates that the selection of genes significantly impacts the outcomes of imputation.
Secondly, what distinguishes dropout data from missing data is that the missing en- tries cannot be precisely located, i.c., the
17

---

Published as a conference paper at ICLR 2024
D indicators are latent. In scRNA-seq data, we only observe zero entries but cannot differentiate between biologically unexpressed genes and technical zeros. Most existing imputation methods treat all zeros "missing holes" resulting from technical dropout and fill them in using information from non-zero entries. However, given that biologically unexpressed genes are ubiquitous in cells, such imputation can introduce false signals and spurious relationships (Andrews & Hemberg. 2018), as illustrated in Figure 7 Notably, many imputation methods only evaluate their accuracy under the sxeeddrrooout-atee scheme (Jiang et al., 2022 as in Example [ where dropout happens completely- at-random (CAR) (Little & Rubin 2019). I.e., Z L D. However, even in this relatively simpler case (without the self- masking as in the first reason), p(Z) remains irrecoverable, as shown by the counterexample Example
And thirdly, even if one is willing to make restrictive assumptions to render p(Z) theoretically recoverable, imputation methods maynot be the most suitable choice for gene regulatory network inference, which is inherently a structure learning task. This because imputation has a different inductive bias compared to structure learning: it focuses primarily on the unconditional (pairwise) dependencies and is better suited for tasks such as differential expression analysis. However, for structure learning, conditional dependencies are more important in order to capture the underlying gene regulatory relationships (Saced et al.] 2020] Gao et al.] 2022). ) As interestingly shown by (Gao) et al. (022), in the linear Gaussian setting where covariance matrix is sufficient statistic, imputation methods and a method specifically for causal discovery are compared on recovering the structure with missing data. Even if imputation methods estimate a more accurate covariance matrix (in terms of a smaller Frobenius norm), by inputting the estimated covariance matrices into causal discovery method e.g., PC, their performance on structure learning (in terms of SHD to the true graph) is worse than the other method specifically designed for structure learning. Therefore, imputation methods may not provide the optimal inductive bias required for accurate gene regulatory network inference.
In light of the demonstrated unidentifiability of p( ZZ) using the decomposition technique from the missing data literature, we now delve into the missing data problem and its relevance to our work. Missing data entries are frequently encountered in real-world data, such as unanswered questions in a questionnaire. These missingness patterns can be categorized into three main types: Missing Completely At Random (MCAR), Missing At Random (MAR), and Missing Not At Random (MNAR). Data are considered MCAR when the cause of missingness is purely random, such as instances where the missing rate is fixed across all samples (as in Example I and entries are deleted due random computer errors. On the other hand, data are deemed MAR when the missingness is independent with the corresponding underlying variables given all other fully observed variables (i.e., no missing entries at all), meaning that the missingness can be completely explained by the fully observed variables. For example, consider the gender wage gap study that measures two variables: gender and income, where gender is always observed and income has missing entries. In this scenario. MAR missingness would occur if men are more reluctant than women to disclose their income. Given the fully observed gender, the missing entries are independent with the income. Lastly, data that do not fall under either MCAR or MAR are classified as MNAR.
To handle missing data, trivial approach is list-wise deletion, where all samples with missing entries are discarded. Another common strategy is data imputation, often performed through expectation maximization (EM) techniques. In recent years, there has been growing interest in understanding missing data from a causal perspective, with contributions from various studics (Enders) 2022 Little| & Rubin 2019 Mohan et al. 2013 [Shpitser 2016). By modeling the missingness mechanisms within a causal graphical model, researchers have explored the conditions for recoverability of the underlying distribution and developed techniques for recovery. It has been shown that list-wise deletion can only recover the true distribution when the missing mechanism is MCAR. EM-based imputation methods can unbiasedly recover the true distribution only under the MAR missingness. For MNAR, some cases can be addressed through techniques like factorization or probability reweighting (Mohan et al. 2013). but certain situations are theoretically irrecoverable, such as the well-known self-masking case (Mohan) 2018). Some other studies propose that even if the true underlying distribution is irrecoverable, the conditional independence (CI) relations can still be identifiable through structure learning. This is achieved by employing test- wise deletion, where incomplete records of variables involved in each CI test are removed (Strobl et al., 2018 Tu et al., 2019).
In our specific task, we do not encounter missing entries explicitly marked as, for example, nan in the data matrix. Instead, we only observe zeros, which can be attributed to either technical or biological reasons that are unknown to If we treat all zeros as missing entries, self masking edges exist
18

---

Published as a conference paper at ICLR 2024
for each variable. Intriguingly, even in this challenging scenario, according to the current biological understanding to the dropout mechanism, our proposed causal dropout model can systematically address the dropout issue based on the current biological understanding of the dropout mechanism. As a result, the true graph structure can be consistently estimated. Lastly, it is worth noting that within our test-wise deletion only removes samples with zeros for the conditioned variables, rather than all variables. This approach deviates from the conventional notion of test-wise deletion, and we will further discuss this distinction in Appendix C.2
proposed framework, as illustrated in Definition [
C ON A GENERAL CAUSAL DROPOUT MODEL WITH RELAXED ASSUMPTIONS
Among the assumptions delineated in $4.1 the structural assumption ( (A3) concerning the dropout mechanism is of particular interest: each gene's 's dropout solely affected by the gene itself and not by any other genes. While this assumption is in line with most of the existing parametric models (2.3)
and plausible in the context of scRNA-seq data, where mRNA molecules are eyverse-ranssern independently across individual genes, it is still desirable to empirically validate this assumption based on the data. In other words, except for discovering the relations among genes, we also aim to discover the inherent mechanism to explain the dropout of each gene from the data. Surprisingly, as in $4.2 we propose a principled and simply approach to do so: the causal structure among Z. representing the GRN, as well the causal relationships from Z to R, representing the dropout mechanisms, can be identified up to their identifiability upper bound. Now here we will first give a detailed motivation and elaboration in Appendices C.1 and C.2 and then give the proofs in Appendix C.3
C.I ELABORATION ON THEOREM
3 WITH AN ILLUSTRATIVE EXAMPLE
Theorem 2 gives a principled approach to identify both the GRN and the dropout mechanisms up to their identifiability upper bound even without assumption Assumptions T Before we move into the proof in Appendix C.2 here we first explain what the identifiability upper bound is, and why.
To check whether gene Zi affects dropouts of another gene Zj. it basically requires to view R2 also as random variables, and to test for conditional independence and see whether there exists any vari- ables set that can d-separate Z, from Rg. However, there exists natural identifiability upper bound: Proposition 2 (Identifiability upper bound of dropout mechanisms). If 74 and Z,j are adjacent in GRN, then whether they affect each other's 's dropout, i.e., whether edges Zi + R2 and Zj + R₁ exist, unidentifiabl.
This is because condition Zj. Ry must also be conditioned but it is already in the estimands. Due to this, we only focus on identifying the non-adjacent pairs of variables in GRN, i.c., to recover the conditional independencies among Z. We have the following results:
Theorem 3 (Observed independencies are correct independencies). Assume (A) A2 (A4)
Vi,j.S.
X, L X X, 1 X, Zs. Rs = 0 2 7, Z, Zs, Rsu((,))
7, 1 Z, Zs
(C.1)
Theorem 3 provides insights that align with the previous Proposition 1 al- though the conditional independencies inferred from observational data with dropouts are "rare", they are accurate. The first three terms are corresponding to Proposition]] [ Definition [ and Definition 2 respectivvel, and the fourth is oracle. When the structural assumption (A3) is in the play, the converse directions of @ and ② also hold, as stated in Theorem [ However, in the general setting without [A3] these implications hold only in one direction. Progressing from left to right, the testing approaches provide increasingly tighter bounds the correct underlying conditional independence relations
Z
Z2
Z3
R3
R2
while also in empirical terms, they are less data-efficient due to the increased deletion of samples.
Figure 8: An illus- trative example for dropout mechanism identification.
Consider when we have identified a non-adjacent pair Z, and 2,. and the question arises: How can we determine the presence of causal relationships
Z; + R2 and Z, + R,? Should we conduct a conditional independence (CI) test between Z; and R, as suggested in Proposition ? However, engaging in such a procedure would lead astray. Let
2 As for how tight it theoretically can be, see Appendix C.2
19

---

Published as a conference paper at ICLR 2024
us examine Figure as example. We have identified Z and Za non-adjacent based on the conditional independence Z 1 Zz| R1,,,3 = 0. Then we test for causes of R3. However, even Z1 + R3 does not exist, we cannot establish their independence: Z, I 22, = 0 due to the presence of the collider R2. So, what should be our course of action? Surprisingly, the answer already lies before us: Z1 - R3 must not exist; otherwise, the independence Z1 1 Zz| Zzz R1,2,3 0 would not have been identified initially, owing to the presence of the collider Rg. In other words, the non-adjacency of 21 + R3 not established through a CI test, but rather through logical deduction.
The sparse causal graphs identified from real-world data also suggests the sparse dropout mechanisms, which, from the empirical view, supports our [A3)
Finally, we have come to the result in Theorem 2
C.2 MOTIVATION ON THE GENERAL CAUSAL DROPOUT MODEL
Generally speaking, we discover the causes of zero indicators R { by also viewing them as random variables and conducting causal structure search among them. As is mentioned in $2.2 the causal model in containss redundant deterministic edges and auxiliary variables just for convenience of notations. Now given the already built basic understanding, we proceed with an equivalent but more compact model, with only two parts of variables, Z and R.
Given that R; = [Xx; = 0) = D, OR 1(2 0). it follows that the edge Z; + R, always exists. Notably, even in the presence of dropout that happens completely at random, i.e., Z, Jl. D, as shown in Example| 1 the edge Zi - R; remains, due to the 1( = 0) term. In this section, we drop (A3) and seek to empirically verify it, while retaining the assumptions (A7) (A2) (A4) and [A5)
(specifically, we continue to assume (A2) i.e., dropouts do not affect genes underlying expressions).
To conduct causal structure search on variables Z U R. we need to identify d-separation patterns by performing CI tests. However, due to the definition of this causal model. some d-separation patterss naturally cannot be read off from the data. Specifically, the following rules should be followed:
Proposition 3 (Testability of CI queries), A CI query in form A LBBIC is testable if and only if:
1. (R, E (p) s.I. Z4 E (A,B) 2. A B and A. C,
BY C) C C,, and
where 1. means that any underlying variables involved the test is testable (accessible) only when the corresponding R variable is also conditioned, i.e., keep only the non-zero samples. 2., logether with I., rules out the queries e.g.. g. R, L2R, Zj or Zj L2 (i.e., R2 cannot be further conditioned).
Proposition 3 gives sufficient and necessary condition for a CI query to be testable from observational data X. Now we further give the criteria asymptotically estimate these CI queries from data X. For any CI query that is valid according to Proposition 3 it can be estimated by:
Proposition 4 (Testing CI queries from data). For each variable involved in a CI query A 13B C:
1. If it is a Zi variable, either in the bivariate estimands (A, B} or the conditioning set C, the non
non-zero zero samples of the corresponding observational values are used in the test.
2. If it is an R² variable,
a) If it is conditioned (i.e., R, E C) and the corresponding true variable is involved (i.e.,
2, E (A, B) U C), then it is conditioned by selecting non-zero sumples, ie., R; 0. b) Otherwise, all samples of the corresponding binary values of the variable R2 are used.
Specifically, regarding 2.b) of Proposition| 4 it coincides with a recent proposal to use only binarized counts, though with different purpose, cell clustering QQuu 2020). Note that in our CI estimation, R, can also be equivalently represented by full samples of the corresponding observation with zeros, X,, without binarization. This is because of the faithfulness in [AT] and the self-maste masking edges Zi + X,, which renders the identical structures, and thus d-separation patterns, among R and X. However, empirically, it is better to use X, as magnitudes of non-zero values may capture more fine- grained dependencies. E.g., Z1 ~ Unif[- æ]. Ze cos Z1, and D, 1 1(R, < 0). Then X1 L X2 (faithfulness holds), but this dependence in only captured by magnitudes of X2 when R3 = Ra = 0, not by the zero patterns represented binary indicators, i.e., R2 1L R2 (faithfulness violated).
With the essential infrastructure in place for reading off d-separations from data, we are now prepared to develop algorithms to recover the causal structure among Z U R. In addition to discovering the gene regulatory network within Z, our objective extends to uncovering the causes for each R,, i.., to discover the dropout mechanisms. Since R Dj OR 1 (Z+ 0), an edge 7, -> R₁ exists if and
20

---

Published as a conference paper at ICLR 2024
only if Z, - D.D..i.. i..., gene i's expression influences the dropout of gene j. However, we first notice that not all such dropout mechanisms are identifiable. See the natural upper bound in Proposition 2
Proposition 2 can be readily derived from Proposition to determine the existence of Zj + R. it is necessary to condition on at least Z, as otherwise the path 2, Z; - R; is left open. However, since R, is part of the bivariate estimands, it is impossible to condition on Z;, leading to a contradiction. Consequently, we can only identify edges of the form Z; - R, where Z, and Zj are non-adjacent, meaning there exists subset S [p] such that Z, 1 2z Zs.. The question then arises: since Z cannot be directed test on, how can we initially identify these non-adjacent pairs based only conditional independencies from the data with dropouts? Interestingly, we have the result in Theorem 3 observed independencies are correct independencies, though correct independencies may not always be observed.
Theorem 3 (Observed independencies are correct independencies). Assume (A7)
AA2
(A4)
Vi,j,S
X; T XX 2 x, 1 X, Zs, Rs = Zi 1 2, Zs, Rsu(()) = Z, 1 Z, Zs. (C.1)
Theorem 3 provides insights that align with the previous Proposition D although the conditional independencies inferred from observational data with dropouts are "rare", they are accurate. The leftmost hand side of Equation (C.1) represents the direct testing of conditional independence using complete observational data, as mentioned in Proposition[] The second term is precisely the approach proposed in Theorem T involving the removal of samples containing zeros in the conditioning set. The third term is more stringent, as requires the elimination of samples with zeros in all variables involved in the test, not just the conditioning set. When the structural assumption [A3]is is in the play, the converse directions of @ and @ also hold, as stated in Theorem ] However, in the general setting without (A3) these implications hold only in one direction: the observed independencies must be correct, but the correct independence may not be observed.
For instance, consider the graph Z, Za -> Z3 with [2, R.) and our focus is on 21 IL Za Z2. we introduce the additional edges Z1, 23 + R2, it serves as a counterexample to the inverse of 3 Similarly, if we add edges R2 Rl + R2g contradicts the inverse of @ Progressing from left to right, the testing approaches provide increasingly tighter bounds on the correct underlying conditional independence relations, while also in empirical terms, they are less data-efficient due to the increased deletion of samples.
Z
Zn
Z
> Z2
Z
> Z2
R22
R1
R22
R1
R2
(a)
(b)
(c)
Figure 9: Three exemplar causal dropout graphs that are unidentifiable with each other.
Having established that the influence of dropout mechanisms between adjacent gene pairs in a GRN is unidentifiable, i.e., the presence of an edge Z; + R, can only be identified when Z; and Zj are non-adjacent in the GRN, subsequent question arises: for non-adjacent pair of Z, and Z, but with Z: + R2 can such causal effect through dropout mechanisms be distinguished from the causal effect through gene interactions in GRN? The answer still generally negative. To illustrate this, let us consider the three examples depicted in Figure 9 Cases (b) and (c) are unidentifiable due to the result in Proposition [ while (b) and (a) are also unidentifiable, as all testable CI relations remain the same, e.g., R, 1 RR and Z, L Zz R, R2. This finding aligns with the conclusions presented in Theoremm3ans Appendix B.T] namely that the absence of a causal effect, whether mediated through dropout mechanisms or gene interactions, can be ascertained. However, determining the existence and precise nature of a causal effect generally remains elusive.
C.3 PROOFS OF THE RESULTS IN THE GENERAL CAUSAL DROPOUT MODEL IN $42
We first show the proof of Theorem B
i.e., observed independencies are correct independencies.
Proof Theorem 3 We prove for each implication as follows:
By faithfulness in (AI) X; T X Xs + R, L.R, Rs in the compact causal dropout model with only two parts of variables, Z and R,, as shown in Appendix C.2 By definition of d-separation,
21

---

Published as a conference paper at ICLR 2024
either there is no path connecting R, and Rj, or for any such path, there exists non-collider in Rs, or a collider s.t. it and all its descendants are not in Rg. Now add Zg into the conditioning set. For each path, collider could have itself or its descendants in Zs, as otherwise it is still already opened by the descendants Rg. Therefore, Rd la R2 [Zs, Rs. By Markov condition in [AI) X; 1 Xy|Zs, 0.
8: Show by its contrapositive. By Markov condition in (AI) 7, A. Zy|Zs, Rsu(k) Z: LiZ, Lz, Rs i.e., there must be path connecting 2, and Zj that is open conditioning on Zs. Rsu( 1) If none of such paths has any of R; or R on it, then consider one with two ends extended, Lee, R₁ < Z ..... Zj + R,: for any non-collider, including Z₁ and Zj. it is not in Zg, Rs; for any collider, either it or one of its descendants is in Zs. Rs and not in {R., R3 }. as otherwise there must exist another path with R. or R2 on it, e.g., Zi + R2 t.... collider already opened by R2 and rest variables, contradiction. 2) Otherwise, consider path with R2 or R2 on it. w.I.o.g.. consider the segmentation of this path with one end extended, i.c., R. Zy + Rg- Similarly, it is still open conditioning on Zs. Rs. Therefore, we have R₁ L Zs. Rs. By faithfulness (Af) and faithful observability (A4) we have X; L x Lzs Rs 0.
Z, which is
Also show by its contrapositive. By Markov condition in (AI) Z; 1 2, Ls + ZZ L Vzg By [A2) Rsu(k,) does not contain any non-collider on the open path connecting 7. and Zj. and thus involving them into the conditioning set will still keep the paths open, i.e., Z; , Rsu(t,). By faithfulness in (A7) and faithful observability in (A4) we have Z, 1 Zz Zs. Rsu(i.) = 0. 0
Then, to show the proof of Theorem we will present the sufficient and necessary condition for observing Z, IL Zig Rsu(i.) = 0, and examine the identifiability of the causal effects (either through gene interactions dropout mechanisms) under the condition.
(A4) and [A5) Theorem 2 (Identification of GRN and dropout mechanisms). Assume [AI] (22
In the CPDAG among Z output by Definition 2 if Zi and Zj are non-adjacent it implies that they indeed non-adjacent in the underlying GRN, and for the respective dropout mechanisms, 7e does not cause R, and Zj does not cause R₃, On the other hand, if Zí and Z are adjacent, then in the underlying GRN Zi and Z, are non-adjacent only in one particular case, and for dropout mechanisms, the existence of the causal relationships Zi + R2 and Zj + R& must be naturally unidentifiable.
Proof of Theorem 2 We have shown necessary condition for observing non-adjacent Z,, Zj as in Theorem 3 i.e., they are truly non-adjacent in GRN. Now we give the sufficient condition (stronger): for at least one Zs that d-separates Z, and Zj in the GRN, for any variable involved, i.e., Vk S U {i, j}. J. Rg is neither common children of Z, Zj. nor descendant of a common children of Z,, Z; (if any). In other words, at least one conditional independence between Z,, Z, needs to be observed.
We show this by contrapositive: if Zi, Z, are truly non-adjacent in GRN but for all Zs that d-separates Z, and Z, in the GRN, the corresponding conditional independencies cannot be observed, i.e. Z, L Zg but Z; Z Zs. Rsu(ig)) then for each of such S, there must be an open path connecting Z,, Zj in the form of Z, > W + Zy, where either W (s.j).)) or there exists a k S J {i. j} s.t. W is R2 ancestor. Otherwise, if the open path has more than one variable W. then at least one variable is not collider and can be conditioned on to block the path, contradiction.
Therefore, if Zi and Zj are non-adjacent in the output by Detinition 2, they must be truly non-adjacent and does not cause each other's 's dropout. Otherwise if Z, and Z; are adjacent in the output by Definition they are either truly adjacent, or in a non- adjaacettnt but particular case as shown above. In either case, the existences of 2, + R2 and Z, + RG are unidentifiable, as illustrated in Proposition 2 and Figure a)..
D SUPPLEMENTARY EXPERIMENTAL DETAILS AND RESULTS
D. LINEAR SEM SIMULATED DATA
Here we provide more detailed elaboration on the simulation study and implementation as in 5.1
22

---

Published as a conference paper at ICLR 2024
--
Ornele*
Testwise deletion
Full sumples
MAGIC
mixedCCA
PC
GES
Fibed prop 2006
20-
SHD CPDMG (loverts benee)
Tranntted
20
Logiss
Degna
Degeee
(a)
Z 10; Z
Gaussian.
(b)
Z 10; Z Lognormal.
PC
GBS
PC
GES
1.0
S01D of COPOAO (bower bener)
4-
190
Transttd
192
100
102
Logiete
52-
Dogeze
Degnne
(a) Z 20; Z
Gaussian.
(b) IZ
20; Z ~ Lognormal.
GBS
900
Findd Fi
592
5111 of EXDAA (loweer (otee)
202
Thuncatel
903
302
300
202
Engasti
Degree
Degree
(a) [ 100; Z Gaussian.
(b) 7Z
100; Z Lognormal.
Figure 10: Experimental results (SHDs of CPDAGs) of 10, 20, and 100 variables on simulated data.
Methods We investigate the effectiveness of our method (denoted as testwise deletion method) by comparing it to other methods on simulated data. Specifically. we consider the following baselines: (1) MAGIC (Van Dijk al. 2018) that imputes the dropouts, and (2) mixedCCA (Yoon al., (020) that applies mixed canonical correlation to estimate the correlation matrix of Z. The imputed data and estimated correlation matrix by MAGIC and mixedCCA, respectively, can used as input to suitable causal discovery algorithm. Furthermore, we also consider another baseline that directly applies a causal discovery algorithm to full samples (without deleting or processing the dropouts), which correspond to the method in Proposition [| We also report the results assuming that the complete data Z are available, denoted as Oracle*, in order to examine the performance gap between these methods and the best case (but rather impossible) scenarios without dropouts. We
23

---

Published as a conference paper at ICLR 2024
consider PC Spirtes et al. (0000 and GES (Chickering) 2002) the causal discovery algorithms for these methods. We describe the hyperparameters of these methods in detail in Appendix D.I
Simulation setup We randomly generate ground truth causal structures based on the Erdös Erdös-Rényi model (Erdós et al. 1960) with E {10, 20) nodes and average degree of to 6. On each setting five causal structures are sampled, corresponding to the confidence intervals in Figure 4 For each causal structure, we simulate 10000 samples for the variables Z according to linear SEM, where the nonzero entries of weighted adjacency matrix are sampled uniformly from (-1,- -.225) U (0.25, 1), and the additive noises follow Gaussian distributions with means and standard deviations sampled uniformly from (-1,1) and (1, 2), respectively. In this case, Z are jointly Gaussian; we also consider another case by applying an exponential transformation to each variable Zú, such that Z follow Lognormal distributions. Lastly, we apply three different types of dropout mechanisms on to simulate X: (1) dropout with the fixed rates, (2) truncating low expressions to zero, and (3) dropout probabilistically determined by expression, which are described in Examples D and respectively, For the first two categories of dropout mechanisms, we randomly set 30% to 70% of samples on each gene to get dropped out. For the third mechanism, we set D; ~ Bernoulli(lgiit, (-1.5Z; (logiti-.1ii11
0.5)).
PC We use the implementation from the causal- package Zheng et al.. 2023 For speed consideration, we use FisherZ. as the conditional independence test and set the significance level alpha to 0.05. Note that theoretically FisherZ cannot be used, as after the log transformation, dropout, or zero-deletion distortion, the data is generally not joint Gaussian. However, since the result already good, we did not further use the much slower but asymptotically consistent Kernel-based conditional independence test (Zhang et al., 2011). It is also interesting to see why FisherZ still empirically works here, and examine under which condition the FisherZ test is still asymptotically consistent even after test-wise deletion (one trivial case is Definition 2 where Z follows joint Gaussian and dropout happens with fixed rates). We leave this as an interesting future work.
We use the Python implementatior We modify the local score change by first calculating the test statisties of the corresponding testwise-deletion CI relation, and then rescale the partial correlation that it produces the same test statistics in the full sample case. This partial correlation and the full sample size are then used to calculate the BIC score change. The LO penalty is set to 1.
GES
It is observed that our proposed method (1) has a much lower SHD compared to the baselines across most settings, especially when PC is used as the causal discovery algorithm, and (2) leads to SHD close to the oracle performance. Furthermore, one also observes that directly applying causal discovery to full samples results in a poor performance in most cases, possibly because, as indicated by Proposition I the estimated structures by this method may contain many false discoveries, leading to a low precision. worth noting that mixedCCA is specifically developed for jointly Gaussian variables and truncated dropout mechanism. Therefore, compared to the other settings, it performs relatively well in this setting, despite still having a higher SHD than our proposed method. This further validates the effectiveness of our method that outperforms existing parametric models with no model misspecification. One may wonder why test-wise deletion method outperforms Oracle* in the Lognormal setting. A possible reason is that, after zero deletion, many small values (especially for truncated and logistic dropout mechanisms) are removed; therefore, the resulting distribution is closer to Gaussian, which is more suitable for the specific CI test adopted for these methods.
Experimental results for Figure 4
The results of 30 nodes are provided in Figure 4
On different scales As a complement to Figure 4 (30 nodes), here we also present the results of SHDs of CPDAGs in different scales (small as 10 and 20 nodes, and large as 100 nodes), as shown in Figure 10 Impressively, our proposed test- wise zero deletion still performs best among baselines across various settings.
Additional result 1:
Moreover, to better echo the high dimensionality in real-world scRNA-seq data (e.g., around 20k genes for humans), also conduct experiments in super large scale: The ground-truth graph has 20,000 nodes and an average degree of 5. The sample size is 20,000. The dropout rate is 70%. We used the FGES (Ramsey et al. 2017 implementation on an Apple MI Max with 10 cores. Here are the results:
'https://github... http:/gghhhhucoojjuuaaammmllalaa http:://gituu..oo//muuuu
24

---

Published as a conference paper at ICLR 2024
Omnele
Testwise delction
Full samples
MAGIIC
mixedCCA
PC
GES
PC
GES
1-24
FI-06
RI-BS
1.0
Fixed dropout raile
0.8
0.8
G B 99
0.5
0.-02
0.2
Skeleton precision
1.0
0.8
Truncuted
0.5
0.2
0.0
0.0
1.0
Logistic
0.5
0.2
0.2
0.25
0.6
0.25
0.6
0,25
0.6
0.25
0.6
Skeleton recall
Skeleton recall
(a) [Z|
30; Gaussian.
(b) |2| 30; Lognormal.
Figure 11: Experimental results of 30 variables simulated data (Figure 4 using another metric: precisions and recalls of skeleton edges. Each subplot corresponds to specific algorithm (PC or GES) specific base data distribution (Gaussian Lognormal) with specific dropout mechanism (Fixed dropout rate, Truncated, or Logistic). In each subplot, the 5 lines correspond to the 5 dropout- handling strategies (Oracle* Testwise deletion, Full samples, MAGIC, mixedCCA). On each line, the 6 markers with sizes from small to big correspond to the average graph degrees from to 6.
Table 1: FGES results on a large-scale graph with 20,000 nodes.
Oracle*
Odel
Full samples*
Skeleton Fl score
0.93
0.90
0.14
Time
1421s
1512s
> 24 hrs
As FGES directly on samples with dropouts has not ended within 24 hrs (since there too many false dependencies), we reported the Fl-score (0. 14) given by the graph at the last iteration before we terminate it. We see that: 1) Zero-deletion is effective in recovering the true structures, and 2) Zero-deletion is fast, because it corrects the false dependencies and the graph in search is sparser than that on dropout-come contamintcd data.
Regarding the speed (particularly in the high dimensionality). it worth noting that our proposed zero- deletion method itself is straightforward technique that introduces additional time complexity. As simple and versatile tool, it can seamlessly integrated into various other algorithms, and the whole time complexity is completely contingent on the specific algorithm accompanies (e.g., PC, GES, or FGES, or GRNI-specific methods in $5.2).
Additional result 2: On different metries than SHD As SHD metrics can sometimes be mislead- ing, we also calculate another metric, precision and recall of skeleton edges, the simulation setting
25

---

Published as a conference paper at ICLR 2024
Dropout rate
10
50
90
Method
Full
Full
Odel
Full
0del
Full
0del
Full
0del
Full
Odel
Skel Precis
.98
.58
.98
.54
.98
.57
.98
.54
1.0
.44
1.0
Skel Recall
1.0
.98
.98
.96
.89
.98
.89
.8
.60
.82
.16
Skel
.99
.73
.98
.69
.93
.72
.93
.70
.75
.57
27
PDAG SHD
19
56
23
62
28
55
28
61
28
81
43
(a)
100000
Drspont
F1a0.6
F1
0.8
Fi-04
Sample size with Rs = 0
31%
CI tests precision
0.6
0.4
0000
0.2
0.2
0.4
0.6
0.8
1.0
Conditioning set size [8]
CI tests recall
(b)
(c)
Experimental results to test the validity of assumption (A5) (effective sample size after test-wise deletion) with varying dropouts. On a dataset simulated with 30 variables and average graph degree of 3, truncated dropout mechanism is applied to produce 5 datasets with dropout rates 10%0 30%, 50%, 70%, and 90%, respectively, shown by the colors above. PC with test-wise deletion is run each dataset. (a): The structure identification accuracies w/wo test-wise deletion; (b): The remaining sample size after test-wise zero deletion, for each CI test during the PC run, with conditioning set sizes [S] growing from 0 to maximum of 4. The marker and the band show the mean and the standard deviation of sample sizes under a specific size [S] (c): The precisions and recalls of all CI tests (d-separation conditional independencies the PC run, where each colored line is on the data with a specific dropout rate, and each marker on the line is with specific conditioning size [S). Marker sizes from small to big correspond the conditioning size [S] from 0 to the bigger.
Figure 12:
counted positives) during
of Figure (because PC and GES output edges without strengths, we can only calculate the precision and recall without thresholding curves). The results are shown in Figure II
The results strongly validate our method's efficiency: we see from Figure that the orange lines (Testwise deletion) are often above all other lines, i.e., achieving consistently highest precisions. This aligns with our Proposition I and Theorem [ testwise deletion can correctly recover conditional independencies, i.c., reduce false dependencies, and thus reduce fälse positive estimated edges.
Additional result 3: On varying dropout rates While in $4 we demonstrate the validity of assumption [A5] (effective sample size after testwise deletion) in real-world examplar dataset (Figure 3) here we also conduct set of synthetic experiments with varying dropout rates to test the validity of assumption (A5) On dataset simulated with 30 variables and average graph degree of 3, examined the truncated dropout mechanism with dropout rates varying from 10% to 90%.
Figure 12 depicts the results of the remaining sample sizes in PC runs, the accuracies of CI tests, and the overall structure identification accuracies. These results support the validity of assumption [AS] even when the dropout rate as high as 70%, the remaining samples sizes for CI tests are still considerably high to avoid Type-Il errors, leading consistently better SHDs and skeleton F1-scores.
26

---

Published as a conference paper at ICLR 2024
D.2 REALISTIC SYNTHETIC AND CURATED DATA UNDER THE BEELINE FRAMEWORK
Except for the linear SEM experiments, we also conduct synthetic experiments that are more 'realistic' biological setting (closer to real scRNA-seg data) and more 'competitive' (on algorithms specifically desgined for GRNI, not just PC and GES). Specifically, we use the BEELINE GRNI benchmark framework (Pratapa et al. 2020). For data simulation, we use the BoolODE simulator that basically simulates gene expressions with pseudotime indices from Boolean regulatory models. We simulate all the 6 synthetic and 4 literature-curated datasets in (Pratapa et al. 2020), each with 5,000 cells and a 50% dropout rate, following all the default PC and GES, we have attempted to encompass all 12 algorithms benchmarked in (Pratapa al.) 2020). However, due to technical issues, five of them were not executable, so we focused on the other seven ones: SINCERITIES, SCRIBE, PPCOR, PIDC, LEAP, GRNBOOST2, and SCODE, For each algorithm, five different dropouu-handding strategies are assessed, namely, oracle*, testwise deletion, full samples, imputed, and binarization QQuu 2020
For algorithms, in addition to
Our test-wise zero-deletion strategy, as a simple and versatile tool, is employed in these SOTA algorithms as follows. We divided the seven SOTA algorithms into three categories: 1) For algorithms that explicitly compute (time-lagged) partial correlations, including SCRIBE, PPCOR, and PIDC, the test-wise cero-deletion can be readily plugged in as in PC and GES; 2) For algorithms that entail gene subset analysis (e.g., regression among gene subsets) but without explicit partial information, including SINCERITIES and I.EA,, in each subset analysis we delete samples containing zeros in corresponding genes; and 3) For algorithms that estimate the structure globally without any explicit subset components, including GRNBOOST2 and SCODE, we directly input data with list-wise deletion, i.c., delete samples containing zeros in any gene.
Note that the evaluation setting is not in favor of our implementation (especially for PC and GES with testwise deletion), due to the following three reasons: 1) Model mis-specification for PC and GES, both of which assume i.i.d. samples while BoolODE simulated cells are heterogeneous with pseudo- time ordering. These pseudotime indices are input the SOTA algorithms additional information, when applicable; 2) Model mis-specification for the CT tests. For speed consideration, in PC and GES we simply use FisherZ and BIC score which assume joint Gaussianity, instead of non- non-parametric CI tests and generalized scores; and 3) Other 7 SOTA algorithms output edges with strengths so we report the best Fi-score through thresholding, while PC and GES do not output strengths.
Even though, the results still strongly affirm the efficacy of our method's 's efficacy, as supported by the Fl-scores of the estimated skeleton edges in Figure See analysis in $5.2 Note that in Figure |5| the minimum of the colorbar is anchored to 0.5 in order to visualize the differences among the majority of entries with a finer resolution among all the 450 entries there are 380 ones with values 2 0.5. F1 values <0.5 yields a relatively poor performance and lacks significant reference value.
Notably, one may wonder why there is a seemingly "performance degeneration" " from full data to zero deletion for PC algorithms in the four curated BEELINE datasets. In fact, this performance degeneration is more likely attributed to the PC algorithm itself, rather than the zero-deletion method. Among the 4 curated datasets, there is performance degeneration on two of them: HSC and GSD. Notably however, it is exactly on these two datasets where tull-samples performs even better than oracle*, i.e., dropouts even do "good", instead of harm. This contradiction with theory makes the further analysis of zero-deletion less meaningful on these 2 cases. While on the remaining 2 cases (F1 ooaclee F1-full samples), exactly we also have 1-Odde==F1-full samples. We further extend this insight across all cases:
Among the 65 dataset-algorithn pairs where dropouts indeed do harm (F1-full sample<<F- oracle*), zero-deletion yields efficacy in dealing with dropouts 111-del>-F1-full sample) on 47 (86.2%) of them.
Among the 17 alloshtta pairs where dropouts unreasonably do good (F1-full sample>/F-oorccles zero-deletion yields efficacy (F1-Odel> sample) on only 4 (23.5%) of them.
Interestingly, the above statistics reaffirms the efficacy of zero-deletion (Theorem I]: zero-deletion yields CI estimations (and thus GRNI results) more accurate/faithful to oracle* regardless of whether the oracle* is better or worse than full samples. As for the reasons why in some cases dropouts unreasonably performs even better, will investigate more into details later.
27

---

Published as a conference paper at ICLR 2024
D.3 WORLD EXPERIMENTAL DATA
We conduct experiments on three representative real-world single-cell gene expression datasets to verify the effectiveness of our proposed method, including data from Bone Marrow-deriver Dendritic Cells (BMDC) (Dixit et al. 2016), Human Embryonic Stem Cells (HESC) (Chu et al. 2016) and Chronic Myeloid Leukemia Cells (CMLC) (Adamson et al. (016). The three datasets consist of 9843, 758, and 24263 observational samples respectively. Following the previous exploration of genetic relations, we determine the "ground-truth" GRN using the known human gene interactions from the TRRUST database (Han et al. 2018).
For the BMDC (Dixit et al. 2016) data, we follow the practice in (Dixit al. 2016 Saced et al. 2020) and only consider the subgraph on 21 transcription factors believed to be driving differentiation in this tissue type. As is shown in Figure using the deletion based CI test in Definition | the PC estimated graph is much sparser, i.e., more true conditional independencies are identified, as shown by the red edges. Furthermore, more known interactions identified, shown by the blue edges. Note that for PC's output, for several undirected edges in CPDAG, we manually orient them in the favored direction from known interactions.
Table 2: SHD Results of PC on HESC and CMLC datasets, using different methods as in Figure
Subgraph
Testwise deletion
Full
samp camples
MAGIC
HESC
12 nodes, edges
16
20
23
16 nodes, 13 edges
20
26
23
16 nodes, 11 edges
19
30
25
18 nodes, 17 edges 19 nodes, 16 edges
26
38
29
25
34
28
CLMC
16 nodes, 14 edges 17 nodes, 15 edges
28
49
41
33
65
46
18 nodes, 15 edges 21 nodes, 16 edges
32
70
47
41
98
42
23 nodes, 18 edges
33
106
53
For the rest two datasets HESC (Chu et al., 2016) and CMLC (Adamson et al. 2016), to enable more comprehensive experimental comparison while also considering the computation efficiency of causal discovery algorithms, we randomly we randomly sample five subgraphs on each dataset to evaluate the performances of different methods. The subgraphs criteria: no loops, no hidden confounders, and each gene should have a zero sample size less than 500 for HESC, or zero-rate less than 70% for CMLCC The last condition is for the test power of deletion based CI tests: if gene have too many zeros across all cells, then for any CI relation conditioned on it, the sample size after removing the zero samples will be too small, leading an empirical bias. Thus we do not consider those genes. This is reasonable as if a gene is almost always zero across all cells (either not expressed all dropped out), may not be that important or actively involved in the genes interaction in these cells. On these sub-datasets, we conduct PC on full observational samples, on test-wise deleted samples, and on MAGIC imputed data, respectively. The SHD results are shown in Table We notice that the graphs produced by Testwise tese deletion have the smallest SHDs, i.e. closest to the groundtruths. One may further be curious why the SHDs are generally larger than the number of edges in ground-truth i.c., even empty graph would have smaller SHD. This is because the graph discovered from data is usually much denser than the groundtruths, i.c., there may be additional gene interactions that yet not known not recorded in the TRRUST database (Han) jet al. 2018). To this end, we also examine the number of true positive edges (i.e., the detected true interactions) returned by each method, where restwise deletion is slightly better than Full samples, and significantly better than MAGIC.
sampled with the following
28


Reviews:
[{"review_id": "1Vu5GI8JIJ", "reviewer": "Reviewer_f2vv", "summary": "The authors proposed a causal graphical model, named causal dropout model, to characterize the dropout mechanism in scRNA-seq data. \nThey found that simply ignore the data points in which the conditioned variables have zero values can still lead to consistent estimation of conditional independence (CI) relations with those in the original data.", "questions": "For network inference, they can use some evaluation metrics such as ROC curve or PR curve to assess how well their predicted network recovers the true network. \nThey should conduct more experiments to show the performance gained by using their causal dropout model.\nSeveral gene network inference methods have been designed to handle missing values in scRNA-seq data. Therefore, as a practical analytical framework, the authors should prioritize the comparison of their model with the most advanced existing network inference methods.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The task of inferring gene regulatory network is of interest especially in the bioinformatic domain.\nThe empirical results seem to indicate that the approach can be integrated into existing causal discovery methods to handle dropouts.\nWriting and presentation skill is well.", "weaknesses": "For network inference, they can use some evaluation metrics such as ROC curve or PR curve to assess how well their predicted network recovers the true network. \nThey should conduct more experiments to show the performance gained by using their causal dropout model.\nSeveral gene network inference methods have been designed to handle missing values in scRNA-seq data. Therefore, as a practical analytical framework, the authors should prioritize the comparison of their model with the most advanced existing network inference methods.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}, {"review_id": "SaIrMNNskB", "reviewer": "Reviewer_pofH", "summary": "The authors propose a method for causal discovery of gene regulatory networks (GRNs), i.e., causal gene-gene relationships. Instead of conducting conditional-independence tests on the raw data without consideration of the dropout patterns of single cell sequencing (e.g., doing PC algorithm on the entire data), they propose conducting the tests only on non-zero conditioning variables, by showing that the conditional independence relations of variables conditioned on variables with non-zero values are the same as for data without dropouts. The authors demonstrate that their method outperforms competing state-of-the-art methods in causal discovery of GRNs on several data sets.", "questions": "- The authors state \"while a zero entry of $X_2$ may be noisy (i.e., may be technical), a non-zero entry of $X_2$ must be accurate, i.e., biological.\". As far as I can tell, non-zero values might also be technical due to sequencing/mapping/algorithmic errors, correct? \n- Supposing elevated dropout rates of, say 50%, the method relies on the fact that a conditioning variable can be found which \"breaks\" dependencies. Is this correct?\n- Using Fisher's $z$-test assumes multivariate Gaussianity. Wouldn't a kernel-based independence test be better for the log-normal data?\n- Figure 5 is not very readable. I believe a simple table or something similar could improve the presentation.\n- The source code could be improved and properly documented (What libraries are required? Which versions of the libraries did you use for validation? How do I run all experiments? Etc etc.).", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The method is in my opinion very original and will certainly be of actual use in the field of computational biology.\n- The paper is clearly written and easy to follow. As far as I can judge, the authors demonstrate an excellent grasp of the contemporary literature, both in causal discovery as well as in computational biology.\n- The method outperforms state-of-the-art methods for GRN inference in several benchmarks.\n- The experimental section is convincing, and should be easy to reproduce.", "weaknesses": "I do not have major comments on possible weaknesses.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}, {"review_id": "LiwEGHnsaC", "reviewer": "Reviewer_6osj", "summary": "The authors tackle the task of inferring gene regulatory networks (GRNs) from single-cell RNA sequencing (*scRNA-seq*) data, which is made difficult by the fact that *scRNA-seq* data showcases many zero values, which are either due to technical reasons (dropouts) or biological (no gene expression). This missing data problem can lead to overly-dense graphs being produced by SOTA causal discovery algorithms applied to this task. The authors propose a `Causal Dropout Model` to characterize the dropout mechanism and come up with a simple solution to handle dropouts, namely conditioning on non-zero entries in the conditioning set when performing conditional independence sets. They show that this approach is sound under relatively mild assumptions and performs well on a large number of synthetic, semi-synthetic, and real-world data sets.", "questions": "1. On page 6, when discussing the connection between BIC score and Fisher-Z test statistics, how do you conclude that \"only local consistency is actually needed\"? The reference in Nandy et al. (2018) shows this connection, but does not say anything about the assumptions made in the GES paper. Could you elaborate on this point?\n2. On page 7, what do you mean by \"identifiability upper bound\" for identifying the dropout mechanisms? How is this upper bound characterized?\n3. Perhaps I missed it, but in Theorem 2, what is the \"one particular case\" in which Z_i and Z_j are non-adjacent in the underlying GRN.\n4. In Figure 4(b), how do you explain that testwise deletion performs better than the oracle? Shouldn't that be the best case scenario?\n5. I haven't seen *causal sufficiency* mentioned anywhere, which is an important assumption made when using algorithms like PC and GES. Is it reasonable for this type of data to assume that there are no hidden confounders? In either case, I would say something about this important practical limitation.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The problem tackled is very important, as gene regulatory networks can provide direct insight into the workings of biological mechanisms, and *scRNA-seq* is increasingly available for this task. The authors present their idea in an almost flawless manner, with sufficient attention being given to describing related work, providing relevant examples, and to framing and testing the assumptions made for the `Causal Dropout Model`. Finally, the authors showcase the performance of their approach in an extensive series of experiments, in which they examine different dropout mechanisms, different causal discovery and GRN inference-specific algorithms, on multiple settings and types of data.", "weaknesses": "My only (minor) gripes are that most of the theoretical analysis is deferred to the appendix, which can sometimes lead to questions due to insufficient detail (see below), and that the references are a bit sloppy. \n\nMiscellaneous comments:\n- page 4, Example 4: I think the Bernoulli distributions should be reversed, since there should be a minus sign in the denominator exponential when computing the logistic function.\n- page 10, reference typo: * after \"Tabula Sapiens Consortium\"\n- page 11, formatting error: \"ALBERTS\" should not be capitalized\n- page 11, there seems to be no reference to Gao et al. (2022) in the paper.\n- page 12, duplicate author: \"Paul R. Rosenbaum\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (oral)