Causal discovery with endogenous context variables

Authors: wiebke  günther
oana-iuliana  popescu
martin  rabel
urmi  ninad
andreas  gerhardus
jakob  runge
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 15:49:14

Abstract: Systems with variations of the underlying generating mechanism between different contexts, i.e., different environments or internal states  in which the system operates, are common in the real world, such as soil moisture regimes in Earth science. Besides understanding the shared properties of the system, in practice, the question of context-specific properties, i.e., the change in causal relationships between contexts, arises. For real-world data, contexts are often driven by system variables, e.g., precipitation highly influences soil moisture. Nevertheless, this setup needs to be studied more. To account for such endogenous contexts in causal discovery, our work proposes a constraint-based method that can efficiently discover context-specific causal graphs using an adaptive testing approach. Our approach tests conditional independence on the pooled datasets to infer the dependence between system variables, including the context, to avoid introducing selection bias. To yield context-specific insights, conditional independence is tested on context-specific data. We work out the theoretical framework for this adaptive testing approach and give a detailed discussion of the connection to structural causal models, including sufficiency assumptions, which allow to prove the soundness of our algorithm and to interpret the results causally. A simulation study to evaluate numerical properties shows that our approach  behaves as expected, but also leads to a further understanding of current limitations and viable extensions.

Main paper: 
 introduction: 
 In recent years, causal discovery (CD), that is, learning cause-effect relationships from observational data, has garnered significant interest across various domains, from Earth sciences and genomics to industrial settings (11; 26; 27; 1) . Such complex systems often exhibit a change of causal mechanisms over time or across different environments, including changes in the associated causal graph. An example are the causal relations between surface heat and latent heat and their dependence on soil moisture regimes: When the soil is dry, it triggers a dependence of soil moisture on latent heat due to the dry soil reflecting heat, while environments with wet soil do not feature such a feedback. These changing causal mechanisms can be represented using so-called contexts, where the causal mechanisms and structures vary depending on the states of one or more context variables that augment the assumed underlying structural causal model (7; 17; 15) . Both time-series (29) and non-time-series data can exhibit such behavior. Context-dependent systems have been widely studied; see Mooij et al. (15) for a review. Context variables (15; 34; 18) have also been called selection variables (17) or regime indicators (23; 29) in the literature.Causal discovery on such systems is challenging (7; 15) and methods for CD in case of systems with changing causal mechanisms typically assume that the context variable is exogenous, as will be discussed below and in Sec. 3.2. Yet there exist many real-world examples where the context indicators are endogenous. A regime/context variable such as soil moisture can be influenced by the dynamical system variables -precipitation can make a dry soil wet -leading to joint systems with endogenous context variables, as illustrated in Fig. 1 , where we show a hypothesized causal model of the drivers of soil moisture (SM). Understanding context-specific causal relationships is crucial for applications such as extreme event forecasting and for enhancing our general understanding of dynamical systems.Figure 1 : In this strongly simplified example, the variable soil moisture (SM) is an endogenous context variable influenced by variables such as precipitation (TP) and latent heat flux (LH) in both wet and dry regimes of SM, represented by solid edges. In dry conditions only, a feedback loop from the latent heat flux to the soil moisture is created because the soil additionally reflects the heat, represented by the dashed edge (a feedback loop is a time-delayed cycle, leading to a cyclic summary graph representation while the underlying time-series graph remains acyclic).Our work investigates the assumptions for constraint-based CD on such systems. A core difficulty in treating endogenous context indicators is that context-specific CD methods which resort to masking, i. e., conditioning on the context of interest and ignoring data from other contexts, e.g., (7; 29) , can introduce selection bias if the assumption of exogenous context variables is violated. For example, applying CD only on samples from the dry context in the soil moisture problem of Fig. 1 would introduce a spurious correlation between the precipitation and latent heat flux. Other methods apply CD on the pooled data, i.e., the concatenated dataset containing all contexts. Then allowing for endogenous variables is possible (15) , but does not yield context-specific information, i.e., no distinction between solid and dashed edges in Fig. 1 can be made. We propose an adaptive constraintbased discovery algorithm that allows to distinguish context-specific causal graphs and give a detailed discussion of the connection to structural causal models, including sufficiency assumptions, which allow to prove the soundness of our algorithm and to interpret the results causally.Context-specific graphs contain more information than union graphs on pooled data, i.e., the dashed edges in in Fig. 1 . This information can under suitable assumptions be discovered by additionally taking context-specific independencies (CSIs) into account. CSIs are independence relations that are only present when data of a single context is viewed in isolation. Naively searching through all CSIs would lead to an inflated search space of available independence tests. Further CSIs necessarily use only a reduced sample size which can also increase error rates. Indeed, a simple baseline method we compare against suffers from an increased error rate, whereas our proposed method is able to avoid an inflation of required tests. This is achieved by an easy-to-implement decision rule which allows to adaptively decide for one and only one test to be executed per separating set. Whenever consistent, tests are further executed on the pooled dataset, thus on all available data.We now summarize our contributions: (1) We propose a general framework for adapting existing constraint-based CD algorithms to extract context-specific information, in addition to joint information, about skeleton graphs when contexts are (potentially) endogenous without inflating the number of required independence tests. (2) We give an SCM-based interpretation of this additional information and unveil interesting and little-studied difficulties whose understanding we consider essential to the problem. (3) We work out the theoretical assumptions under which our proposed method soundly recovers the context-specific skeletons despite these difficulties. (4) We conduct a simulation study to numerically evaluate the performance of the proposed method and understand the failure modes of current approaches. 
related work: 
 The problem of CD from heterogeneous data sources has gained interest in recent years. Mooij et al. (14) offers a broad overview of the topic and groups recent works on CD for heterogeneous data into two categories: Works that do CD separately on each context-specific dataset (the data points belonging to a specific context) and works that use the pooled datasets (all data points). Works that obtain context-specific graphs apply CD on the context-specific dataset, where the assignment to contexts is either known (8; 24) , identified before applying CD (12; 36) , or iteratively learned (29) .The information from the context-specific graphs can then be summarized in a single graph, the so-called union graph. Mooij et al. (14) propose the Joint Causal Inference (JCI) framework to deal with multi-context CD by modeling the context variable and using CD on the pooled datasets to uncover the relationships between the context and system variables, where the context can also be endogenous. Several CD methods can be extended for JCI, e.g., FCI (32) . FCI is also applied in (28; 10) on the pooled data to discover causal relationships among all variables. Zhang et al. (39) test for conditional independence between system and context variables on the pooled dataset under the assumption of exogenous context variables. Constraint-based CD for the multi-context case has also been studied specifically for the time-series case, e.g., (29; 33; 5) . Among non-constraint-based methods, Peters et al. (19) ; Heinze-Deml et al. (6) exploit the invariance of the conditional distribution of a target variable given its direct causes across multiple contexts, while Zhou et al. (41) find causal structure by modeling the direct causal effects as functions of the context variables. Although not directly methodically related, CD for cyclical graphs is also relevant to our problem, as cycles can appear between the different contexts (8; 3) . Most above-mentioned methods using both context-specific and pooled data assume, either implicitly or explicitly, that the context variable is exogenous. This is possibly due to the fact that, for context-specific approaches, endogenous variables might lead to selection bias, while for methods that use pooled data it suffices to treat endogenous variables as any other variable in the graph. In our work, we discover context-specific causal graphs for the case of endogenous context variables. Assuming access to the context variable, we obtain information from the pooled and context-specific graphs. We focus on the non-time-series case and leave the extension to time-series for future work.Labeled Directed Acyclic Graphs (LDAGs) (18) are conceptually closest to our work. LDAGs summarize context-specific information in graphical form by assigning a label to each edge according to the context in which dependence arises. Hyttinen et al. (9) propose a constraint-based method to discover LDAGs from data similar to our approach: A variant of the PC algorithm that tests conditional independence of X ⊥ ⊥ Y |R = r within each context R = r instead of testing X ⊥ ⊥ Y |R on the pooled data. We discuss how our definitions of context-specific graphs fundamentally differ from CSIs in Sec. 3.1 and delimit our work from LDAGs in App. A.Question: Does the paper discuss the limitations of the work performed by the authors?Answer: [Yes]Justification: We explicitly state the assumptions of our method, we evaluate our method on a wider set of experimental setups to investigate its robustness, and explicitly discuss the limitations of our work in Section Sec. 6. We discuss computational complexity in Sec. D.9.Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 
formal description: 
 We first briefly outline preliminary definitions (see, e.g., (16) for detailed definitions). We then point out an unexpected difficulty of our approach in the connection to SCMs and its remedy.Preliminaries We work within the framework of structural causal models (SCM) (16) . We define the set of observed endogenous variables {X i } i∈I , which take values in X i , with finite index set I. Let V := {X i | i ∈ I} represent the set of all endogenous variables. The set of (hidden) exogenous noises {η i } i∈I , taking values in N i is denoted U := {η i | i ∈ I}. For any subset A ⊆ V , let X A := j∈A X j . We denote X := X V and N := N U .A set of structural equations F := {f i |i ∈ I} is a set of parent sets {Pa(i) ⊆ V |i ∈ I} together with functions f i : X Pa(i) ×N i → X i . An intervention F do(A=g) on a subset A ⊆ V is a replacement of f j → g j for j ∈ A. We will only need hard interventions, g j ≡ x j = const below.Given a distribution P η of the noises U = {η i } i∈I , a collection of random variables V = {X i } i∈I satisfies the structural equations F on U if X i = f i ((X j ) j∈Pa(i) , η i ). An SCM M is a triple M = (V, U, F), where V solves the structural equations F given U . The intervened model M do(A=g) is the1 SCM M do(A=g) := (V ′ , U, F do(A=g) ). A causal graph describes qualitative relations between variables; given a collection of parent-sets {Pa i } i∈I , on the nodes I add a directed edge p → i for all i ∈ I and all p ∈ Pa i . See next section Sec. 3.1.We observe data from a system with multiple, slightly different SCMs, corresponding to contexts and indicated by a categorical context indicator / variable, as illustrated in example 3.1. This variable R ∈ V is a variable whose value indicates the context, i.e., the dataset to which the data point belongs. Therefore, the value-space of R defines all datasets, and its values associate data points to datasets. Example 3.1. Given a binary context indicator variable R and a multivariate mechanism of the formY := f Y (X, R, η Y ) = 1(R)g(X) + η Y , the dependence X → Y is present in the context R = 1, but absent for R = 0. This entails a context-specific independence (CSI) X ⊥ ⊥ Y |R = 0 and X ̸⊥ ⊥ Y |R = 1.In order to create a connection between CSIs and SCMs, we define multiple graphical objects. We motivate their necessity and exemplify them in Sec. 3.2.Given a set of structural equations, F, we define the mechanism graph G[F], constructed via parent sets Pa i specified via F (see "causal graphs" above). Given a pair consisting of a set of structural equations F, and the (support of) a distribution of the observables P (V ) we define the observable graph G[F, P ], constructed via parent sets Pa ′ i ⊆ Pa i (of X i ), such that j / ∈ Pa ′ i if and only if for all values pa of Pa i -{j} the mapping x j → f i | supp P (Pai) (x j , pa, η i ) is constant (where defined).M = (V M , U M , F M ), where V M is distributed according to P M , is the observable graph G union [M ] := G[F M , P M ]which consolidates the information from multiple contexts into a unified graphical representation.If F = F min satisfies standard minimality assumptions (2, Def. 2.6), then G[F min ] = G union [M ](given faithfulness, Sec. 4.2) thus corresponds to the causal graph in the standard sense. Given an SCM M with a context indicator R, there are multiple meaningful sets of structural equations and multiple meaningful associated distributions on observables, which result in multiple observable graphs associated to the SCM M ; an overview is provided in Table 1 in the appendix, see also Fig. 2 and the next subsection. We thus distinguish multiple context-specific graphs. The descriptive graphḠdescr R=r [M ] := G[F M do(R=r) , P M (V |R = r)],captures relations present in the intervened model's mechanisms F M do(R=r) (the SCM in the context r) and encountered in observational data together with this context; the last point is formalized by restricting our attention to the support of P M (V |R = r). For edges not involving R, the same information is contained in G[F M , P M (V |R = r)] using the unintervened F, which supports the intuition that the context-specific graph captures the overlap of observations and the context of interest. Since R is a constant per-context, we add edges involving R from the union graph to Ḡdescr G phys R=r = G[F M do(R=r) , P M (V )]Importantly, the context-specific graphs defined above depend only on the SCM, without making explicit or implicit reference to independence. This sets them apart from graphical representations of independence structure like LDAGs (18) : LDAGS are graphical independence models describing context-specific independence (CSI) structure of a dataset. However, a causal analysis requires an understanding of interventional properties, which, in turn, requires inferring knowledge about the causal model properties. In the single-context case, under the faithfulness assumption, knowledge about the causal properties of models is directly connected to the independence structure. As will be explored in Sec. 2 ). This is a subtle issue: From an independence structure viewpoint, the dependence between X and T exists yet does not seem to make any intuitive sense and does not admit causal semantics. A formal description in which the absence of this edge is considered correct -as it should be for a causal analysis -requires the underlying model as a starting point. Indeed the graphs defined above from the SCM directly (also shown in Fig. 2 ) consider the absence of this edge correct.The support problem A second problem is encountered, because CSI can arise in two different ways from a given model: Either by an explicitly disabled dependence as in the first example, or indirectly through the lack of observational support in a specific context. Note that, in Fig. 2 , our graphical models capture this difference: Once G descr = G phys (left-hand side), while in the other case G descr ̸ = G phys (right-hand side).Concerning the second example in Ex. 3.2 in particular, we make the following observations: (a) There is no evidence for the absence of the link T → Y , e. g. when considering interventional questions. (b) The absence of the link does correctly describe an independence, e. g. edge orientations, cycle resolution, or function fitting on the respective data. Thus, for these problems, the link T → Y is not required for a meaningful description in the context R = 0. (c) R is adjacent to T , but we could add an intermediate node T → M → R without changing the problem, so adjacency of R to the changing link is not guaranteed. (d) We lack observations of another context variable R ′ , which indicates whether T > T 0 , and would suffice to fully describe the qualitative structure of the mechanism at Y . Observations (a) and (b) demonstrate that there is not a single "correct" graph to learn, and the additional flexibility of observable graphs introduced above is necessary. (a) and (b) also make it clear, that it should be carefully specified what should be learned for a particular application and what could be learned (see Sec. 3.1) , and what a specific algorithm will learn (see Sec. 4.2). Finally, CD heavily relies on restricting the search space for separating sets to adjacencies, which becomes non-trivial in sight of (c). We will return to (d) momentarily.In our setup, the context indicator is given by the user. The context indicator can either arise naturally from the problem setting, or it can be the result of preprocessing using anomaly or regime detection methods (see also Sec. F.2). Since the indicator is chosen, it becomes important to understand which context indicator makes for a good choice in view of the problems illustrated above. A core contribution of our work is to understand which assumptions about the model, including the choice of context indicator, prevent drawing wrong conclusions. We give suitable assumptions to make the problem tractable and the result interpretable (in terms of SCMs) in Sec. 3.3 based on observation (d).We now define the assumptions that a proper context indicator must fulfill such that our causal discovery method yields meaningful results. These regime sufficiency assumptions are described below and further formalized in Sec. C; the general approach and the naming are inspired by observation (d) made above and its resemblance to causal sufficiency. anymore. These assumptions mirror the causal sufficiency assumption: They exclude hidden contexts and lead to fewer CITs and a stronger interpretation of results as physical changes. We further discuss the hierarchy of these assumptions in Sec. C.1.We assume causal sufficiency (Def. C.4, no hidden confounders). Furthermore, we require at least weak (descriptive) context-acyclicity (Def. C.5), i.e., G descr R=r is a DAG (acyclic) for all r. In practice, we assume strong (descriptive) context-acyclicity (Def. C.5), meaning there are further no cycles involving any ancestors of R (including R itself) in the union graph. Efficient discovery with our method requires directed cycles of length ≤ 2, as discussed in Sec. D.7, Sec. D.8. We assume that R is a categorical variable which takes finitely many values r with P (R = r) > 0.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?[Yes] Justification: We explicitly state the assumptions our approach made. We add proof sketches for all the claims and theorems in the main paper, and offer the complete proof in appendix. Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
causal discovery with endogenous context variables: 
 We propose a method to discover context-specific graphs and their corresponding union graph which extends the PC(-stable) algorithm (31; 4) with a novel adaptive testing strategy. In Sec.4.2 we will connect the output of our method to the graphical objects defined in 3.1. In particular, under suitable sufficiency assumptions the output is precisely G phys R=r . The approach is not specific to PC but could also be used with other constraint-based algorithms.Formal details are given in Sec. 4.2. First observe that for physical changes (edges in G union not in G phys R=r ) or under weak context-sufficiency (Def. C.3), the only changes in G phys R=r relative to the union graph occur downstream of the context indicator, more precisely in the adjacencies of its children (lemma C.4). On the other hand, spurious links due to selection bias occur upstream of the context indicator (cf. Sec.3.2). By either of the context-acyclicity assumptions introduced above, these two cases are mutually exclusive. Skeleton discovery has to find a valid separating set Z to remove any edge X to Y which is not in the context-specific G phys R=r . By assuming an appropriate version of faithfulness, the decision between independence testing on the context-specific data or on the pooled data does not matter for all other candidates S for Z, as long as we decide correctly for (one) valid S. On the other hand, under causal sufficiency a valid S 0 ⊆ Pa(Y ) exists by a suitable Markov-property. Thus R ∈ S 0 implies R is a parent of Y and we are in the first case, thus should test on the context-specific data. By mutual exclusiveness, we can consistently do so without the risk of selection bias. Conversely R / ∈ S 0 implies R is not a parent of Y , and there is no context-specific information to be obtained here; thus testing on the pooled data is consistent (it does not loose information), avoids selection-bias (X and Y could be ancestors of R) and uses all available data. Thus for the valid S 0 we will by this decision rule always use the correct test. Initialize G R=r as fully connected graph 5:for sepset size = 0 to D -1 do 6:for each j ∈ {1, . . . , D, R} do 7:for each i ∈ Adj(X j ) and S ⊆ Adj(X j ) -{i} with #S = sepset size do 8:TestX i ⊥ ⊥ X j |S \ R, R = r if R ∈ S, else test X i ⊥ ⊥ X j |S on pooled data 9:if independence is found then 10:Remove edge X i -X j from G 11:Orient edges as in PC algorithm 12: return {G R=r | r ∈ {1, ..., n r }} Discovering context-specific graphs First, we discover G R=r which, under the assumptions of Thm. 1, agrees with both the descriptive and the physical graph. For a fixed context value R = r, we adapt the skeleton phase of the standard PC algorithm described in (4) to test on the pooled data whenever R is not in the conditioning set. Whenever R is in the conditioning set, the algorithm testsX ⊥ ⊥ Y |Z, R = r instead of X ⊥ ⊥ Y |Z, R, as described in Algorithm 1.Obtaining the union graph To obtain the union graph, we combine the information about the edges in all context-specific graphs as follows: If the edge from X to Y is found in none of the context-specific graphs, then no edge is added to the union graph. If an edge X → Y was found in any of the context-specific graphs, then the directed edge is added to the union graph. The resulting graph may contain both the edge X → Y and the edge X ← Y . See also Sec. D.14.We study univariate context indicators, but our approach allows multiple indicators. A discussion on the time complexity of our algorithm for this case can be found in Sec. D.9.The Markov property ensures that links that should be removed, are removed from suitable independence relations. Lemma 4.1. Assume strong context-acyclicity, causal sufficiency and single-graph-sufficiency. If X and Y , both X, Y ̸ = R, with Y / ∈ Anc descr R=r (X) (this is not a restriction by context-acyclicity, rather fixes notation), are not adjacent in G descr R=r , then there is Z ⊆ Adj descr R=r (Y ) such that: If Y is neither part of a directed cycle in the union graph (union cycle) nor a child of R in the union graph (union child), then (a) X ⊥ ⊥ Y |Z, otherwise (b) X ⊥ ⊥ Y |Z, R = r. Remark 4.1. If either X = R or Y = R,then there is no per-context test available. In this case, we fall back to the standard result: If X and Y are not adjacent in ACycl(G union ), then there isZ ⊆ Adj descr R=r (Y ) or Z ⊆ Adj descr R=r (X) such that X ⊥ ⊥ Y |Z (see Sec. D.4).The full proof can be found in Sec. D.3. The idea is as follows: Case (a) is relatively standard and follows by a "path-blocking" argument (37) . We focus on case (b), proved in three steps:(i) Assuming single-graph-sufficiency, G descr R=r agrees with the graph that would have been observed if R had been intervened to r. This graph is, in a standard way, associated to the intervened SCM. Thus, the usual Markov property holds (20, Prop. 6.31 X r ⊥ ⊥ Y r |Z r , R r = r implies X ⊥ ⊥ Y |Z, R = r.The faithfulness property ensures that links that should remain in the graph, do remain in the graph and are not erroneously deleted due to ill-suited independence relations. A probability distribution P is faithful to a DAG G if independence X ⊥ ⊥ P Y |Z with respect to P implies d-separation X ⊥ ⊥ G Y |Z with respect to G. This means that, if G ′ ⊆ G is (strictly) sparser, then faithfulness to G ′ is (strictly) weaker than faithfulness to G. Now, G descr R=r ⊆ G union , so "P M (. . .) is faithful to G descr R=r " is weaker than the standard assumption "P M (. . .) is faithful to G union ". Similarly (excluding links involving R), ḠdescrR=r is sparser than what one would expect for a graph of the conditional model (there are no edges induced by selection bias Ḡdescr R=r ), so "P M (. . . |R = r) is faithful to ḠdescrR=r " is also weaker than expected. The hypothesis of the following lemma is thus unconventional but not particularly strong. This lemma is proved in Sec. D.5. Lemma 4.2. Given r, assume both P M is faithful to G descr R=r and P M (. . . |R = r) is faithful to Ḡdescr R=r (we will refer to this condition as r-faithfulness, or R-faithfulness if it holds for all r). Then: a ij X j + c j η j , where i = 1, . . . , D + 1, at a desired sparsity level s. One of the nodes is randomly selected as the context indicator R. To obtain categorical context indicators for the n contexts -1 context-specific SCMs, thresholding is applied by taking the n contexts -1 quantiles of the continuous values of R, adjusted using a data imbalance factor b. To generate a context-specific SCM, n change links of the base graphs are edited sequentially by selecting a variable randomly (excluding the context indicator), adding, removing, or flipping an edge to this variable and ensuring it is a child of R by adding an edge from R to it (and in the case of edge flips, also to the node at the other end of the edge). These operations lead to a non-linear relationship between the context indicator and the other variables. We enforce the presence or absence of cycles in the union graph. The data generation approach is detailed in Sec. E.1.∃Z s. t. X ⊥ ⊥ Y |Z or X, Y ̸ = R and X ⊥ ⊥ Y |Z, R = r ⇒ XWe compare our method, named PC-AC (AC is short from "adaptive context"), to a variant where PC-stable is applied on each context-specific dataset separately by masking, denoted by PC-M. We expect PC-M to have a higher false positive rate (FPR) due to selection bias, and a lower true positive rate (TPR), as links involving the context indicator cannot be found using the masked data. We also compare our algorithm to PC-stable applied on the pooled dataset, named PC-P, which only recovers joint information, i.e., union graphs. We expect PC-P to have a slightly higher TPR, as it tests all links on all data points. PC-P should not suffer from selection bias, thus, the FPR of PC-P should be lower than that of PC-M for acyclic union graphs. For cyclic union graphs, PC-P only finds the acyclification, thus, the FPR should asymptotically remain finite. We also compare to a baseline method, named PC-B: For each context, we compute the intersection of the links found by PC-M and PC-P and then add the links found using PC-P for R. PC-B should converge asymptotically toward the true graph (see Sec. D.10). However, as PC-B runs more tests overall, PC-AC should have better finite-sample properties. As we do not control for cycle length, the performance of PC-AC might be affected by large cycles, as discussed in Sec. D.7 and Sec. D.8. We obtain union graphs for PC-M and PC-B using the unionization approach from Sec. 4. We also compare union graphs from two methods which test on pooled data: FCI-JCI0 (32; 15) and CD-NOD (39) . We do not compare to other methods that obtain context-specific graphs, as they assume exogenous context indicators. A fair comparison to our algorithm is hard to realize due to our assumptions and a lack of available methods suitable to this problem. We compare computational runtimes in App. E.4.We present results where D + 1 = 8 and density s = 0.4. For all j, η j ∼ N (0, 1). We draw the non-zero a i,j randomly from {1.8, 1.5, 1.2, -1.2, -1.5, -1.8}, and set c j = 1 for all j. For the context indicator R, we set c R = 0.2. We discuss results with and without enforced cycles. For the configurations without cycles, we apply only the remove operation (see Sec. E.1). The dataset presented here is balanced in the number of samples for each r, i.e., b = 1. A challenge of our experimental setup is that non-linear relationships between continuous variables and the categorical R must be tested (see also App D.11). Non-parametric permutation-based tests for mixed-type data with non-linear relationships exist (13; 38; 21) , but are computationally expensive. Therefore, we apply a parametric mixed-type data test (35) for most experiments. For a fair comparison, we run versions of the algorithms with so-called link assumptions, where the adjacencies of R are either partially or fully known, i.e., we input the (directed) links to or from R (adj. R-all), or the (directed) links from R to its children (adj. R-ch.), or test without known links (adj. none). For all tests involving only continuous variables, we use partial correlation. When testing pooled data, the link from a regressed-out variable might be missing in one of the contexts, but partial correlation tests should still work, as discussed in Sec. D.12. We show results for "adj. R-ch." and "adj. none" and postpone results with "adj. R-all" to Sec. E.3. Here, we discuss the TPR and FPR for the union and context-specific graphs (averaged across contexts) but also evaluate edgemark precision and recall in Sec. E.3. Error bars are generated using Bootstrap with 200 iterations. We repeat each experiment 100 times, however, the graph generation is not always successful. We evaluate how well our method scales, e.g., with different D, s and b in Sec. E.3. The code, based on causal-learn (40) and Tigramite (25) , and details for replication (random seeds) can be found here: https://github.com/oanaipopescu/adaptive_endogenous_contexts. The left plot of Fig. 3 shows results for the experimental setup described in Sec. 5.1 without cycles. For individual contexts, our method outperforms PC-M in terms of FPR for both link assumption cases. Without any link assumptions, the FPR for both PC-B and PC-AC increase significantly. We test whether this is due to the testing problem (Sec. D.11) by using a non-parametric CIT (21) . Results in the right plot of Fig. 3 show that the FPR is significantly lowered for all methods. For PC-M, the union graph FPR is high compared to the context-specific FPR, as, for the individual context graphs, the links between context indicator and system variables cannot be found (as R is a constant). Still, false positives can appear in the union graph due to selection bias. Our method performs best in terms of TPR and FPR on union and context-specific graphs with context children link assumptions and without link assumptions. PC-B performs slightly worse due to its finite-sample properties. However, with fully known context-system links (Fig. 4 in Sec. E.3.1), masking finds more links overall, having the best TPR, but also higher FPR, followed by PC-AC, PC-B, and PC-P. Results for the same setup as in Sec. 5.1 with enforced cycles in the union graph (Fig. 6 in Sec. E.3.2) show similar behavior of the TPR and FPR, albeit some differences are less pronounced, possibly since the complexity of the problem increases. As expected, the FPR for PC-P, which detects only the acyclification, increases. We observe that for all versions without link assumptions, the FPR increases significantly compared to the setup without cycles. Larger differences in FPR between PC-M and PC-AC can be seen for samples up to 2000. Regarding TPR, PC-AC and PC-B perform best for the context-specific graphs, especially without any link assumptions. Results for the comparison with FCI and CD-NOD for the union graphs, detailed in Sec. E.3.3, show that our method slightly outperforms CD-NOD, which assumes acyclicity. Our method outperforms FCI in terms of TPR, but FPR is lower for FCI.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: We propose a version of the PC-algorithm, for which we offer pseudocode in the paper. We generate our own synthetic datasets. We discuss the data generation model at length in Sec. E.1, and report the parameters used for our experiments in Sec. 5.2. For exact reproducibility, we submit the code as a supplementary .zip file. The details for reproducing the experiments are presented in the readme file. Upon acceptance, we will publicly release the code. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
open access to data and code: 
 Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: We have submitted the code as a supplementary .zip file, and will release it publicly upon acceptance. In the .zip file, all code for the following steps needed to reproduce the results is available:• the implementation of all algorithms discussed in the paper • the data generation • the parameters for the different experimental setups (as a YAML file per setup)• the code for running the results both in a sequential and parallel (cluster) setup • the code for obtaining the metrics and plotting the results Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for a new open-source benchmark). • The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https: //nips.cc/public/guides/CodeSubmissionPolicy) for more details. • The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc. • The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only a subset of experiments are reproducible, they should state which ones are omitted from the script and why. • At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable). • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 
discussion and conclusion: 
 In this work, we have introduced a novel method for causal discovery (CD) of context-specific graphs and their union in the presence of a (possibly) endogenous context indicator. We propose a simple modification of the PC algorithm, but our method can be combined with any other constraint-based CD algorithm. Our approach provides an elegant and efficient solution to the problem of selection bias. Further, we discuss how to formally link context-specific independence (CSI) relations to structural causal models (SCMs); to this end we introduce novel graphical objects and assumptions and prove the soundness of our algorithm. We discuss the theoretical and practical limitations to recovering the true graphs, e. g. in presence of large cycles (D.7). Our method has a lower FPR than applying CD on each context individually, which can suffer from selection bias; it displays better finite-sample performance than a proposed baseline method. Our method is comparable with CD-NOD (39) and FCI (28) on the union graphs, but additionally provides context-specific information.Limitations Our method requires that the context indicator R is observed and chosen appropriately, as discussed in Sec. 3.2; it can handle some union cycles, but not reliably so for large union cycles, see Sec. D.8. Without prior knowledge about adjacencies of R, our method has difficulties to significantly outperform pooled methods and to keep up with the baseline in numerical tests; this is likely due to the CIT problem discussed in Sec. 5.2, an interpretation supported by a smaller experiment with a non-parametric CIT (Fig. 3 , right panel). We thus believe that the problem can be solved using appropriate CITs. Finite sample errors incur further propagated errors similar to standard CD; errors involving links to R additionally lead to a regression to standard CD concerning context-specific information (see App D.13). In this work, we have assumed causal sufficiency, however, real-world data often contains hidden confounders.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results? Answer: [Yes] Justification: We have presented the parameters used for the data generation model in Sec. 5.2. Additionally, all parameters for each experimental setup can be found in individual YAML files in the .zip file accompanying this paper, with additional information in the readme file. Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about statistical significance of the experiments? Answer: [Yes] Justification: Yes, the paper reports error bars. As described in Sec. 5.2, we generate error bars for all results presented in this paper using bootstrapping with 200 samples. For some of the results, especially for the true positive rate (TPR), the error bars are small, and thus cannot be distinguished from the dots indicating the result. Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: We report the computational setup on which we ran our experiments in Sec. E.3. Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 
code of ethics: 
 Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: We mention the licesenses under which the code is released in the code files provided in the .zip file. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [NA] Justification: As our work is not tied to particular applications, we do not believe that this question applies to it. Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)? Answer: [NA] Justification: As our work only uses synthetical datasets, we do not believe this question applies to it. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: We have mainly focused on skeleton discovery and omitted a detailed study of orientation rules. Algorithms such as FCI and CD-NOD already perform well for edge orientations of the union graph, as small experiments in Sec. E.3.3 indicate. But there also seem to be additional edge orientations available only by combining per-context information (see Sec. F.1 for an example), a possibility also mentioned in (9) ; these orientations can be beyond what can be concluded even with JCI arguments. In practice, contexts are often assigned by thresholding, so an extension to contexts given as deterministic functions of observed variables (see Sec. F.2)is of interest. Similar remarks apply to anomaly-detection outputs, defining normal and anomalous contexts. The assumption of causal sufficiency simplifies the arguments, but does not seem strictly necessary; it likely can be relaxed in future work. Finally, our method is in principle extendable to the time-series case. Particularly when allowing for contemporaneous relations of the context indicator, selection bias is a non-trivial problem for time-series.A Delimitation from LDAGs LDAGS (18; 9) are conceptually closest to our work. We thus delimit our work from LDAGs as follows: First, we derive a connection from context-specific graphs to SCMs; also for identification from data, only CSIs are insufficient in our case, a combination of CSIs and independence-statements on the pooled data is required. Second, LDAGs have been studied only under the assumption of acyclicity, while we allow for cycles in the union graph. Third, the method of Hyttinen et al. (9) only accepts categorical variables, while we enable continuous and discrete variables. Fourth, our method reduces the number of tests compared to Hyttinen et al. ( 9), as we do not test each CI relation conditioned on the context. Instead, we adaptively select whether to test on the context-specific or the pooled data. Finally, we distinguish system variables from the context indicator, whereas in LDAGs every variable is treated as a potential context indicator. From a theoretical perspective, this is advantageous, as it increases information content. However, from a practical perspective, the LDAG approach increases the number of required tests.As described near the end of Sec.3.2, in our case, the context variable is chosen from the set of variables. This choice is primarily guided the regime-sufficiency assumptions. This allows for a lot of flexibility: The context variable can be selected from or introduced into the dataset by expert knowledge or any regime detection algorithm; it simply has to behave like a random variables from the perspective of the CD-algorithm. Context variables as such are an orthogonal concept: Given multiple distinct datasets over the same set of variables, the context injects the information about the specific associations in each dataset to the analysis by appending a pseudo-variable to the dataset. The value of this pseudo-variable typically is simply a dataset index. This context could be functional but in general both contexts and functional nodes co-exist independently of each other. Now, we are interested in context-specific changes similar to Example 3.1, which (see above) can be represented as soft interventions, thus as F-variables. This means, we can formally describe these as context-variables which are F-variables. This allows a different interpretation of the strong regime-sufficiency assumption (cf. App.C): It demands, that the context variable belongs to a specific class of functional nodes, which implement a soft intervention removing some of the parents' effects on certain variables.In comparison to JCI (15) , our method reveals not only the position of this context F-variable in the graph; it also reveals some of the internal structure of some functional nodes (those satisfying the reformulated strong sufficiency assumption): It additionally uncovers the graphical changes induced by the implemented soft intervention.Finally, temporal regimes were an important motivation for our study of context-specific behavior. To make the connection to non-stationary time series, these can be formulated using context variables. For the context to be categorical (which we and others assume), the non-stationarity cannot be a drift etc. but must be driven by a regime change, i.e., there must be temporal regions sharing a model or context. By finding the contexts in which the time series is stationary and describing the causal relationships for those stationary parts of the time series, we can describe the non-stationarity system as a multi-context system.Here, we provide details on the three context-sufficiency assumptions discussed in Sec. 3.3 in the main text. Then, we show that these are indeed hierarchical: Strong context-sufficiency ⇒ singlegraph-sufficiency ⇒ weak context-sufficiency. Finally, we give examples of when the other direction is violated.First, we quickly summarize the problem to solve and the ideas behind the different conditions: If the observational support given R = r of, say, X, is restricted to a region where a mechanism f Y (X, . . .) is constant in X, then this can make the link from X to Y disappear in independence testing. Excluding CSIs arising like this requires a way of saying what "constant in a region" means.The single-graph-sufficiency assumption will capture this requirement on a graphical level. This is powerful for proofs, but does not provide much in terms of explanation and might be hard to verify given a specific model. Example 3.1 shows that there clearly is a rather large and intuitive class of models which look like f Y (X, R, . . .) = 1(R)g(X, . . .) + h(. . .). Clearly, one still has to ensure that g or h are "constant in a region", so that only the indicator function can remove links. The strong sufficiency condition will capture this intuitive notion, however, we note that this is only one possible formalization (see Rmk. C.1).As indicated above, the strong sufficiency condition requires a formal way of saying a map is not constant on any region. One possible way to do so is to require it to be injective (in every argument, after fixing other arguments): Definition C.1. We call an indicator R strongly context-sufficient for an SCM M , if both (a) the f i are continuous, such that:(i) If R /∈ Pa(X i ) in the mechanism graph, then f i injective in each argument, after fixing the other arguments to an arbitrary value.(ii) If R ∈ Pa(X i ) in the mechanism graph, then f i is of the form f i (X, X ′ , R, η i ) = 1(R ∈ A) × g i (X, X ′ , R, η i ) + 1(R / ∈ A) × g ′ i (X, R, η i ), with A ⊆ X R and g i , g ′i injective in each argument, after fixing the other arguments to an arbitrary value. This could also be generalized to a sum over multiple indicator functions without changing the results.(b) for every set of variables V not containing R, the distributions P (V ), P (V |R = r) and P (V | do(R = r)) have continuous densities. Remark C.1. The notion of "non-constant on any region" could be formalized differently, for example through smooth mechanisms and a requirement on gradients. Further, this is not supposed to be particularly general, but rather to be intuitive. More general conditions to ensure that detected changes are physical is single-graph-sufficiency or that results are (descriptively) complete is weak context-sufficiency (see below The efficiency of the algorithm proposed in Sec. 4 comes from the restriction that context-specific changes must occur in children of the context R. This requires a weaker condition which we formally describe as weak regime-sufficiency. In this case, the difficult and interesting question to ask is: How should the output of the algorithm be interpreted under this assumption? This is discussed in Rmk. D.2, which can be interpreted as a robustness result: Even when violating the context-sufficiency conditions, the output of the algorithm is descriptively correct, and, given weak context-sufficiency, it is still complete. Definition C.3. We call a context indicator R weakly context-sufficient for an SCM M , ifR / ∈ Pa union (Y ) ⇒ Pa descr R=r (Y ) = Pa union (Y ). Remark C.3.Under this condition, changes may be descriptive only, but the results are complete (all edge removals are found), see Rmk. D.2.Here, we show that above-mentioned assumptions are indeed hierarchical: Strong context-sufficiency ⇒ single-graph-sufficiency ⇒ weak context-sufficiency. We will start from Def. C.1. The underlying idea for requiring (a) is rather direct:Lemma C.1. If X ∈ Pa(Y ) in G[F], and f Y ∈ F is injective in X after fixing all other arguments, then given S ⊆ X Pa(Y ) with x ̸ = x ′ ∈ X X and pa -∈ X Pa(Y )-{X} such that (x, pa -), (x ′ , pa -) ∈ S, then f Y | S is non-constant in X.Proof. f Y is injective in X after fixing its other arguments to pa -, thus f Y | S is injective in X after fixing its other arguments to pa -, so (x, pa-) ̸ = (x ′ , pa -) implies f Y | S (x, pa -) ̸ = f Y | S (x ′ , pa -).Thus, as long as the support S is never singular, i.e., as long as x ̸ = x ′ as above exists, restricting to S does not change the parent sets. The continuous densities required in (b) are an intuitive possibility to ensure this:Lemma C.2. If X ∈ Pa(Y ) in G[F] and a distribution Q is given with Q(X) and Q(Pa(Y )) possessing continuous densities, then ∃x ̸ = x ′ ∈ X X and pa -∈ X Pa(Y )-{X} such that (x, pa -), (x ′ , pa -) ∈ supp(P (Pa(Y )|R = r)).Proof. If Q(X) has a continuous density q x , it is not a point measure, thus there are x ̸ = x ′ with q x (x) ̸ = 0 and q x (x ′ ) ̸ = 0. With X ∈ Pa(Y ) and Q(Pa(Y )) possessing a continuous density q pa , we know q x = Pa(Y )-{X} q pa (pa -, -)d pa -is the marginalization, so 0 ̸ = q x (x) = Pa(Y )-{X} q pa (pa -, x), which implies ∃ pa ′ with q pa (pa -, x) ̸ = 0. By continuity thus (x, pa -), (x ′ , pa -) ∈ supp(P (Pa(Y )|R = r)).Requiring strong context-sufficiency ensures that the only non-injectivities come from R, such that after fixing R in F do(R=r) and ignoring edges involving R (i. e., using the barred Ḡ versions):Lemma C.3. If M = (F, U, V ) is strongly context-sufficient and R = r fixed, then the following are equivalent for X ̸ = R ̸ = Y : (i) X ∈ Pa(Y ) in G[F do(R=r) ] (ii) X ∈ Pa(Y ) in G descr R=r [M ] (iii) X ∈ Pa(Y ) in G phys R=r [M ] (iv) X ∈ Pa(Y ) in G CF R=r [M ] Corollary C.1. If M is strongly context-sufficient, then it is single-graph-sufficient.Proof of the Corollary. By the lemma, Ḡdescr R=r = Ḡphys R=r = ḠCF R=r . The unbarred versions insert the same links involving R from the union graph, thus also G descr R=r = G phys R=r = G CF R=r . Recall that this is the definition of single-graph-sufficiency.Proof of the Lemma. First note that (ii) ⇒ (i), (iii) ⇒ (i) and (iv)⇒ (i) are trivial (if a restriction of f i is non-constant in X, then f i is non-constant in X).Further, f i which do not depend on R are already injective, and f i which depend on R are of theform f i (X, X ′ , R, η i ) = 1(R ∈ A) × g i (X, X ′ , R, η i ) + 1(R / ∈ A) × g ′ i (X, R, η i ),such that the g i , g ′ i are injective. Thus, depending on the value r of R we fixed, eitherf i = g i or f i = g ′ i is in F do(R=r) , both of which are injective. So all f i ∈ F do(R=r) are injective. "(i) ⇒ (ii)": Apply lemma C.2, to Q = P (V |R = r) (as a distribution over V ), to get ∃x ̸ = x ′ ∈ X X and pa -∈ X Pa(Y )-{X} with (x, pa -), (x ′ , pa -) ∈ S := supp(P (Pa(Y )|R = r)). Thus lemma C.1 applies to yield f Y | S is non-constant in X, thus by definition of G descr R=r [M ], also X ∈ Pa descr R=r (Y ). "(i) ⇒ (iii)": As (i) ⇒ (ii) above, but use Q = P (V ) in lemma C.2. "(i) ⇒ (iv)": As (i) ⇒ (ii) above, but use Q = P (V | do(R = r)) in lemma C.2.Next, recall that we had defined R to be (weakly ) context-sufficient: if Y is a not a child of R, then X / ∈ Pa descr R=r (Y ) ⇒ X / ∈ Pa union (Y ).The assumption about distributions seems, especially in light of Lemma C.2 (which is the only implication of this assumption that is needed), relatively weak. However, there are some intuitive examples where we observe that it can be easily violated, and why this is a problem: Example C.2. Assume a root variable X has a noise term, i.e., mechanism kernel (f X ) * η X , with randomness arising (within measurement precision) from a process involving sunlight, e. g., X might be the photosynthesis rate or the solar power provided. Further, assume that R indicates whether it is day or night. Then P (X|R = r) is, for some values of R during the night, singular at a specific value, violating part (b) of the strong sufficiency assumption.On the pooled data, this may not pose an issue. However, on the night-time values of R, nothing will depend on X. This is descriptively correct, but it holds for all children Y of X, irrespective of whether Y are also children of R. Therefore, R is not even weakly context-sufficient. Our algorithm, which relies on weak context-sufficiency, may not remove these links. This kind of problem may be less severe than a violation of part (a) of strong context-sufficiency because it may be testable beforehand.Single-graph-sufficient, but not strongly context-sufficient: Part (b), as stated cf. Lemma C.2, however, can be violated simply by constructing densities with "jumps," which does not interfere with any of the other conditions. More interesting seem violations of part (a). However, this can also easily be done because the f i need not be injective, only "sufficiently injective" to be non-constant on the supports. For example, Y = X 2 + η Y is non-injective but is unlikely to be constant (indeed, this would only happen for P (X) having singular support exactly at two points x 0 and -x 0 , thus violating part (b)). R=r (Y ) (as before), while (also as before) T / ∈ Pa descr R=r (Y ). In the intuitive sense, weak context-sufficiency means that descriptive changes only occur at mechanisms that involve R, but R may, as in the example above, not directly cause the change of the mechanism.Definition C.4 (Causal Sufficiency). A set of variables V is causally sufficient relative to an SCM if, within the model, every common cause of any pair of variables in V is included in V , i.e., there are no hidden common causes or confounders outside V that affect any pair of variables in V . Here, we give proofs and additional information concerning the theoretical properties of the algorithm proposed in Sec. 4.As a systematic overview over the graphical objects introduced in Sec. 3.1, these objects are summarized and compared in Table 1 .The proof of the Markov property in Sec. D.3 will require an additional graphical object, the counterfactual graph. As it is only needed as formal objects, it is introduced here rather than in the main text. Furthermore, the proof of the Markov property should be comprehensible without detailed knowledge of these graphical objects beyond the understanding that there is, in general, more than one (cf. points (a) and (b) in Sec. 3.2). We advise the reader to first read Sec. D.3, and refer here for formal details if required.The counterfactual graph, defined as follows:G CF R=r [M ] := G[F M do(R=r) , P M (V | do(R = r))] = G union [M do(R=r) ]The counterfactual graph describes relations that we would have observed had R been set to r under the same noises. Refer to Sec. D.2 for a more detailed discussion. Generally G do(R) [M ] ̸ = G[F M do(R=r) , P M (V | do(R = r))], even for edges not involving R. This is because the support of P M (V | do(R = r)) may be different (potentially even larger) than for P M .The definition of the intervened model in Sec. 3 is the SCM M do(A=g) := (V ′ , U, F do(A=g) ) with U distributed according to P η of M . At a per-sample level, i.e., drawing a value ⃗ η from P η and running it through both M and M do(R=r) to obtain X i (⃗ η) and (X i ) r (⃗ η) yields exactly the definition of potential responses given in (16, Def. 7.1.4 (p. 204)). Often, the joint distribution P (X i , (X i ) r ) (induced from P η by the procedure described above) is of interest. This would mean counterfactually observing both M and M do(R=r) at the same time. Here we only need X i or (X i ) r variables separately, so the "distribution-level" observables (defining M and M do(R=r) on the same P η , but not on the same realizations ⃗ η) suffice.From the definition of X i (⃗ η) and (X i ) r (⃗ η) the consistency property defined in (16, cor. 7.3.2 (p. 229)) is relatively easy to see: If, for a particular ⃗ η, one obtains R(⃗ η) = r from M , then solving the equations of M do(R=r) yields the same result (they only replace the equation for R(⃗ η) by R r = r, which for this ⃗ η makes no difference). So P ((X i ) r |R r = r) = P (X i |R = r). By definition R ∈ Z ′ , so we can rewrite this as Z ′ = Z ∪ {R} with R / ∈ Z. Now this reads X r ⊥ ⊥ Y r |Z r , R r . Generally, P (X r |R = r ′ ) is not identifiable, however, if r = r ′ it is simply P (X|R = r) by the consistency property of potential responses (16, cor. 7.3.2 (p. 229)) (intuitively, under conditions which lead to R = r anyway, if we had intervened to set R = r, that would not have resulted in any changes in the values of the other variables, and therefore, the conditional independencies). ThusX r ⊥ ⊥ Y r |Z r , R r ⇒ X r ⊥ ⊥ Y r |Z r , R r = r ⇒ X ⊥ ⊥ Y |Z, R = r.We use the definitions of σ-separation (2, A.16) and the induced Markov property (note that in the reference Markov property refers to σ-separation implying independence, while our definition refers to non-adjacency implying independence or CSI) (2, A.21). The following path-blocking argument is rather standard (see e. g. ( 37)), however, we formulate the hypothesis slightly different for ease of application. Proof. Let γ be an arbitrary path from X to Y in G. If the last node W n along γ before reaching Y itself is in Z, then γ is blocked, because the last edge is then the edge W n → Y from a parent W n to Y pointing at Y , such that W n is not a collider along γ. Thus γ is blocked in the d-separation sense, but W n also points at Y along γ with Y in a different strongly connected component as W n (because Y is not part of a directed cycle, its strongly connected component is sc(Y ) = {Y }), so it is also blocked in the σ-separation sense.Otherwise (if the last node W n / ∈ Z), W n must be a child of Y , andγ ends as γ = [X . . . W n ← Y ]. Since Y /∈ Anc G (X), the edges along γ cannot all point towards X. Let W k be the node to the right of the last (closest to Y ) arrow pointing to the right γ= [X . . . → W k ← . . . ← Y ], then W k is a collider along γ.But the part of γ between W k and Y is a directed path from Y to W k (because W k was chosen after the last arrow pointing to the right), so W k ∈ Desc G (Y ) and so are all descendants of W k , Desc G (W k ) ⊆ Desc G (Y ). Thus, Desc G (W k ) ∩ Z = ∅ and γ is blocked both in the d-separation and σ-separation sense.The hypothesis of the following lemma is motivated in the main text in Sec. 4 s. t. X ⊥ ⊥ Y |Z or X, Y ̸ = R and X ⊥ ⊥ Y |Z, R = r⇒ X and Y are not adjacent in G descrWe will refer to this condition as r-faithfulness, or R-faithfulness if it holds for all r.Proof. The statement is symmetric under exchange of X and Y , so it is enough to show X / ∈ Pa descr R=r (Y ). We do so by contradiction: Assume X ∈ Pa descr R=r (Y ) and let Z be arbitrary. Z can never block the direct path X → Y , so they are never d-separated X ̸⊥ ⊥ G descr R=r Y |Z. By (the contra-position of) the faithfulness assumptions, X ̸⊥ ⊥ P Y |Z and if X, Y ̸ = R also X ̸⊥ ⊥ P Y |Z, R = r (the second statement is by definition the same as X ̸⊥ ⊥ P (...|R=r) Y |Z).Thm. 1. Assume causal sufficiency, single-graph-sufficiency, r-faithfulness and strong contextacyclicity with minimal union cycles of length ≤ 2 (edge-flips). In the oracle case (if there are no finite-sample errors during independence testing), the algorithm recovers G descr R=r = G phys R=r .Proof R=r is acyclic, at least one of the edges in a cycle must vanish for each context. Because each edge must appear in at least one context to appear in G union , for a cycle of length ≤ 2, each directed edge in the cycle must vanish in some context. Thus, already by weak context-sufficiency (implied by the hypothesis, Cor. C.2), all (i. e. both) nodes in the cycle are (union-)children of R. Thus, under the restriction to "small" cycles in the hypothesis, part (a) of the Markov property (Lemma 4.1) applies if R is not a parent of Y .If part (a) applies for a Z with R ∈ Z, then writing (ii) We show: The proposed algorithm runs at least all tests required for (i).Z ′ := Z -{R} by X ⊥ ⊥ Y |Z ′ , R ⇒ ∀r X ⊥ ⊥ Y |Z ′ , R = r,The PC-algorithm (and its derivatives like PC-stable or conservative PC) tests for a link X -Y in an efficient way, i.e., at least all conditional independencies X ⊥ ⊥ Y |Z with Z ⊆ Pa X or Z ⊆ Pa Y and deletes an edge (starting from the fully connected graph) if this independence cannot be rejected. By the same argument, our proposed algorithm tests at least all the conditioning sets required for (i).(iii) We show: If X and Y are adjacent in G descr R=r , then the edge between X and Y remains in the graph.By faithfulness (Lemma 4.2), whenever we find a such an independence, the edge is not in G descr R=r , i. e. irrespective of which (context-specific) independence tests we actually perform, we never delete an edge that is in G descr R=r . Remark D.2. Using the argument of (22, Sec. 4) one can supplement the result of Thm. 1 showing robustness against assumption-violations: Replacing the assumption of "single-graph-sufficient" by "weakly context-sufficient", the algorithm still recovers all descriptive information G descr R=r .R=r with G descr R=r ⊆ G detect R=r ⊆ G phys R=r , i. e. all edge removals are descriptively correct, and at least all physically changing links are removed.Proof. The arguments in (i) for showing that we only have to test on context-specific data if R is adjacent still apply (they only use weak context-sufficiency). Further, by weak context-acyclicity together with strong context-acyclicity, links that are in the G union but not in G descr R=r are never between two ancestors of R in the union graph (union-ancestors). So by the Markov property (22, Prop. 4.3), also the remainder of (i) still works. Part (ii) is still as in the standard PC-case. Part (iii) still works, because our faithfulness result (Lemma 4.2) is already sufficiently general.In this section, we discuss how large cycles, which involve more than two edges, i. e. not arising directly from an edge flip, can lead to problems, as exemplified below. However, this problem may not occur directly as discussed in Sec. D.8).Context R = 0: A valid separating set for X and Y is {W 5 } (the only parent of Y ), so X ⊥ ⊥ Y |W 5 , R = 0. But R is not adjacent to Y (in the union-graph, see below), so this would not be tested by our algorithm (however, see Sec. D.8 for a discussion): Candidates for separating sets are the subsets of adjacencies, and context-specific tests are only considered if R is in the candidate set.W 0 W 1 X W 3 W 4 W 5 Y W 7 Z Z ′ Context R = 1: A valid separating set for X and Y is {W 1 } (the only parent of X), so X ⊥ ⊥ Y |W 1 , R = 1.But R is not adjacent (see Sec. D.8 however) to X, so this would not be tested by our algorithm.W 0 W 1 X W 3 W 4 W 5 Y W 7 Z Z ′Ground-truth labeled union graph Note: σ-separation implies d-separation. So it is enough to exclude d-separation.Evidently the empty set is not a separating set X ̸⊥ ⊥ Y . parents of either X or Y would only work in one context, in the other they become dependent, so X ̸⊥ ⊥ Y |W 5 and X ̸⊥ ⊥ Y |W 1 on the pooled data. For this simple example (without Z, Z ′ ), it seems like conditioning on all adjacencies could work. However, this would lead to other problems. For example, with the additional variable Z conditioning on Adj(X) = {W 1 , W 3 } would open the path [X → W 3 ← Z → W 4 → W 5 → Y ] in context R = 0. (Similar for Y and Z ′ .) W 0 1 W 1 X W 3 W 4 0 W 5 Y W 7 R Z Z ′There is another problem with larger union cycles: Links from R to a union cycle cannot be tested using context-specific data (because links involving R cannot be tested for R = r, as R is in that case a constant), but tests performed on the pooled data will only delete edges not in the acyclification of the union graph. Thus, if there is a child of R in a union cycle -which is always the case given weak context-acyclicity and weak context-sufficiency -then the edges from R to all nodes in the cycle cannot be deleted. Thus, R will always point at (be adjacent to, for skeleton discovery) all nodes in a union cycle, irrespective of whether or not they are children of R in the ground truth.While this is a problem in itself, it seemingly removes the problem of Sec. D.7, so one might be inclined to claim soundness for larger cycles. While this is not entirely wrong, we strongly advise against it, for the following reasons:(i) Adding correct constraints (e. g. from expert knowledge) about missing edges from R to system nodes could break the output.(ii) In practise (on final data) error-propagation behavior is poor.(iii) Large cycles can be treated consistently (if they are likely to occur) with the intersection graph method described in Sec. D.10, at the cost of requiring more tests and therefore generally worse finite-sample performance.It might be possible that (iii) can be used "sparingly": One may observe that the above problem should leave a suspicious trace in discovered graphs: There should be many edges from R to the region of the graph containing the cycle, and all nodes in this region not adjacent to R should be heavily interconnected (because they are not tested on context-specific data, so only the acyclification is found). Thus encountering such a structure may be an indication to test the edges of this region again using the intersection graph method.Pairwise tests are always run on the pool in our algorithm, and all context-specific test are "localized" to near the indicator. This means, our algorithm scales like standard PC for increasing variable count but fixed context-variable count. But also for large (or increasing) context-variable count only the possible values of context-variables "nearby" are taken into account as context, leading to an (approximately) linear scaling: The precise number of tests in each iteration of our algorithm depends on the number links remaining after the previous iteration and is thus difficult to estimate precisely. However, on finite data, dependencies are typically only detectable for a finite number H of hops along a true causal path. Given that the union graph has a finite graph degree D G , each variable is found dependent (is within H hops) of a finite number β ≤ β 0 (D G , H) of other nodes (upper bounded by a function of graph degree D G and maximum number of hops H independent of the number of context variables K and the total number of variables D, i. e. β 0 = O K,D (1)). Every variable is thus after the pairwise tests (pairwise testing is always done on the pool) adjacent to less than β 0 other variables. We assume2 that the density K /D of indicators is small compared to β 0 and each indicator takes at most In comparison, a masking or intersection graph approach, as described in Sec. D.10 has to run all tests for each combination of indicator values, thus scales O(exp(K)), as each indicator can take at least two values, i.e., there are at least 2 K combinations. Thus, our method could scale comparatively well to a larger setups involving many variables and context indicators.EQUATIONFrom the perspective that information in causal graphs is in the missing links, one would expect that running causal discovery once on the masked / context-specific data (all data points with R = r) to obtain Ḡmask R=r (containing edges from selection bias; we exclude R from the graph, as it is a constant) and once on the pooled data (all data points) to obtain G pool (which in the case of an acyclic union graph is not the union graph) one should find all necessary information. Thus, we would expect that the intersection graph obtained by intersecting the obtained sets of edges returns the correct context-specific graphs. Indeed we now also have the theory in place to verify this: Lemma D.3. The intersection of the edges obtained from the pooled and masked results yield the correct descriptive graph:Ḡmask R=r ∩ Ḡpool = Ḡdescr R=r .Proof. This follows directly from the Markov property (Lemma 4. Remark D.3. Generally, the information argument about Ḡmask R=r ∩ Ḡpool works as is, but our results show that this has indeed a meaningful connection to the SCM. This approach even avoids the problems of Sec. D.7. However, it has some severe disadvantages on finite data: Both discovered graphs are individually less sparse than Ḡdescr R=r , and thus more CITs are necessary (and especially more CITs with large separating sets, which are particularly problematic for both precision and runtime). Furthermore, all tests are executed once on the pooled data for G pool and once on the masked data ḠmaskR=r . Thus, twice as many tests would be necessary, even without the lack of sparsity. With multiple context indicators, the intersection graph approach would scale exponentially while our method scales linearly, as discussed in Sec. D.9.We use this approach as a baseline to validate the finite-sample properties of our algorithms. We also discuss how to use the intersection graph approach when large cycles may exist in the union graph in Sec. D.8. While large cycles can possibly pose problems for our method, combining our approach with the intersection graph method could yield the correct result in a more efficient way than only using the intersection graph approach.Testing R ⊥ ⊥ Y is not only difficult due to non-linearities, but also to the following problems:We focus on the case where R is binary for the example Y = 1(R) × X + η Y , and X = η X with both η X , η Y ∼ N (0, 1) standard normal. Here Y ̸⊥ ⊥ R (with the empty conditioning set). This is essentially a two-sample test, i. e., we are asking if P (Y |R = 0) ̸ = P (Y |R = 1). This difference is only visible in variance (higher moments) and not in mean value. P (Y |R = 0) is simply η Y and thus standard-normal N (0, 1), while P (Y |R = 1) is the sum of two independent standardnormal distributed variables (η X + η Y ), so it is N (0, σ 2 = 2). Many independence measures like correlation or generalized covariance measure (30) are insensitive to changes in higher moments only (e. g., testing the correlation of R and Y checks for the non-zero slope of the linear fit through (0, E[Y |R = 0]) and (1, E[Y |R = 1])). Mutual information or distance correlation based tests can gain power against this alternative. Generally, using parametric tests for this scenario (even simple workarounds like testing correlation on squared or absolute values) should also be possible and might provide better power. Some care should be taken (30, Example 1) to formulate the model-class we actually want to test on, as we also have to be able to account for the synergistic (the model is not additive in X and R) dependence on R.Conditional independence testing often relies on regressing out conditions and testing on residuals (e. g. partial correlation). This is usually justified for causal discovery roughly as follows: Only the "true" adjustment set has to work -"erroneous" independence is assumed impossible by faithfulness -and only if X and Y are "truly" non-adjacent. The "true" adjustment-set for X, Y is w. l. o. g. Z = Pa(Y ), with P (Y |Z) = P (Y | do(Z)), so the regression function approximates the structural f Y which exists and is well-defined as a mapping (in particular single-valued); so for additive noise models testing on residuals should work.One might, at first glance, consider that this may fail to hold if some links from Z to Y are contextdependent (thus multi-valued). However, this is not the case: If a link from Z to Y vanishes by a change of mechanism, then R ∈ Z = Pa(Y ) and the multi-variate f Y is not multi-valued (on its multi-dimensional input). If a link from Z to Y vanishes by moving out of the observational support, the previous argument still works, i.e., f Y does not change.Moreover, if the per-context system (dependencies on variables ̸ = R after fixing R = r) is additive (e. g. linear), then the fit can be done additively because, again, if a link from Z to Y vanishes by a change of mechanism, then R ∈ Z = Pa(Y ). Since our method tests on data with R = r only, it learns a consistent restriction of the multi-variate f Y where it is additive. For example, if ∀r : M do(R=r) is linear with Gaussian noise, then conditional independence between two variables X, Y , both ̸ = R, can be tested by partial correlation, and the method remains consistent.Errors involving links to R affect the performance of the algorithm as follows:For false negatives (a link involving R is not found), no CSI is tested, and thus it does not find additional context specific information, but it also does not incur further errors beyond those inherent to PC.In the presence of false positives (a link involving R is detected to a non-adjacent Y ), our methodlike PC -executes more tests than strictly necessary, because more adjacencies are available. If R is erroneously found as adjacent, and is added to the conditioning set, in our case the tests will be run as CSI (for R = r), i.e., using the context-specific samples only. This can lead to testing on a reduced sample size resulting in further finite sample effects. It should not induce selection bias, because we only have to find (at least) one valid separating set, but if R is non-adjacent such a separating set, not containing R, exists, still allowing to correctly delete the link.R=r for all r, we can also construct a labeled union-graph similar to LDAGs (or more precisely -context-specific LDAGs (18, Def. 8) ). This follows the philosophy that one obtains strictly more information than just the union graph, and wants to represent this information as an extension of the graph rather than in separate graphs.There are two viewpoints on how to symbolize multiple graphs in one: The information in the graph (about factorization of the distribution) is in the missing links, hence, link labels indicate when they "disappear" (stratified graphs and LDAGs follow this convention) -or the edges represent mechanism relationships of the underlying SCM and labels represent when mechanisms are "active." The second option seems more straightforward to read (from a causal model perspective), so we follow this convention below:• If there is an edge X → Y in G descrR=r for all r, insert an (unlabeled) edge X → Y into the labeled union graph.• If there is never an edge X → Y in G descr R=r for any r, than there is no edge from X to Y in the labeled union graph.• If there is an edge X → Y for r ∈ R, and no edge for r / ∈ R, with R ̸ = ∅ a proper (i. e. not equal to) subset of the values taken by R, then add an edge labeled by R, e. g.------→ Y , to the labeled union graph.Here, we describe the data generation approach for the numerical experiments in detail. Moreover, we present further results which aim to explore the robustness of our method under different scenarios.We generate the simulation data using SCMs. We start by constructing what we will call the base SCM / graph, a linear SCM with up to D + 1 nodes, and with acyclic (mechanism) graph, of the formEQUATIONwhere i = 1, . . . , D + 1, Pa(X i ) are the parents of X i , c j are constants, and η j are exogenous noise terms. The function f i (Pa(X i )) is defined as:EQUATIONWe use a desired sparsity level s ≤ 1, which leads to the number of total links L (non-zero coefficients a ij ) to be found as the fraction s of the total possible number of links D(D + 1):EQUATIONAs a next step, we select a variable to become the context indicator R. To this end, we randomly select an index k ∈ {1, ..., D + 1}. Then, we assign the variable X k as the context indicator, i.e., R = X k . For example, if k = 2, then R = X 2 .As discussed in Sec. 3, R is a categorical variable that indicates which of the n contexts contextspecific SCMs the respective sample belongs to. Thus, we must generate categorical values from the continuous values of X k . To achieve this, we define n contexts -1 quantiles {q 1 , . . . , q ncontexts-1 } as linearly spaced values between 0 and 1. These quantiles are then adjusted using a balance factor b to compute the final quantiles as {q b 1 , . . . , q b ncontexts-1 }. The balance factor b affects the distribution of data points across contexts:• If b = 1, the data points are balanced across contexts.• If b < 1, the balance tips towards the left tail of the distribution.• If b > 1, the balance tips towards the right tail of the distribution.Using the adjusted quantiles, we compute the values for R as follows:R =         1 if X k < q b 1 2 if q b 1 ≤ X k < q b 2 . . . . . . n contexts if X k ≥ q b ncontexts-1Thus, the continuous variable X k is converted into a categorical variable R indicating the context based on the thresholds defined by the adjusted quantiles. This results in R having a non-linear relationship to all other variables.Using R, we can now proceed to generate context-specific SCMs which have n change changed links compared to the base graph. For this, we edit the base graph sequentially using the allowed operation set S ops , as follows:1. Randomly select a variable X c (excluding R),2. Randomly select an operation from the set S ops , 3. Apply the change to the graph:• Add an Edge: Add an edge from X c to another randomly selected variable X c ′ , if an edge between X c and X c ′ does not already exist. • Remove an Edge: Remove an existing edge from X c to a randomly selected adjacent variable X c ′ . • Flip (Reverse) an Edge: Reverse an existing edge from X c to a randomly selected adjacent variable X c ′ , making X c ′ the parent of X c .4. Add a link from R to X c ′ . If the selected operation is flip, then add a link from R to X c and X c ′ .We note that performing the selected operation on the randomly selected node is not always possible. For example, the node might not have any edges to remove, or adding an edge would lead to the graph being cyclic. We, therefore, check for possible inconsistencies and repeat the process several times (in the case of the experiments presented here, 10 times). If an accepted graph is not obtained after these repetitions, we interrupt the process and continue with a new base graph. Therefore, not all repetitions of an experiment may be successful.As described in Sec. 5.2, we can either enforce cycles between the context-specific graphs (and thus, cycles in the union graph) or we do not allow them for a specific configuration. If the configuration must contain cycles, for each proposed edit operation, we check whether doing this operation would result in cycles and accept it if so. Otherwise, we reject the operation and start a new trial. This also leads to some repetitions of an experiment not being successful. We mention that we do not enforce a specific cycle length.Once a context-specific SCM is created, we generate the data according to this SCM. For each variable X i in the context-specific SCM, we determine its value based on its parents and the structural equations of the SCM. Formally, for a variable X i with context-specific parents for context R = r, Pa R=r (X i ), its value is generated as:X i = f i (Pa R=r (X i )) + c i η iwhere f i is defined as in Equation 2. By following this process, we ensure that each variable affected by the context is generated according to the modified SCM, incorporating the context-specific changes introduced during the editing phase.As mentioned in Sec. 5.2, we repeat each experiment 100 times. However, as discussed above, due to the random data-generating process described above, not all graph generation trials are successful. The accompanying 'failed.txt' file in the .zip file reports how many graph generation trials have failed for each configuration. Here, we mention that, for the results presented in Sec. 5.2, 10 trials have failed. also report the number of failed trials in the following subsections which present additional results.We now present results for further configurations, which vary D, s and b, as well as detailed discussions of the experiments mentioned in the main paper. All results were obtained by conducting trials in parallel across a cluster of 116 CPU nodes, each equipped with 2x Intel Xeon Platinum 8260 processors. The supplementary .zip file contains the original code used for running the experiments on the cluster. Additionally, we provide an option to run the experiments in sequential mode. Results with context-system link assumptions Fig. 4 presents the results from the setup presented in Sec. 5.1 without cycles where all links to and from R are known. Knowing these links significantly improves the TPR for all methods. Notably, PC-M has the highest TPR and FPR, indicating that only knowing about the context-system links in the endogenous context indicator setting does not suffice for PC-M to discover the correct graphs. The FPR is generally lower for the other methods, indicating fewer false discoveries due to the reduced uncertainty about context-system connections. Edgemark precision and recall Fig. 5 shows the precision and recall values for both the union and the context-specific graphs (averaged over contexts) for PC-AC, PC-M, and PC-B and the union graphs for PC-P for the different link assumptions. On the context-specific graphs, we observe that assuming all links to and from the context indicator are found (adj. R-all), PC-AC has the lowest precision, i.e., lower risk of falsely oriented edges, and highest recall, i.e., finding the correct orientations, followed PC-B. For the other link assumption cases, our method outperforms all other methods for the context-specific graphs. For the union graph, PC-P scores best, likely because it finds the most links in the union case, followed by our method. PC-M scores last.Fig. 6 presents results for the experimental setup as described on Sec. 5.2 with the difference that cycles between the context-specific graphs are enforced here. For this experiment, a total of 29 trials have failed (see Sec. E.2). From the total of 71 samples, 69 samples have cycle length 2, and 2 samples have cycle length 3. We observe similar behavior of the TPR and FPR for the non-cyclic results presented in Sec. 5.2. However, some differences are less pronounced, possibly due to the increasing complexity of the problem. As expected, the FPR for PC-P, which only detects the aycyclification, increases. We observe that for all versions without link assumptions, the FPR increases significantly compared to the setup without cycles. The largest differences in FPR between PC-M and PC-AC are for a smaller number of samples, up to 2000. For 3000 samples, the FPR of PC-M is lower than the FPR of PC-AC, possibly due to overall fewer links being tested compared to PC-AC: Since PC-AC also finds context-system links, it tests more adjacent pairs and possibly finds more links. However, in terms of TPR, PC-AC, and PC-B perform best for the context-specific graphs, especially without any link assumptions. Skeleton discovery As mentioned in Sec. 5.2, we compare the results of our algorithm for the discovery of the union graph to the Fast Causal Inference (FCI) algorithm (32) , discussed as JCI-FCI0 in Mooij et al. (15) . We apply JCI-FCI0, which does not assume exogenous context variables (we name it only FCI from now on) for the multi-context setup, and CD-NOD (39), a causal discovery algorithm for heterogeneous data, which assumes exogenous context variables. Both algorithms are applied to the pooled datasets. We use the causal-learn implementation (40) for both algorithms. Since, for the available implementation, it is not possible to introduce link assumptions from the beginning of the discovery phase (while it is possible to introduce link assumptions, at the moment, they are only used to refine the graph after the discovery phase), we only evaluate the algorithms for the case without link assumptions (adj. none). As the algorithms only use pooled datasets, we do not present any result for the context-specific graphs.Fig. 8 presents the results for the setup in Sec. 5.2 of the main paper, without cycles (left-most two plots) and with cycles (two right-most plots). We observe that both FCI and CD-NOD behave similarly to PC-AC. We observe minimal differences between PC-AC and CD-NOD for smaller samples, with our method performing slightly FCI performs best across all methods in terms of FPR. Nevertheless, FCI can only recover the union graph, and thus, our method still has the advantage of gaining more descriptive information than FCI. However, our method may perform poorly due to the CIT, and we believe that using an adequate CIT, our method can improve over FCI and CD-NOD.Edge orientations We observe that FCI and CD-NOD perform best when comparing the precision and recall values for the edge orientations between the different methods. This result is expected: Our method focuses on skeleton discovery and does not propose any special rules to orient edges obtained from the union of the context-specific graphs, which is left for future work, as discussed in Sec. 6. CD-NOD and FCI develop specific orientation rules involving the context indicator, and thus are expected to perform better in terms of orientations at the moment. However, as we consider that there is additional information from the context-specific graphs that can be used to orient more edges, as exemplified in Sec. F.1, we believe that the performance of our method regarding edge orientations can be significantly improved. Here, we explore the scenario where the number of changing links between contexts varies. We keep the experimental setup as in Sec. 5.2, only changing the number of changing links n change to 2 and 3. Figure 9 shows that as the number of changing links increases, the TPR generally improves, suggesting that the presence of more changes provides more signals that the causal discovery process can detect. Conversely, a higher number of changing links leads to an increase in FPR, especially for the cases without link assumptions. However, we believe this is due to the CIT, which generally leads to a higher FPR, as the results in Sec. 5.2 of the main paper indicate. Our method performs similarly or outperforms PC-P, PC-M, and PC-B for cases with known links from R to its children. We mention that, for this experiment, a total of 18 trials have failed for n change = 2 and 25 trials have failed for n change = 3 (see Sec. E.2). Here, we examine the impact of increasing the density of the base graph on the performance. We keep the experimental setup as in Sec. 5.2, and only the graph density to s = 0.6. As Fig. 10 shows, with higher graph density, the TPR decreases slightly compared to the results in the main paper, while the FPR increases. This is likely due to the increased complexity of the problem and possible confounding. The denser the graph, the more candidate links the algorithm needs to evaluate, which can increase the likelihood of falsely identifying spurious links as causal. Our algorithm still outperforms PC-M for known links between R and its children. PC-B scores best overall in terms of FPR, but its finite-sample issues lead to lower TPR results, indicating that it finds fewer links overall.For this experiment, we mention that a total of 12 trials have failed (see Sec. E.2).Fig. 11 present the performance of the algorithms across graphs with different N = 6 (left plot) and N = 12 (right plot), while the rest of the experimental setup remains as described in Sec. 5.2. Generally, TPR decreases as the number of nodes increases, possibly due to the increased complexity and the higher dimensional space in which the causal discovery process operates. Larger graphs have more potential causal connections, making it harder to distinguish true causal links from noise. The FPR slightly increases with the number of nodes. However, our method (PC-AC) still outperforms PC-M when no links are assumed, or the links from R to its children are known. For this experiment, we mention that a total of 14 trials have failed for N = 6 and 12 trials have failed for N = 12 (see Sec. E.2). Results investigating the effect of data imbalance, presented in Fig. 12 show that our method, PC-AC, still performs better in terms of FPR or TPR compared to PC-M, despite a data imbalance factor of b = 1.5, if links from R to its children are known or without any link assumptions, similar to the results in the main paper. While the overall FPR increases slightly from the results presented in Fig. 3 , for PC-P, this increase seems insignificant. We mention that, for this experiment, a total of 10 trials have failed (see Sec. E.2).We evaluate the computational runtimes of the different methods evaluated in Sec. 5.1 and App. E.3 and present these in Table 2 . PC-AC has a slightly higher runtime than PC-M due to the additional tests involving the context variable. Among the methods that generate the union graph, the baseline PC-B has, as expected, the longest runtime, as it uses the information from both the PC-M results and the results of PC on the pooled data (PC-P). We also observe that CD-NOD seems to be the most time-efficient algorithm on pooled data, while PC is the most computational intensive. This difference in particular might arise from implementation details. In the Discussion Sec. 6, we have outlined a few ideas for extensions and future directions for the algorithm presented in this work. Here, we discuss these ideas in more detail.In some cases, additional edge orientations beyond JCI arguments (15) are possible by using contextspecific information. For instance, consider the following example:Example F.1. Orientations Beyond JCI: In context 0, we have X → Z ← Y . In context 1, we have X → Y , and no edges at Z. Context 0:X Y Z Context 1: X Y Z JCI / Union: X Y Z R ?Context 0 can be oriented by the unshielded v-structure, context 1 can be oriented by a JCI argument (15) . In the JCI case, both X → Z ← R and X → Y ← R are unshielded v-structures, thus leading to orientations. The v-structures X → Z ← Y and R → Z ← Y are both shielded. Furthermore, the v-structures that would occur for Z → Y , i. e. X → Y ← Z and R → Y ← Z are both shielded and thus could not have been detected, and there is no conclusion that could be drawn from not detecting them. Thus the (pooled) independence structure is symmetric in Y and Z. Therefore, also no graphical consistency criterion (avoid directed cycles etc.) can lead to an orientation of the edge between Y and Z and the graph obtained by JCI is also completely symmetric under exchange of Y and Z, missing the orientation of the edge Y → Z (which in context 0 is found immediately by our method).For the following example,X Y Z R A Ba slightly different problem occurs. Orientations are obtained by that X → Y ← R is unshielded, thus detectable, and in context A, Y ← Z can be excluded, because it would introduce the unshielded v-structure X → Y ← Z (which is used for orientations in context B). This argument works only for the sparse union graph (as opposed to its acyclification) and only per context (the orientation is only well-defined per-context). The latter point makes it fundamentally inaccessible to JCI arguments, the former point also leads to the problem of Example F.1 (where it occurs independently of the second point).In realistic examples, context indicators are often obtained by thresholding. An extension to such deterministic functions of observables seems possible as follows:Given a deterministic indicator C = g(W ) for an observable W (g could be multi-variate, see below), add C as a node to the output graph and draw a special edge W → C (indicating that this edge is deterministic and provided by expert knowledge). Then execute a standard constraint-based causal discovery method on the variables (excluding C) with these modifications for all c values of C occurring in data: EQUATIONY = α + 1(C) × (1 + βW ) × X + η Y .If α = 0, then the mechanism depends qualitatively on C: If C = 0, then Y does not depend on X.If β ̸ = 0, then the mechanism depends quantitatively on W . Justification: The adaptive nature of our method is explained in Sec. 4, the assumptions and connection to SCMs are in Sec. 3, theoretical results are in Sec. 4 and finally, the numerical experiments are reported on in Sec. 5.2.Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.We assume all SCMs are uniquely solvable, i. e. given F and U , there is a unique V solving these structural equations. This is not a restriction for the per-regime acyclic models considered later.Then within the β many adjacencies there is typically at most one indicator; generally we could replace C by C β 0 , as there are at most β0 many adjacent indicators with at most C β 0 many combinations of values, which is still OK,D(1). 


Decision: The paper addresses the challenge of causal discovery in scenarios where causal relationships vary across different contexts, particularly when the context variable is endogenous, meaning it depends on other variables within the system. The authors propose a modified version of the PC algorithm that adaptively chooses between performing conditional independence tests on either pooled data or context-specific data, depending on the context. This approach is designed to avoid the selection bias that can occur when using traditional methods in the presence of endogenous context variables. The proposed algorithm is theoretically grounded, with proofs of its completeness under certain assumptions, and its effectiveness is demonstrated through empirical studies.

The reviewers are unanimously in favor of accepting the paper, though their support is somewhat tentative. There are several issues with clarity (definitions, assumptions, introduction, etc.), which I believe the authors can address during the revision. However, I am also a bit concerned about whether it may require another round of review, as it could involve a nontrivial amount of changes. One of several concerns is about the known context indicator. Given that this paper is the first to address an endogenous context variable, it is reasonable to assume a known context indicator. In the rebuttal, the authors claimed that a regime detection algorithm can be used to identify a context indicator, but no reference or clear explanation is provided on how this could be applied. If it is feasible, adding experiments with a learned context indicator would significantly strengthen the paper.

Given that the paper is well-motivated and the contribution is novel enough, I lean towards recommending the acceptance of the paper. The authors are strongly encouraged to revise the paper for better readability and also to add a discussion on learning a context indicator.

Reviews: review: paper_summary: paper_summary: The authors address the problem of causal discovery in situations where causal relationships change across different contexts. They propose a modified version of the PC algorithm, which either performs a conditional independence test (CIT) on pooled data or context-specific data, depending on the scenario. The paper presents several theoretical results for this modified algorithm and evaluates its performance through various empirical studies.
strength_weakness: strength_weakness: The organization of the paper makes it difficult to follow. For instance, numerous graphical definitions are presented without examples or figures, even in the appendix, which would aid understanding. I strongly suggest at least a running simple graphical example to show different graphs. Additionally, the definitions of the main assumptions are not included in the main text, further complicating comprehension.

The results rely on the assumption of causal sufficiency, but it is unclear how these results could be generalized to scenarios with hidden variables. Another significant limitation of the proposed framework, as acknowledged by the authors, is the requirement for knowledge of the context indicator.

The proposed algorithm require $Adj(X_j)$ and whether R belongs to this set for all $j$. This knowledge is obtain by a standard PC algorithm using pooled data. What happen if the output of this PC algorithm is noisy, e.g., it returns that R belongs to the conditioning set but in reality it is not?

While the proposed algorithm clearly outperforms the baselines in certain settings, the margin of improvement is relatively small, particularly when the conditional independence test is non-parametric.
questions: questions: Please see above.
limitations: limitations: Please see above.
paper_summary: paper_summary: The authors consider an SCM $M$ with a labelled and observed context variable $R$ that is endogenous, i.e., causally depends on other variables in the model. In this setting it is generally not true that $P_{M}(\ldots \mid R=r) = P_{M}(\ldots \mid \mathrm{do}(R=r))$. Therefore, selecting only data for which $R=r$ and applying constraint based causal discovery algorithms - as one might do with exogenous context variables - will introduce selection bias when estimating context-specific graphs $G^{\text{phys}}_{R=r} := G[M \mid \mathrm{do} (R=r)]$. 

The authors propose an algorithm that adaptively selects whether to do a conditional independence test (CIT) using context specific or pooled data. Their method is provably complete (recovering the population graph of the underlying SCM with oracle CITs) under a set of assumptions that relate both to the underlying SCM and to the (context specific) observational distribution.
strength_weakness: strength_weakness: (Presentation) The authors should present their results more carefully by using terms only once they are defined and being consistent with notation. For instance, the motivation behind a "$G^{\text{phys}}_{R=r}$" is considered informally in the discovery goals section (a few pages before the notation is used), but the notation itself is defined neither formally nor informally in the main text. The definition is left to appendix but is not cross-referenced in the main text. Similarly, I have no idea what the "AR" in their method "PC-AR" actually stands for.

The structure of the paper unfortunately makes simple typos (e.g., "$G^{\text{descr}}_{R=r} \subset G^{\text{union}} = G^{\text{union}}$" on line 262) difficult to interpret or correct by the reader. If the authors could restructure the paper to assist the reader in understanding their technical contributions this would greatly impact my evaluation of the submission's presentation. For instance, the physical and descriptive context graphs could have simply been defined in the discovery goals section.

(Interpretation of strong context-sufficiency) It is not clear to me what the authors mean by "injective" in the definition of strong context-sufficiency. I will assume that "$f_i$" refers to the structural equation for variable $X_i$ (the authors did not define this in the text) and take "injective" to literally mean $f_i$ has at most one input that corresponds to any output in $\mathrm{dom}(X_i)$. This seems to me to be a very restrictive assumption: say $f_i$ had one parent $A$ and a noise term $e$ and both were, e.g., normally distributed, wouldn't the "injective" property preclude any model like $f_i(A, e) = A + e$? 

(Experimental design and results) In the experiments, cycle length in the union graph was not controlled despite this being required for the authors' method to recover $G^{\text{phys}}_{R=r}$. This complicates comparison with the brute-force baseline "PC-B". Furthermore, details on the time taken to run PC-B versus PC-AR should have been included to understand the computational advantage of the authors' approach.
questions: questions: (Algorithm 1) There is no definition of the dummy variable $i$ in Algorithm 1. Is it meant to be looped over in line 4? Is line 12 in Algorithm 1 rather meant to say "remove $X_i - X_j$ from $G_{R=r}$"?
limitations: limitations: Yes - the authors make good effort to document the challenges to their approach and motivate potential future work.
paper_summary: paper_summary: This paper proposes a constraint-based algorithm for context-specific causal discovery, which accommodates endogenous context variables.
strength_weakness: strength_weakness: 1. The proposed method assumes that the context variable is both observed and known. In real-world applications, we may not know which observed variable is the context variable. Furthermore, the context variable might be a latent variable.

2. This paper lacks readability. It entails too many definitions and notations. It might be better if the authors could summarize these things in Appendix for reference. Besides, the authors would better defer more details to appendix rather than present almost all results in the main text.
questions: questions: In line 80, the authors claim that "where the context can also be endogenous (in [17])", but in line 89, the authors also claim that "all above-mentioned methods assume that the context variable is exogenous.", which is quite confusing. If [17] allows endogenous context variables, please detail the difference between [17] and this work.
limitations: limitations: Please see weaknesses.
paper_summary: paper_summary: The authors tackle the challenging causal discovery task, namely, causal discovery from the pooled dataset collected under different environments, where the environment (a.k.a., the context) can be dependent on the system (endogenous) variables (i.e., the variables whose causality we are interested in). The proposed algorithm is simple and intuitive: It employs a modified PC algorithm using context-dependent independence tests. Although I am not so familiar with this topic, overall, I have enjoyed reading the paper. However, there are several clarity issues, so I hope that my comments will be helpful for paper revision.
strength_weakness: strength_weakness: Below I will enumerate several clarity issues.

Section 1

* While the first and the second paragraphs in Section 1 are very clear, the third one was a bit disappointing. The authors suddenly introduce several technical terms, such as context-specific models, single causal union graph, and cyclic union graph, etc, without definitions or explications. Readers cannot understand, for instance, which graph in Figure 1 is a « single causal union graph » because there is no description about it.

* « as exemplified in Example 3.1 » Example 3.1 is very far away from Introduction. Thus, this paragraph is not well structured.

Section 3

* Relationship between the notion of intervention and the context indicator is unclear. I am not expert of causal discovery from multiple datasets like JCI algorithm, but I understand that the context is a more general notion than intervention: Some interventions can be represented using context variables, but not vice versa. Am I correct? For instance, how can non-stationary time series be formulated using context variables? Please elaborate the context variables more (in Appendix, if there is not enough space).

* Related to the above, Example 3.1 can be regarded as soft intervention on $Y$? 

* Section 3.1 involving Example 3.2 is very hard to follow, although the example is intuitive. Lines 141-149 are very difficult to follow. Please introduce each variable (e.g., ice cream sale) together with notation (e.g., $Y$). Notation is inconsistent: Sometimes $R=0$, and sometimes $R=« ice »$. This paragraph should be clearer to clarify the significance of this work.

* In Section 3.2, each graph is not clearly explained. For instance, is the mechanism graph « a causal graph in a standard sense »? If so, state it clearly. 

* Many notations are introduced without definition. Examples include $f_i|_{support P(Pa_i)}$, $\mathcal{F}^M$, $P_M$, and $G_{R=r}$ (Section 4.1). This makes it hard to follow the paper. Please proofread the paper before submission.

Section 4

* I am not so familiar with the statistical test for context-specific independence, but could you elaborate the statistical procedure for the testing? Are you simply picking up only the subset of data instances with $R=r$? Is it statistically reliable? Is the test statistic identical with the usual independence testing?
questions: questions: NA
limitations: limitations: NA

Comments: title: Response to Comment by J2yL
Response
Response to authors
comment: We thank the reviewer for pointing out their reference [1], henceforth referred to as [ref1].

To the best of our understanding, the referenced paper approaches the problem from a similar perspective as the work of reference [22] from our paper on LDAGs: The authors of [ref1] seem to be primarily interested in independence structures. The basic idea of the staged trees they use seems to be that for a purely categorical model independence $P(Y=y|X)=P(Y=y)$ simply means that \emph{all} conditionals $X=x$ lead to the same probability, $P(Y=y|X=x_i)=P(Y=y|X=x_j)$. Therefore, one can generalize independence by looking at the equalities of values in conditional independence tables (replace "all" by "some" in the previous sentence). This generalization includes CSI
as a special case. 

While [ref1] has valuable contributions, there are a number points where the approach substantially differs from ours:

* Their core discussion is on (generalized) independence structure. 
* Their approach is inherently tailored to purely categorical systems.
* The relation to causality of the staged trees approach is indirect: In [ref1], SCMs are never discussed. Instead, the authors use a global causal order (which is assumed to be known) to obtain interventional distributions directly from a Markov factorization (a tree-like one along the total order, rather than the one along a graph customary with SCMs).
* Their approach to finding this total causal order (via BIC) should not work generally (eg with continuous variables, see Prop. 4.1 (p.44) from Peters, Janzing, and Schölkopf. Elements of causal inference: foundations and learning algorithms. 
* They cannot use sparsity of the graph for efficiency, while our PC-like approach does.
* They discuss neither causal (SCM-based) interpretations of the results nor the assumptions required to this end, while our work focuses specifically on those challenges.
* They do not specify whether acyclicity is assumed, but the discussion around equation 4 seems to implicitly assume acyclicity.


While we agree that there are similarities to our work, the approach of [ref1] seems less related to our work than LDAGs. This is why we we have discussed LDAGs in detail in the main paper. Furthermore, we would like to point out again that a core difficulty and insight within our work is the study of observational-support-related problems.
We believe treatments focused on independence structures and typical (SCM-based) treatments of causality have very different underlying ideas, strengths and assumptions. To some degree, our results can bridge the gap between SCM based approaches and independence models, which we consider an important contribution of our work.
Many thanks to the authors for their responses, their clarifications and their revisions.

I find the authors' approach quite exciting, especially their insights into challenges for large cycles and unknown endogenous context variables. However, it is difficult to understand from rebuttals alone what impact their (quite heavily) revised presentation will have on the overall clarity of their work. For these reasons, I am willing to revise my score for the soundness of their approach to 3: Good, and my overall score to 6: Weak Accept.
Dear Reviewers, 

Thank you to those who engaged in the discussion with the authors. 

Please confirm that you have read the authors' rebuttal and outline any concerns that remain (`9Xqc`, `DjMB`). Are there aspects of the rebuttal where the explanation is insufficient? If so, please specify. Additionally, please update your score to reflect your current view of the paper, even though you can change it later. This will give the authors a better sense of how your perspective has evolved.

Best,
AC
Thank you for your response. 

I roughly understand what the authors mean, but still, I cannot clearly understand the definition of context variables and related graphical concepts. Since they are strongly relevant to the study, their poor clarity will reduce the impact of this work. For this reason, I will maintain my overall rating, and I hope that the authors can clearly explain these notions more.
Dear Reviewers,

The authors have submitted the point-to-point response. Currently, the reviews are all around borderline accept/reject. 
Could you all read the rebuttal and see if it addresses the main issues raised by your review?  (Thanks J2yL for the engagement.)
Do your ratings change? Are there still some unresolved issues? 

AC
I think the contribution of this work is somewhat limited because

- This work is not the first one investigating causal discovery with context-specific variables. Some previous works such as [1] has also discussed this problem.

- Both detection of context variables and latent context variables are out of your scope.

[1] Context-specific causal discovery for categorical data using staged trees. AISTAT 2023.
We are sorry we could not change your mind. Nevertheless, we do not want to leave the new input unaddressed: After the introduction, we give a list of four main points where the paper makes novel contributions. Indeed one of those points is "[...] adapting existing constraint-based CD algorithms to extract context-specific information, in addition to joint information [...]". We are unsure which previous work the reviewer is referring to. Assuming this extends on the original question about [17], we want to clarify again that [17] is concerned with joint information and does not consider context-specific information. They employ the same setup to study a different problem.
Thanks for your hard work. I have read your rebuttal and I'm now sure that I have no misunderstanding about this paper. I agree that this paper has provided some novel insights. But considering that this work is not the first one investigating the problem of causal discovery with endogenous context variables and this work remains at a relatively elementary stage (both detection and latent context variables are out of your scope), I decide to maintain my borderline accept score.