iSCAN: Identifying Causal Mechanism Shifts among Nonlinear Additive Noise Models
Authors: Tianyu Chen
Kevin Bello
Bryon Aragam
Pradeep Kumar Ravikumar
Venue: NeurIPS 2023 poster (2023)
Keywords: distribution shifts
heterogeneous data
feature-shift
structural causal models
additive noise models
causality
root-cause analysis
TL;DR: Novel approach to localize causal mechanism shifts given data from multiple environments. The method can be applied to root cause analysis, fault localization, domain adaptation, explaining distribution shifts, etc.

Abstract: Structural causal models (SCMs) are widely used in various disciplines to represent causal relationships among variables in complex systems.
Unfortunately, the underlying causal structure is often unknown, and estimating it from data remains a challenging task. 
In many situations, however, the end goal is to localize the changes (shifts) in the causal mechanisms between related datasets instead of learning the full causal structure of the individual datasets. 
Some applications include root cause analysis, analyzing gene regulatory network structure changes between healthy and cancerous individuals, or explaining distribution shifts. 
This paper focuses on identifying the causal mechanism shifts in two or more related datasets over the same set of variables---*without estimating the entire DAG structure of each SCM*.
Prior work under this setting assumed linear models with Gaussian noises; instead, in this work we assume that each SCM belongs to the more general class of *nonlinear* additive noise models (ANMs).
A key technical contribution of this work is to show that the Jacobian of the score function for the *mixture distribution* allows for the identification of shifts under general non-parametric functional mechanisms.
Once the shifted variables are identified, we leverage recent work to estimate the structural differences, if any, for the shifted variables.
Experiments on synthetic and real-world data are provided to showcase the applicability of this approach.
Code implementing the proposed method is open-source and publicly available at  https://github.com/kevinsbello/iSCAN.

iSCAN: Identifying Causal Mechanism Shifts among Nonlinear Additive Noise Models
Tianyu Chen"
Kevin Bello
Bryon Aragam?
Pradeep Ravikumar
DDepartment of Statisties and Data Science, University of Texas at Austin
Booth School of Business, University of Chicago Machine Learning Department, Carnegie Mellon University
Abstract
Structural causal models (SCMs) are widely used in various disciplines repre- sent causal relationships among variables in complex systems. Unfortunately, the underlying causal structure often unknown, and estimating it from data remains challenging task. In many situations, however, the end goal is to localize the changes (shifts) in the causal mechanisms between related datasets instead of learn- ing the full causal structure of the individual datasets. Some applications include root cause analysis, analyzing gene regulatory network structure changes between healthy and cancerous individuals, explaining distribution shifts. This paper focuses on identifying the causal mechanism shifts in two or more related datasets over the same set of setrables-withone esfimating the entire DAG sfructure of each SCM. Prior work under this setting assumed linear models with Gaussian noises; instead, in this work we assume that each SCM belongs the more general class of nonlinear additive noise models (ANMs). A key technical contribution of this work is to show that the Jacobian of the score function for the mixture distribution allows for the identification of shifts under general non-parametrie functional mechanisms. Once the shifted variables are identified, we leverage recent work to estimate the structural differences, if any, for the shifted variables. Experiments on synthetic and real-world data are provided showcase the applicability of this approach. Code implementing the proposed method is open-source and publicly available at http:://ggguu...ookevinneeeeeeee
1 Introduction
Structural causal models (SCMs) are powerful models for representing causal relationships among variables in complex system [54, 58). Every SCM has an underlying graphical structure that is generally assumed to directed acyelic graph (DAG). Identifying the DAG structure of SCM is crucial since it enables reasoning about interventions [54]. Nonetheless, in most situations, scientists can only access observational or interventional data, or both, while the true underlying DAG structure remains unknown. As result, in numerous disciplines such computational biology [66, 30, 20], epidemiology [64], medicine [61, 62], and econometrics [34, 28, 18), is critically important to develop methods that estimate the entire underlying DAG structure based on available data. This task is commonly referred to as causal discovery or structure learning, for which a variety of algorithms have been proposed over the last decades.
Throughout this work, we make the assumption of causal sufficiency (i.e., non-existence of unobserved confounders). Under this condition alone, identifying the underlying DAG structure is not possible in general, and remains worst-case NP-complete [14, 16]. Indeed, prominent methods such as PC [71]
"Wor done while at the Department of Statistics at the University of Chicago. Correspondence kbello8cs cmu, edu.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

and GES [15] additionally require the arguably strong faithfulness assumption [77] to consistently estimate, large samples, the Markov equivalent class of the underlying DAG. However, these methods are not consistent in high-dimensions unless one additionally assumes sparsity or small maximum-degree of the true DAG [36, 50, 78]. Consequently, the existence of hub nodes, which is well-known feature in several networks [5, 6, 7], significantly complicates the DAG learning problem.
In many situations, however, the end goal is to detect shifts (changes) in the causal mechanisms between two (or more) related SCMs rather than recovering the entire underlying DAG structure of each SCM. For example, examining the mechanism changes in the gene regulatory network structure between healthy individuals and those with cancer may provide insights into the genetic factors contributing to the specific cancer, within biological pathways, genes could regulate various target gene groups depending on the cellular environment or the presence of particular disease conditions [32, 60]. In these examples, while the individual networks could be dense, the number of mechanism shifts could be sparse [69, 74, 55]. Finally, in root cause analysis, the goal is to identify the sources that originated observed changes in a joint distribution; this is precisely the setting we study in this work, where we model the joint distributions via SCMs, as also done in [52, 33].
In more detail, we focus on the problem of identifying mechanism shifts given datasets from two or more environments (SCMs) over the same observables. We assume that each SCM belongs to the class of additive noise models (ANMs) [29], i.e., each variable is defined as a nonlinear function over a subset of the remaining variables plus a random noise (see Section 2 for formal definitions). Importantly, we do not make any structural assumptions (e.g.. sparsity, small maximum- degree, or bounded rre-widthh) on the individual DAGs. Even though ANMs are well-known to be identifiable [29, 59], we aim to detect the local distribution changes without estimating the full structures individually. See Figure I for a toy example of what we aim to estimate. A similar setting to this problem was studied in [82, 23] albeit in the restrictive linear setting. Finally, it is worth noting that even with complete knowledge of the entire structure of each SCM, assessing changes in non-paraettric functions across different groups or environments remains very challenging problem [see for instance, 44].
Contributions. Motivated by recent developments causal structure learning of ANMs [65], we propose a two-fold algorithm that (1) Identifies shifted variables (i.e., variables for which their causal mechanism has changed across the environments); and (2) If needed, for each shifted variable, estimates the structural changes among the SCMs. More concretely, we make the following set of contributions:
To identify shifted variables (Definition 3), we prove that the variance of the diagonal elements of the Hessian matrix associated with the log-density of the mixture distribution unveils information to detect distribution shifts in the leaves of the DAGs (see Theorem 1). Due to this result, our algorithm (Algorithm 1) iteratively chooses a particular leaf variable and determines whether or not such variable is shifted. Importantly, this detection step does not rely on any structural assumptions on the individual DAGs, and can consistently detect distribution shifts for non-parametric second-order differentiability.
functtttriccuunciooono. functionals under very mild conditions such as
To identify structurally shifted edges (Definition 4), we propose nonparametric local parents recovery method (Algorithm 2) based on a recent measure of conditional dependence [3]. In addition, based on recent results in [4], we provide a theoretical justification for the asymptotic consistency of Algorithm 2 in Theorem 2. Importantly, since structural changes can only occur on shifted nodes, this second step can be conducted much more efficiently when the sparse mechanism shift hypothesis [69] holds, which posits that only a small subset of the causal model's mechanisms change.
We empirically demonstrate that our method can outperform existing methods such as DCI, which is tailored for linear models, as well as related methods for estimating unknown intervention targets such UT-IGSP [72]. See Section 5 and Appendix C for more details. Moreover, in Section 5.2, we provide experiments on an ovarian cancer dataset, thus, showcasing the applicability of our method.
2

---

X1 (N1) def N1.
X2 f2 N2) def Lanh(Xs) N2.
X3
X
Xg = (x1,, 11,NN No. def sinc( X1) + No.
del Xi X1 No. X4 (X1, NA)
X5 f3 (X3, No) def X4 sin(Xs) + Ns.
(a) Underlying SCM M
(b) Structural equations of the unknown SCM M
X1
X3
X Peree
X3
X
(ib(3r))
(e) Shifts in mechanisms between (c) Env. Ei originated by interventions on Xs and Xs- E, and &g. known intervention on Xs.
Env. &2, originated by unknown
Figure 1: Illustration of two different environments (see Definition 2) in (1c) and (1d), both orig- inated from the underlying SCM in (Ia) with structural equations given in (lb). Between the two environments, we observe change in the causal mechanisms of variables X3 and X5- red nodes in (1e). Specifically, for Xs. we observe that its functional dependence changed from X₂ in & to X3 in &g. For Xs, its structural dependence has not changed between & and &2, and only its functional changed from sinc(X,) in ε to the sigmoid function (Xx) in Ez. Finally, in (le), the red edges represent the structural changes in the mechanisms. The non-existence of edge from X X3 indicates that the structural relation between X, and X3 is invariant.
2 Preliminaries and Background
In this section we introduce notation and formally define the problem setting. We use [d] to denote the set of integers {I...., d]. Let G [[d], E) be DAG with node set [d] and of directed edges E C [d] x [d]. where any (i,j) E indicates and edge from i to j. Also let X = (X1,.... Xa) denote d-dimensional vector of random variables. An SCM M (X, Px) over d variables is generally defined as collection of d structural equations of the form:
X,- (PA,, Nz), Vj E [d].
where PA, S (X Xa) (X,) are the direct causes (or parents) of X -= ( a set of functional mechanisms f, I + R: and PN is a joint distribution" over the noise variables Nj. which we assume to be jointly independent'. Moreover, the underlying graph G' of an SCM is constructed by drawing directed edges for each X E PA, to X,. We henceforth assume this graph to be acyclic, i.e., a DAG. Finally, every SCM M defines unique distribution Ixx over the variables X [Proposition 6.3 in 58], which by the independence of the noise variables (a.k.a. the Markovian assumption), Px admits the following factorization:
P(X)
II P(X, PAj), 3=1
(2)
where P(X, PAy)) is referred as the causal mechanism of X,-
2 We will always assume the existence a density function w.r.t. the Lebesgue measure. Note that this implies that there is no hidden confounding.

---

The model above is often too general due to problems of identifiability. In this work we will consider that the noises are additive.
Definition 1 (Additive noise models (ANMs)). An additive noise model is an SCM M = PN, as in (1), where each structural assignment has the form:
Xj (PA,) + N,, Vj E [d].
Depending the assumptions on /g and Nj. the underlying DAG of ANM can be identifiable from observational data. E.g., when fj is linear and N, is Gaussian, in general one can only identify the Markov equivalence class (MEC) of the DAG, assuming faithfulness [54]. For linear models, an exception arises when assuming equal error variances [56, 78, 47], or non-Gaussian errors [70]. In addition, when fj is nonlinear on each component and three times differentiable then the DAG is also identifiable [59, 29]. Very recently, Rolland et al. [65] proved DAG identifiability when fj is nonlinear on each component and N₁ Gaussian, using information from the score's Jacobian.
2.1 Data from multiple environments
Throughout this work we assume that we observe a collection of datasets, D {
h== from H
(possibly different) environments. Each dataset X = (x**.* mn from environment contains my (possibly non-independent) samples from the joint distribution Pk i.c., X 1 We consider that each environment originates from soft interventions' [54] of an unknown underlying SCM M with DAG structure G and joint distribution (*(X) = I IP" (X, PPA)). Here PA denotes the parents (direct causes) of Xj in G Then, an environment arises from manipulations or shifts in the causal mechanisms of subset of variables, transforming from (Xx; PAA) to P(X, PAA)). Throughout, we will make the common modularity assumption of causal mechanisms [54, 38], which postulates that an intervention other mechanisms PPX;; PA:), for i & j. remain unchanged.
node Xj only changes the mechanism P(X, PA2), while all
Definition 2 (Environment). An environment En fh, ), with joint distribution Ipb. and density phd independently results from an SCM M* by intervening of causal mechanisms, that is, we can factorize the joint distribution Pa (X) as follows:
unknown subset S S [d]
Ip (X)
II P PAA) JE[d
II [I' PAA) jES
II (', PAJ).
(3)
jESN
h where PA, is a (possibly empty) subset of the underlying causal parents PAJ, iee, PA, C PAA: and, P.t PX(x; PAA)) are the invariant mechanisms.
Remark 1. In the literature [e.g., 55], it is common to find the assumption that in a soft intervention the direct causes remain invariant, i.e., PA, = PAA for all j E Sh h E [H]. In this work we consider more general setting where none, some, all of the direct causes of an intervened node removed, i.e., PA PPA.fo for all j e Shh [H].
h
We next define shifted nodes (variables).
Definition 3 (Shifted node). Given environments {Er, = fh fh,P originated from an ANM M node j is called a shifted node if there exists h, E [H] such that:
PM PPA)) * Tpt (X, PPA)).
To conclude this section, we formally define the problem setting.
Jカ=ı {X"} where X Pt consists of my (possibly non- independent) samples from the environment distribution 1o originated from underlying ANM M estimate the set shifted nodes and structural differences.
Problem setting. Given H datasets
We note that [82, 23] have study the problem setting above for H 2, assuming linear functions f/h and Gaussian noises No. In this work, we consider a more challenging setting where f/t are nonparametric functions (see Section 3 for more details).
"These types of interventions are more realistic in practice than "hard" or perfect interventions. However, note that we allow soft intervention on a variable to remove some or all of its causes, where the latter is also known as stochastic hard intervention.

---

2.2 Related Work
First we mention works most closely related to ours. The problem of learning the difference between undirected graphs has received much more attention than the directed case. E.g., [88. 46, 85, 19] develop algorithms for estimating the difference between Markov random fields and Ising models. See [87] for recent developments in this direction. In the directed setting, [82, 23] propose methods for directly estimating the difference of linear ANMs with Gaussian noise. More recently, [67] studied the setting where a dataset is generated from a mixture of SCMs, and their method is capable of detecting conditional distributions changes: however, due to the unknown membership of each sample, it is difficult to test for structural and functional changes. Moreover, in contrast to ours, all the aforementioned work on the directed setting rely on some form of faithfulness assumption.
Causal discovery from a single environment. One way to identify mechanism shifts (albeit inefficient) would be to estimate the individual DAGs for each environment and then test for structural differences across the different environments, A few classical and recent methods for learning DAGs from single dataset include: Constraint-basse algorithms such as PC and FCI [71]: in score-based methods, we have greedy approaches such as GES [16], likelihood-basee methods [56, 47, 59, 2, I, 29], and methods that test for asymmetries [70, 12]. and hybrid methods [50, 76]. Finally, note that even if we perfectly estimate each individual DAG (assuming identifiable models such as ANMs), applying these methods would only identify structural changes. That is, for variables that have the same parents across all the environments, we would require an additional step to identify distributional changes.
constrrinnd learning [89, 51, 39, 8]. Order-basee methods [75, 41, 24, 65, 48],
Testing functional changes in multiple datasets. Given the parents of a variable X,, one could leverage prior work [44, 25, 9, 26] on detecting heterogeneous functional relationships. However, we highlight some important limitations. Several methods such as [25, 9, 26] only work for one dimensional functionals and assume that the datasets share the exact same design matrix. Although [44] relaxes this assumption and extends the method to multivariate cases, the authors assume that the covariates (i.c.. PAG) are sampled from the same distribution across the environments, which is a strong assumption in our context since ancestors of X, could have experienced mechanism shifts. Finally, methods such as [53] and [11], although nonparametric, need knowledge about the parent set PAy for each variable, and they assume that PA, is same across different environments.
Causal discovery from heterogeneous data. Another well-studies problem is to learn the underlying DAG of the SCM M that originated the different environments. Under this setting, [83] provided a characterization of the I-MEC, a subset of the Markov equivalence class. [55] provided DAG- identifiability results by leveraging sparse mechanism shifts and relies on identifying such shifts, which this work aims to solve. [10] developed an estimator considering unknown intervention targets. [79] primarily focuses on linear SEM and does not adapt well to nonlinear scenarios. Also assuming linear models, [22, 21] applied ideas from linear invariant causal prediction [ICP, 57J and ICM to identify the causal DAG. [72] proposes nonparametric method that identify the intervention targets; however, this method relies on nonparametric CI tests, which can be time-consuming and sample inefficient. [49] introduced the joint causal inference (JCI) framework, which can also estimate intervention nodes. However, this method relies on an assumption that the intervention variables are fully connected, condition that is unlikely to hold in practice. [31] introduced two-stage approach that removes functional restrictions. First, they used the PC algorithm using all available data to identify the MEC. Then, the second step aims to orient the remaining edges based on novel measure of mechanism dependence. Finally, we note that a common assumption in the aforementioned methods is the knowledge of which dataset corresponds to the observational distribution; without such information, their assumptions on the type of interventions would not hold true. In contrast, our method does not require knowledge of the observational distribution.
3 Identifying Causal Mechanism Shifts via Score Matching
In this section, we propose iSCAN (identifying Shifts in Causal Additive Noise models), a method for detecting shifted nodes (Definition 3) based only on information from the Jacobian of the score of the data d
In this work, the score of a pdf p(x) means V log)(z)
5

---

(X") E R Let X be the row concatenation of all the datasets X", i... x [ - where Σ my- The pooled data X be interpreted mixture of data from the H different environments. To account for this mixture, we introduce the probability mass Wh. which represents the probability that an observation belongs to environment h, i.c., Σ Wh 1. Let Q(X) denote the distribution of the mixture data with density function g(z), i.e. g(x) Σ A== (x)
In the sequel, use s" (x) V log log (x) to denote the score function of the joint distribution of environment h with density p" (ar). Also, we let s(x) V Vlog (ax) to denote the score function of
the mixture distribution with density ((x). We will make the following assumptions on f and Nţ.
[d], the functional mechanisms (PA)) are assumed to be Assumption A. For all h E [H].j E
non-linear in every component.
Assumption B. For all E [d], E [H). the pdf of the real-valued noise Nţ denoted by p'x, satisfies log PN, (n²)) =c where d) is non-zero constant. Moreover, B[N
82 (ong))
For an ANM, Rolland et al. [65] showed that under Assumption A and assuming zero-mean Gaussian noises (which satisfies Assumption B), the diagonal of the Jacobian of the score function reveals the leaves of the underlying DAG. We next instantiate their result in our context.
Proposition (Lemma in [65, 68]). For an environment En with underlying DAG Gh and pdf (xr), let (xr) v log (x)) be the associated score function. Then, under Assumptions A and B, for all j E [d], we have:
(x) 82,
Node j is leaf in G = <Vaxx
0.
Motivated by the ideas of leaf-identiffabilityy from the score's ' Jacobian single ANM, we next show that the score's 's Jucobian of the mixture distribution can help reveal mechanism shifts among the different environments.
Theorem 1. For all h e [H] let G and (2x) denote the undertying DAG structure and pdf of environment En. respectively, and let q(x) be the pdf of the mixture distribution of the H environments such that q(z) = EH whpp" (x). Also, let s(x) = V log q(z) be the associated score function. Then, under Assumptions A, and B, we have:
(i) i is leaf in all DAGs G then j is shifted node if and only if Var (s,(X) dx,
> 0.
(ii) j is not leaf in at least one DAG Gh, then Var
(X) 0.
Theorem I along with Proposition 1 suggests a way to identify shifted nodes. Namely, to use Proposition to identify common leaf, and then use Theorem test if such leaf is a shifted node or not. We then proceed to remove the leaf and repeat the process, See Algorithm 1. Note that due the fact that each environment is result of an intervention (Definition 2) on an underlying ANM M it follows that the leaves in G* will remain leaves in each DAG Gh.
Algorithm iSCAN ISCAAN-Leentffying Shifts in Causal Additive Noise models.
Input: Datasets X x.....,X'
Output: Shifted variables set 3, and topological sort f.
Initialize $ 5= O) = (). f1, [l.....d] .d]
2: Set X [(X2)
(X TJ RRdd
3: while N 7 2 do
4: vh E H,, Var + Varx [diag( V (olgg"
5:
Var + Varx [daag( 2 log q(x))
L Nhe(m {i Var 0,j E [d]
6:
Identify leaves. Identify shifted nodes.
$ - $ U {i Var, + 0.j E L}
7:
8:
N+N- N (L)
9:
Vl L, remove the 1-th column of h Wh E [M], and X.
10:
# t))

---

See Appendix A for a practical implementation Alg. 1. Finally, note that Alg. also estimates a valid topological sort for the different environments by leveraging Proposition 1.
Remark 2.
3.1 Score's Jacobian estimation
Since the procedure to estimate Var,l [w) is similar for estimating Var in this section discuss the estimation for Var, de, \ which involves computing the diagonal
8s (x) dd,,
for each E [H],
of the Hessian of log q(x). To estimate this quantity, we adopt a similar approach to the method in
[45, 65]. First, we estimate the first-order derivative of log g(z) by Stein's ' identity [73]:
Bg V log (=)) + Vh(x)) 0,
where h IR + Rd is any test function such that limg- (x)) = 0. Once we have estimated V log (<x), we can proceed estimate the Hessian's 's diagonal by using second-order Stein's identity:
Bgg (x))(d)daag( V log (2))) [] = VaggV h(x) h( V log q(x) V log q(x) (5)
Using (4) and (5), we can estimate the Hessian's 's diagonal at each data point. Thus allowing us to obtain an estimate of Var, duy
[x,(z)) See Appendix A.1 for additional details.
Remark 3 (Consistency of Algorithm 1). The estimators (6) and eq. (7)(7), given in Appendix
A. correspond to Monte Carlo estimators using eq.(4) and (5), respectively, then the error of the estimators tend zero as the number of samples goes infinity. See for instance the discussion in Section 3.1 in [45]. We empirically explore the consistency of Algorithm 1 in Figure 2.
Remark 4 (Computational Complexity). Since adopt the kernel-basee estimator, SCORE, from [65]. The computational complexity for the estimation of the score's Jacobian in a single environment O (ddm)) In Algorithm I, computation is dominated by the SCORE function applied to the pooled data X E Rmaxd Therefore, the overall complexity of Algorithm is Odm2)) See Figure 2.
Gaussian
Gaussian
100
Runtime (seconds)
30
F1 Score (%)
80-
20-
P-28
10
0
1000
1500
2000
20
40
m (number of samples)
(number of nodes)
m
Figure 2: (Left) Fl score of the output of Alg. w.r.t. to the true set of shifted nodes. For different number of nodes, we observe how iSCAN recovers the true set of shifted nodes the number of samples increases, thus empirically showing its consistency. (Right) Runtime number of nodes for different number of samples. We corroborate the linear dependence of the time complexity on d.
4 On Identifying Structural Differences
After estimating the set of shifted nodes Š through Algorithm 1, it is of high interest to predict which causal relations between a shifted node and its parents have undergone changes across the environments. The meaning of change in causal relationship can vary based on the context and the estimation objective. This section primarily centers on structural changes, elaborated further below, while additional discussion about other types of changes is available in Appendix D.
Definition (Structurally shifted edge). For a given shifted node X,, an edge X; + X5 is called a structurally shifted edge 3,, h E [H] such that X; PA and X; $ PAY
In other words, a structurally shifted edge is an edge that exists in one environment but not in another, indicating a change in the underlying structure of the causal mechanism. To detect the structurally shifted edges, we will estimate the parents of each shifted node in S for all environments En.
7

---

Remark 5. Note that under the sparse mechanism shift hypothesis [69], i.e., |S| & d. estimating the parents of each shifted node much more efficient than estimating the entire individual structures.
Kernel regression and variable selection. A potential strategy to estimate structurally shifted edges involves employing the estimated topoloical order obtained from Algorithm If this estimated topological order remains valid across all environments, it can serve as guide for the nonparametric variable selection process to identify the parents of shifted node Xj- Specifically, we can regress the shifted node Xj on its predecessors (re(X,) and proceed with a nonparametric
variable selection procedure. Here Pro(X;) consists of the set of nodes that appear before X; in the estimated topological order ît To achieve that, there exist various methods under the hypothesis testing framework [42, 17, 63], and bandwidth selection procedures [40]. These methods offer consistency guarantees, but their time complexity might be problematic. Kernel regression, for example, has time complexity of O(m²), and requires an additional bandwidth selection procedure, usually with time complexity of ( Consequently, it becomes imperative to find a more efficient method for identifying parents locally.
Feature ordering by conditional independence (FOCI), An alternative efficient approach for identifying the parents is to leverage the feature ordering method based on conditional independence proposed by Azadkia and Chatterjee [3]. This method provides a measure of conditional dependency between variables with a time complexity of O(m log m). By applying this method, we can perform fast variable selection in nonparametric setting, See Algorithm 4 in Appendix A.3.
Theorem 2 (Consistency of Algorithm 4). Under Assumption C, given in Appendix B.2, if the estimated topological order ît output from Algorithmm I is valid for all environments. then the output
h
PA of Algorithm 4 is equal to the true parents PA of node X, with high probability, for all h [H].
Motivated by Theorem 2, we next present Algorithm 2, a procedure to estimate the structurally shifted edges. Given the consistency of Alg. and Alg. 2, it follows that combining both algorithms will correctly estimate the true set of shifted nodes and structural shifted edges, asymptotically.
Algorithm 2 Identifying structurally shifted edges
Input: Data (X*) JN((M topological order म, shifted nodes $
Output: Structurally shifted edges set E
I: Initialize =
2: for X, in S do
3: for h in [H] do
4: Estimate PA, h
from Alg. 4 (FOCI) with input (Pre(X))), X/}} -
5: if EX, h ht' such that Xa PA,. X2 & PA, then
È + EXE,
6:
5 Experiments
We conducted comprehensive evaluation of our algorithms. Section 5.1 focuses on assessing the performance of iSCAN (Alg. 1) for identifying shifted variables. In Section 5.2, we apply iSCAN for identifying shifted nodes along with FOCI (Alg. 2) for estimating structural changes, on apoptosis data. Also, in App. C, we provide additional experiments including: (i) Localizing shifted nodes without structural changes (App. C.I), and where the functionals are sampled from Gaussian processes (App. C.1.1); (ii) Localizing shifted nodes and estimating structural changes when the underlying graphs are different; and (iii) Evaluating iSCAN using the elbow method for selecting shifted nodes (see App. C.3 and Remark 6). Code is publicly available https///ghuub.coo/eevissseelllle
5.1 Synthetic experiments shifted nodes
Graph models. We generated random graphs using the Frdōs-Rény (ER) and scale free (SF) models. For given number of variables d, ERK and SFk indicate average number of edges equal to kd.
8

---

Data generation. We first sampled DAG, G', of d nodes according to either the ER or SF model for env. &j. For env. Ez. we initialized its DAG structure from env. & and produced structural changes by randomly selecting 0.2 d nodes from the non-root nodes. This set of selected nodes S, with cardinality |S] 0.2d, correspond to the set of "shifted nodes". In env. &2, for each shifted node X, E S, we uniformly at random deleted at most three of its incoming edges, and use D, to denote the parents whose edges to Xj were deleted; thus, the DAG G² is subgraph of G Then, in E,, each X, was defined as follows:
Xy
Σ iEPAJADJ
sin(X?}) +
Σ IED
cos 3X,)) + N
In &2, each Xj was defined as follows:
x,= Σ tEPA}
sin X?) + N₁
Experiment details. For each simulation, we generated 500 data points per environment, i.c., mi 500, my 500 and m 1000. The noise variances were set to We conducted 30 simulations for each combination of graph type (ER or SF), noise type (Gaussian, Gumbel, and Laplace), and number of nodes (d E {10, 20, 30, 503). The running time was recorded by executing the experiments on an Intel Xeon Gold 6248R Processor with 8 cores. For our method, we used (see Alg. 3).
? 0.05 for eq-(6) and (7) and a threshold
Evaluation. We compared the performance of iSCAN against several baselines, which include: DCI [82], the approach by [11], CITE [79], KCD [53], SCORE [65], and UT-IGSP [72]. Figure illustrates the results for ER4 and SF4 graphs. We note that iSCAN consistently outperforms other baselines in terms of F1 score across all scenarios. Importantly, note how the performance some baselines, like DCI, CITE, Budhathoki's and SCORE, degrades faster for graphs with hub nodes, contrast, iSCAN performs similarly, it not dependent on structural
property of SF graphs. assumptions on the individual DAGs. Additionally, it is worth noting that our method exhibits faster computational time than KCD, Budhathoki's. and SCORE, particularly for larger numbers of nodes.
In Appendix C.I, we provide experiments on sparser graphs such as ER2/SF2, and denser graphs such FR6/SF6. We also include Precision and Recall in all plots in the supplement.
Gaussion
Gumbol
Laplace
Goussian
Gumbel
Leslace
400
F1 Scor 1%6
Runtime (5)
4000
3000
60
2000
PJ5
40
30
50
(number of nodes)
(number of nodes)
Methods
Budhathoki
CITE
DC1
KCD
ISCAN
SCORE
UT-IGSP
Figure 3: Experiments on ER4 and SF4 graphs. See the experiment details above. The points indicate the average values obtained from these simulations. The error bars depict the standard errors. Our method iSCAN (light blue) consistently outperformed baseline methods in terms of F1 score.
5.2 Experiments on apoptosis data
We conducted an analysis on an ovarian cancer dataset using iSCAN (Algorithm 1) to identify shifted nodes and Algorithm 2 to detect structurally shifted edges (SSEs). This dataset had previously been analyzed using the DPM method [88] in the undirected setting, and the DCI method [82] in the linear setting. By applying our method, we were able to identify the shifted nodes and SSEs in the datasct (see Figure 4a). Our analysis revealed the identification of two hub nodes in the apoptosis

---

pathway: BIRC3, and PRKAR2B. The identification of BIRC3 as a hub node was consistent with the results obtained by the DPM and DCI methods. Additionally, our analysis also identified PRKAR2B as hub node, which was consistent with the result obtained by the DCI method. Indeed, BIRC3, in addition to its role in inhibiting TRAIL-induced apoptosis, has been investigated as a potential therapeutic target in cancer treatment including ovarian cancer[35, 81]; whereas PRKAR2B has been identified as an important factor in the progression of ovarian cancer cells. The latter serves as key regulatory unit involved in the growth and development of cancer cells [84, 13].
IIIII
FAS
II.RI FAS
CSF2RB
PRKAR2B
CSF2RB
PRKAR2B
PIK3R1
PIK3R1
BIRC3
BIRC3
TP53
AIFMI
TP53
AIFMI
TNESFIO ENDOG
TNFSFI0 ENDOG
(a) The red nodes the shifted nodes estimated by iSCAN (Alg. 1). The edges are the structurally shifted edges estimated by FOCI (Alg. 2).
(b) Undirected difference network estimated by DPM
[88]. The red nodes indicate hub nodes, however, it not clear which node mechanisms have changed.
Figure 4: Results on apoptosis data.
Conclusion
In this work, we showed a novel connection between score matching and identifying causal mech- anism shifts among related heterogeneous datasets. This finding opens up a new and promising application for score function estimation techniques.
Our proposed technique consists of three modules. The first module evaluates the Jacobian of the score under the individual distributions and the mixture distribution. The second module identifies shifted features (variables) using the estimated Jacobians, allowing us to pinpoint the nodes that have undergone mechanism shift. Finally, the third module aims to estimate structurally shifted edges, a.k.a. the difference DAG, by leveraging the information from the identified shifted nodes and the estimated topological order. It is important to note that our identifiability result in Theorem is agnostic to the choice of the score estimator.
The strength of our result lies in its capability to recover the difference DAG in non-linear Additive Noise Models (ANMs) without making any assumptions about the parametric form of the functions or statistical independencies. This makes our method applicable in a wide range of scenarios where non-linear relationships and shifts in mechanisms are present.
6.1 Limitations and future work
While our work demonstrates the applicability of score matching in identifying causal mechanism shifts in the context of nonlinear ANMs, there several limitations and areas for future exploration:
Extension to other families of SCMs: Currently, our method primarily focused on ANMs where the noise distribution satisfies Assumption B, e.g., Gaussian distributions. It would be valuable to investigate the application of score matching in identifying causal mechanism shifts in other types of SCMs. Recent literature, such [48], has extended score matching to additive Models with arbitrary noise for finding the topological order. Expanding our method to accommodate different noise models would enhance its applicability to wider range of real-world scenarios.
Convergence rate analysis: Although the score matching estimator is asymptotically consistent, the convergence rate remains unknown in general. Understanding the convergence properties of the estimator is crucial for determining the sample efficiency and estimating the required number of samples to control the estimation error within desired threshold. Further theoretical developments, such as [37], on score matching estimators would provide valuable insights into the performance and sample requirements of iSCAN.
10

---

Acknowledgments and Disclosure of Funding
K.B. was supported by NSF under Grant 2127309 to the Computing Research Association for the CIFellows 2021 Project. B.A. was supported by NSF IIS-1956330, NIH R01GM140467, and the Robert H. Topel Faculty Research Fund at the University of Chicago Booth School of Business. P.R. was supported by ONR via N00014181 181661, and NSF via IIS- 16, NIS-195533,, IIS AS-22111077 We are also grateful for the support of the University of Chicago Research Computing Center for assistance with the calculations carried out in this work.
References
[1] Aragam, B., Amini, A. and Zhou, Q. [2019), 'Globally optimal score-based learning of directed acyclic graphs in high-dimensions' Advances in Neural Information Processing Systems 32.
[2] Aragam, B. and Zhou, Q. [2015], 'Concave penalized estimation of sparse Gaussian Bayesian networks', The Journal of Machine Learning Research 16(1), 2273-2328
[3] Azadkia, M. and Chatterjee, S. [2021], 'A simple measure of conditional dependence', The Annals of Statistics (9),), 3070-3102.
[4] Azadkia, M., Taeb, and Bühlmann, [2021], 'A fast non-parametric approach for local causal structure learning', arXiv e-prints pp. arXiv-2111.
[5] Barabási, A.-L. and Albert, R. [1999], 'Emergence of scaling in random networks', science 286(5439), 509-512.
[6] Barabási, A.,L., Gulbahce, N. and Loscalzo, J. [2011], 'Network medicine: a network-based approach to human disease', Nature reviews genetics 12(1), 56-68.
[7] Barabasi, A.-L.. and Oltvai, Z.. N. [2004], 'Network biology: understanding the cell's functional organization', Nature reviews genetics 5(2), 101-113.
[8] Bello, K., Aragam, B. and Ravikumar, P. [2022], 'Dagma: Learning dags via m-matrices and a segedettrmmnant acyclicity characterization' arXiv preprint arXiv:2209.0800 08037.
[9] Boente, G. and
C. [2022], 'Robust testing to compare regression curves',
arXiv preprint arXiv:2205. /2065
[10] Brouillard, P., Lachapelle, S., Lacoste, A., Lacoste-Julien, S. and Drouin, A. [2020], 'Differen- tiable causal discovery from interventional data', Advances in Neural Information Processing Systems 33, 21865-21877.
[11] Budhathoki, K., Janzing, D., Bloebaum, P. and Ng. H. [2021], Why did the distribution change?, in 'International Conference on Artificial Intelligence and Statistics', PMLR, pp. 1666-1674.
[12] Bühlmann, P., Peters, J. and Ernest, J. [2014), 'Cam: Causal additive models, high-dimensional order search and penalized regression' The Annaly of Statistics 42(6), 2526-2556.
[13] Chiaradonna. Balestrieri, C., Gaglio, D. and Vanoni, M. [2008], 'Ras and pka pathways can- cer: new insight from transcriptional analysis'. Frontiers in Bioscience- 5278.
13(14), 5257-
[14] Chickering, M. [1996], Learning bayesian networks is np-complete, 'Learning from data', Springer, pp. 121-130.
[15] Chickering, D. M. [2003], 'Optimal structure identification with greedy search', JMLR
3, 507-554.
[16] Chickering, D. M., Heckerman, D. and Meek, C. [2004), 'Large-sample learning of Bayesian networks is NP-hard', Journal of Machine Learning Research 5, 1287-1330.
[17] Delgado, M. A. and Manteiga, W. G. [2001]. 'Significance testing in nonparametric regression based on the bootstrap', The Annals of Statistics 29(5), 1469-1507.
[18] Demiralp, S. and Hoover, K. D. [2003], 'Searching for the causal structure of a vector autore- gression', Oxford Bulletin of Economics and statistics 65, 745-767.
[19] Fazayeli, and Banerjee, A. [2016], Generalized Direct Change Estimation in Ising Model Structure, in 'International Conference Machine Learning'.

---

[20] Friedman, N., Linial, M., Nachman, I. and Pe'er, D. [2000], Using bayesian networks to analyze expression data, in Proceedings of the fourth annual international conference on Computational molecular biology', pp. 127-135.
[21] Ghassami, A., Kiyavash, N., Huang, B. and Zhang, K. [2018), 'Multi-domain causal structure learning in linear systems" Advances neural information processing systems 31.
[22] Ghassami, A., Salchkaleybar, S., Kiyavash, N. and Zhang. K. [2017), 'Learning causal structures
using regression invariance', Advances Neural Information Processing Systems 30.
[23] Ghoshal, A., Bello, K. and Honorio, [2019], 'Direct learning with guarantees of the difference dag between structural equation models', arXiv preprint arXiv: 1906.
[24] Ghoshal, A. and Honorio. J. [2018], Learning linear structural equation models in polynomial time and sample complexity. in Proceedings of the Twenty- -Fistt International Conference on Artificial Intelligence and Statistics', Vol. 84 of Proceedings of Machine Learning Research, PMLR, pp. 1466-1475
[25] Hall, P. and Hart, J. D. [1990], 'Bootstrap test for difference between means in nonparametric regression', Journal of the American Statistical Association 85(412), 1039-1049.
[26] Hardle, W. and Marron, S. [1990], 'Semiparametric comparison of regression curves', The Annals of Statistics pp. 63-89.
[27] Hastic, T. J. [2017]. Generalized additive models, in 'Statistical models in S', Routledge, pp. 249-307.
[28] Hoover, K. D., Demiralp, S. and Perez, S. J. [2009], 'Empirical identification of the vector autoregression: The causes and effects of us m2', The methodology and practice of econometrics: Festschrift in honour of David F. Hendry pp. 37-58.
[29] Hoyer, P., Janzing, D., Mooij, J. M., Peters, J. and Schölkopf, B. [2008], 'Nonlinear causal discovery with additive noise models', Advances in neural information processing systems 21.
[30] Hu, P., Jiao, R., Jin, L. and Xiong, M. 8), 'Application of causal inference to genomic analysis: advances in methodology'. Frontiers in Genetics 9, 238.
[31] Huang, B., Zhang, K., Zhang. J., Ramsey, J., Sanchez-Roment R., Glymour, C. and Schölkopf.
B. [2020], 'Causal discovery from beccrogencouu/hoostttiinayy data', The Journal of Machine Learning Research 21(1). 3482-3534.
[32] Hudson, N. J., Reverter, A. and Dalrymple, B. P. [2009], 'A differential wiring analysis of ex- pression data correctly identifies the gene containing the causal mutation', PL.S computational biology 5(5), e1000382.
[33] Ikram, A., Chakraborty, S., Mitra, S., Saini, S., Bagchi, S. and Kocaoglu, M. [2022), 'Root cause analysis of failures in microservices through causal discovery', Advances in Neural Information Processing Systems 35, 31188-31170.
[34] Imbens. G. W. [2020], 'Potential outcome and directed acyclic graph approaches to causality: Relevance for empirical practice in economies', Journal of Economic Literature 58(4), 1129- 1179.
[35] Johnstone, R. W., Frew, J. and Smyth, M. [2008], "The trail apoptotic pathway in cancer onset, progression and therapy', Nature Reviews Cancer 8(10), 782-798.
[36] Kalisch, M. and Bühlman, P. [2007], 'Estimating high-dimensional directed acyelic graphs with the pe-algorithm.', Journal of Machine Learning Research 8(3).
[37] Koehler, F., Hcckett, A. and Risteski, A. [2022], 'Statistical efficiency of score matching: The view from isoperimetry', arXiv preprint arXiv:2210.001 10.0072.
[38] Koller, D. and Friedman, N. [2009], Probabilistic graphical models: principles and techniques,
MIT press.
[39] Lachapelle, S., Brouillard, P., Delcu, Т. and Lacoste-Julien, S. [2019]. 'Gradient-based neural dag learning'. arXiv preprint arXiv: 99660002266
[40] Lafferty, J. and Wasserman, I. [2008), 'Rodeo: Sparse, greedy nonparametric regression'.
[41] Larranaga, P., Kuijpers, C. M., Murga, R. H. and Yurramendi, Y. [1996], 'Learning bayesian network structures by searching for the best ordering with genetic algorithms', IEFE transactions on systems, man, and cybernetics-part A: systems and humans 26(4), 487-493.
12

---

[42] Lavergne, P. and Vuong, Q. [2000]. Nonparametrie significance testing", Econometric Theory 16(4). 576-601.
[43] Li, C., Shen, X. and Pan, W. [2023], 'Nonlinear causal discovery with confounders' Jourrnal of the American Statistical Association pp. 1-10.
[44] Li, X.. Jiang, B. and Liu, J. S. [2021], "Krnel-based partial permutation test for detecting heterogeneous functional relationship' Joounnal of the American Statistical Association pp. 1-19.
[45] Li, Y. and Turner, R. [2017], 'Gradient estimators for implicit models', arXiv preprint arXiv:1705.0770
[46] Liu, S., Suzuki, T., Relator, R., Sese, J.. Sugiyama, M., Fukumizu, K. et al. [2017]. 'Support consistency of direct sparse-change learning markov networks', The Annals of Statistics.
[47] Loh, P.-I.. and Buhlmann, P. [2014], Highh-Diensiooal Learning of Linear Causal Networks via Inverse Covariance Estimation' Journal of Machine Learning Research.
[48] Montagna, F., Noceti, N., Rosasco, L.., Zhang, K. and Locatello, F. [2023], 'Causal discovery with score matching on additive models with arbitrary noise', arxiv:2304.03220
[49] Mooij. J. M., Magliacane, S. and Claassen, T. [2020], 'Joint causal inference from multiple contexts', The Journal of Machine Learning Research 21(1), 3929
[50] Nandy, P., Hauser, A. and Maathuis, M. H. [2018], High-dimensional consistency in score- based and hybrid structure learning', The Annals of Statistics 46(6A), 3151-3183.
[51] Ng. I., Ghassami, A. and Zhang, K. [2020], 'On the role of sparsity and dag constraints for learning linear dags', Advances in Neural Information Processing Systems 33, 17943-17954.
[52] Paleyes, A., Guo, S., Scholkopf, B. and Lawrence, N. D. [2023], Dataflow graphs as complete causal graphs, in '223 IEEE/ACM 2nd International Conference on Al Engineering-Softve Engineering for AI (CAIN), IEEE, pp. 7-12.
[53] Park, J., Shalit, U., Schölkopf, B. and Muandet, K. [2021]. Conditional distributional treatment effect with kernel conditional mean embeddings and u-statistic regression, in 'International Conference on Machine Learning.. PMLR. pp. 8401-8412.
[54] Pearl, J. [2009], CAUSALITY: Models, Reasoning, and Inference, 2nd edn, Cambridge University Press.
[55] Perry, R., Von Kügelgen, J. and Schölkopf. B. [2022]. 'Causal discovery in heterogencous environments under the sparse mechanism shift hypothesis', arXiv preprint arXiv:2206.028 02013.
[56] Peters, J. and Bühlmann. P. [2014], 'Identifiability of gaussian structural equation models with equal eτror variances', Biometrika 101(1), 219-228.
[57] Peters. J., Bühlmann, P. and Meinshausen, N. [2016]. 'Causal inference by using invariant
prediction: identification and contidence intervals', Journal of the Royal Statistical Society. Series B (Statistical Methodology) pp. 947-1012.
[58] Peters, J. Janzing, D. and Schölkopf. B. [2017], Elements of causal inference: foundations and learning algorithms, The MIT Press.
[59] Peters, J., Mooij, J. M., Janzing, D. and Schölkopf, B. [2014], 'Causal discovery with continuous additive noise models'.
[60] Pimanda, J. E., Ottersbach, K., Knezevic, K., Kinston, S., Chan, W. Y., Wilson, N. K., Landry, J.-R., Wood, A. D., Kolb-Kokocinski, A., Green, A. R. et al. [2007], 'Gata2, flil, and sel form a recursively wired gene-regulatory circuit during early hematopoietic development', Proceedings of the National Academy of Sciences 104(45), 17692-
[61] Plis, S. M., Calhoun, V. D., Weisend, M. P., Eichele, T. and Lane, T. [2010]. 'Meg and fmri fusion for non-linear estimation of neural and bold signal changes', Frontiers in neuroinformatics 4, 114.
[62] Plis, S. M., Weisend, M. P., Damaraju, E., Eichele, T., Mayer, A., Clark, P., Lane, T. and Calhoun, V. D. [2011], 'Effective connectivity analysis of fmri and meg data collected under identical paradigms', Computers in biology and medicine 41(12), 1156-1165.
[63] Racine, J. [1997], "Consistent significance testing for nonparametric regression', Journal of Business & Economic Statistics 15(3), 369-378.
13

---

[64] Robins, J. M., Hernan, M. A. and Brumback, B. [2000]. 'Marginal structural models and causal inference in epidemiology'. Epidemiology pp. 550-560.
[65] Rolland, P., Cevher, V., Kleindessner, M., Russell, C., Janzing, D., Schölkopf, B. and Locatello,
F. [2022], Score matching enables causal discovery of nonlinear additive noise models, in "International Conference on Machine Learning', PMI.R, pp. 18741-18753.
[66] Sachs, K., Perez. O., Pe'er, D., Lauffenburger, D. and Nolan, G. P. [2005], 'Causal protein- signaling networks derived from multiparameter single-cell data', Science 308(5721), 523-529.
[67] Saced, B.. Panigrahi, S. and Uhler, C. [2020], Causal structure discovery from distributions arising from mixtures of dags, in 'International Conference on Machine Learning', PMLR, pp. $336-8345.
[68] Sanchez, P., Liu, X., O'Neil, A. Q. and Tsaftaris, S. A. [2022], 'Diffusion models for causal discovery via topological ordering', arXiv preprint arXiv:2210.06201.
[69] Schölkopf, B., Locatello, F., Bauer, S., Ke, N. R., Kalchbrenner, N., Goyal, A. and Bengio, [2021], 'Toward causal representation learning', Proceedings of the IEEE 109(5), 612-634.
[70] Shimizu, S., Hoyer, P. O., Hyvärinen, A., Kerminen, A. and Jordan, M. [2006], 'A linear non-gaussian acyclic model for causal discovery.' Journal of Machine Learning Research.
[71] Spirtes, P., Glymour, C. N., Scheines, R. and Heckerman, D. [2000], Causation, prediction, and search, MIT press.
[72] Squires, C., Wang, Y. and Uhler, C. [2020], Permutation-based causal structure learning with unknown intervention targets, 'Conference on Uncertainty in Artificial Intelligence', PMLR, pp. 1039-1048.
[73] Stein, C. [1972], A bound for the error in the normal approximation to the distribution of a sum of dependent random variables, in 'Proc. Sixth Berkeley Symp. Math. Stat. Prob.', pp. 583-602.
[74] Tanay, A., Regev, A. and Shamir, R. [2005], 'Conservation and evolvability in regulatory networks: the evolution of ribosomal regulation in yeast', Proceedings of the National Academy of Sciences 102(20), 7203-7208.
[75] Teyssier, M. and Koller, D. [2012], 'Ordering-based search: A simple and effective algorithm for learning bayesian networks', arXiv preprint a2X7 1429
[76] Tsamardinos, I., Brown, I.. F. and Aliferis, C. F. [2006]. 'The max-min hill-climbing Bayesian network structure learning algorithm', Machine Learning 65(1), 31-78.
[77] Uhler, C., Raskutti, G., Bühlmann, P. and Yu, B. [2013]. 'Geometry of the faithfulness assump- tion in causal inference', The Annals of Statistics pp. 436-463.
[78] Van de Geer, S. and Bühlmann, P. [2013), 'o-pennliized maximum likelihood for sparse directed acyclic graphs', The Annals of Statistics 41(2), 536-567.
[79] Varici, B., Shanmugam, K., Sattigeri, P. and Tajer, A. [2021]. 'Scalable intervention target estimation in linear models', Advances in Neural Information Processing Systems 34, 1494- 1505.
[80] Voorman, A., Shojaie, A. and Witten, D. [2014], 'Graph estimation with joint additive models', Biometrika 101(1), 85-
[81] Vucic, D. and Fairbrother, W. J. [2007], 'The inhibitor of apoptosis proteins as therapeutic targets in cancer', Clinical cancer research 13(20), 5995-6000.
[82] Wang, Y., Squires, C., Belyaeva, A. and Uhler, C. [2018], "Direct estimation of differences in causal graphs', Advances in neural information processing systems 31.
[83] Yang, K., Katcoff, A. and Uhler, C. [2018), Characterizing and learning equivalence classes of causal dags under interventions, in "International Conference on Machine Learning'. PMLR, pp- 5551-5
[84] Yoon, H., Jang, H., Kim, E.-Y., Moon, S., Lee, S., Cho, M., Cho, H. J., Ko, J., Chang, E. M., Lee, K.-A. et al. [2018]. 'Knockdown of pr results in the failure of oocyte maturation', Cellular Physiology and Biochemistry 45(5), 2009-2020.
[85] Yuan, H., Xi, R., Chen, C. and Deng, M. [2017], 'Differential network analysis via lasso penalized D-trace loss', Biometrika.
14

---

[86] Zhang, K., Peters, J. Janzing, D. and Schölkopf. B. [2012], 'Kemel-based conditional indepen- dence test and application in causal discovery'.. arXiv preprint arXiv:1202.3775.
[87] Zhao, B., Wang, Y. S. and Kolar, M. [2020], 'Fudge: Functional differential graph estimation with fully and discretely observed curves', arXiv preprini arXiv:2003. 05542.
[88] Zhao, S. D.. Cai, T. T. and Li, H. [2014], 'Direct estimation of differential networks', Biometrika 101(2), 253-268.
[89] Zheng, X., Aragam, B., Ravikumar, P. K. and Xing, E. P. [2018), 'Dags with no tears: Continu- ous optimization for structure learning', NeurIPS.
15

---

SUPPLEMENTARY MATERIAL iSCAN: Identifying Causal Mechanism Shifts among Nonlinear Additive Noise Models
A Practical Implementation
In this section, we present a more practical version of Alg. 1 that considers estimation errors, see Alg.
3. First, provide more details of the score's Jacobian estimation.
A.1 Practical Version of SCORE
Let = _='"...., xn a dataset m possibly non-independent but identically distributed samples. From Li and Turner [45], we next present the estimator for the point-wise first-order partial derivative, corresponding to eq.(
G - + -1) (V,)
(6)
where H = (h(="),...,,hl
E gd gd'xm Xm
Vh = 1 Σn Vh((z), K = H H, Kij
h(z'))((>') k(x', 2 = h(x') h(*)), (V. K) = Vh. ( K)
m V V,n((z' 2²)) and 17 2 0 is
regularization parameter. Here & is used to approximate G (V log log" Rmxdd
From [65], present the estimator for the diagonal elements of the score's Jacobian the sample points, i.e. J = (diag( V log ((.")))))..(llll (diag( (V2 )o))("")))) E Raxdd the estimator of is:
j = -diag (GGT) + +K + -1) (V2
ding K)
where H (h(x*), h(x=))) E Rdd'xm, xm, Valag h 1 Ei=1
V2 ding h(**), (Vaing
ding h(x))s
m
K = H H. Kij = 1z(z x²) h(r') h(r²), h(r')
(V K) = mH diag h.
(V2 Viing ding K)ij -==
)a(atta*)) and rj > 0 is regularization parameter. (Ou))²
In the sequel, we use SCORE(X) to denote the procedure to compute the sample variance for the estimator of the diagonal of the score's Jacobian via eq.
A.2 Practical Version of Algorithm 1
Let Var h vector rank", which represents the index of each element in Var after non-decreasing sorting. - For example, if Var (5.2,3.1,44 4.5, 1.6), then rank = (3, 1, 2, 0). Furthermore, we define d-dimensional vector rank as the element-wise summation of rank over all h E [H]. In other words, rank is calculated as rank = Σne(H] r
be d-dimensional vector, where d is the number of nodes. We introduce d-dimensional ---h
Recall that in Section 3. we remarked that we leverage the SCORE approach from Rolland et al. [65] for estimating diag( 2 log p(a:)) at each data point. Recall also that identifiability result (Theorem 1) depends on determining whether a leaf node has variance Var, deæ(x)) = 0. In practice, it is unrealistic to simply test for the equality Var, 0 since Var,, carries out errors due to finite samples. Instead, we define the following statistic for each estimated leaf node L (Line 10 in Algorithm 3);
Varz.
stats,, min, Varl +
(8)

---

Algorithm 3 Practical version of Algorithm
Input: Datasets X ..... X threshold
Output: Shifted variables set S, and topological sort
1: Initialize 5 0, = ().W= = (I,...,d}
2: stats <- (0. (0......) 0) E Rd
3: Set X [(X -.... (X RRaxd
4: while N 7 0 do
5: vh E [H], Var + (C
D Estimate Varx* [diag] V (x)))
6: vh [H], rank arg sort (Var)
7: rank + ΣnE(H) rank"
8: DE arg min, ranky
Estimate a leaf node
9: Var + SCORE(X)
> Estimate Varx [diag( (oggg(logggg (z)))
stats Varz
minA Varz h
11:
N (L)
12: Remove the L-th column of X", vh [H], and X.
13:
* (I.
S {i stats, > Vj [d}
The intuition behind this ratio is that if the leaf node is shifted node then we can expect ming Var to be large since Varz 0 (by Theorem 1), and Varž = 0 (by Proposition 1). On the other hand, if the leaf node L is not shifted node then we can expect ming Var to be small. This is due to the fact that, given a consistent estimator, Var, would converge towards (by Theorem 1) at a faster rate than Varž since we utilize larger amount of data for estimating Var,. Finally, e in the denominator very small value, e.g. 10-9. and acts safeguard against encountering division by zero".
Vary
Vary
Then, given the statistic eq.(8), we can set threshold t and define the set of shifted nodes by all the nodes j such that stats; t (Line 14 in Algorithm 3).
Remark 6 (The elbow strategy). Alternatively, we can employ an adaptive approach to identify the set of shifted nodes by sorting stats non-increasing order, and for the "elbow" point. For example, Figure 5 illustrates the variance ratio in (8) for each node sorted in non- non-inoneenn con-increasing order. In this case, node index corresponds to the elbow point, allowing us to estimate nodes 5 and 8 as shifted nodes. Identifying the elbow point has the advantage to detect shifted nodes without relying on a fixed threshold.
A.3 Algorithm details for FOCI
Algorithm 4 Feature ordering by conditional independence (FOCI)
Input: Data Pre(X))) X²
Output: Estimated parents of X,, PA
1: P+( (
2: Let T (i
X X/)
Tm is the estimator in Azadkia and Chatterjee [3].
3: while max X,EP, XxEFre(X,) Ten (i, P)
do
4: P PU arg max EP X.EPre (m(
5: PA b P
'In practice, one could omit from the denominator as it is unusual to obtaín mins Varî 0 from finite samples and computational precision.
17

---

60
50
variance rabio
40
30
20
10
node index
Figure 5: Statistic in eq.(8) for each node sorted in non-increasing order. In this case, node index corresponds to the elbow point, allowing us to estimate nodes 5 and 8 as shifted nodes.
B Detailed Proofs
B.1 Proof of Theorem
To prove Theorem we will make use of the following lemmas.
Lemma 1. Let {an) H and (bn Hi be two sequences of real numbers, where on > 0, Vh. Then
have:
H
Σ anbi
Σ Ob
Σ
anhn
2 0,
hml
with equality if and only 6; bj Vj tiiie [M].
Proof. We can invoke the Cauchy-Schwarz Schwarz inequality with vectors a (Vai,.... vaz). and (bi by waw.aaann /an), then we have:
(u T)) s | |2/22111
which proves the inequality. The equality holds if and only if u and v are linearly dependent, i.e., when b, by for all # j [H].
Lemma 2. For any if (X, PAA5) 1p PAA ), then df =cty.
Proof. Denote the associated density of (X, PAA) when X, = 2y as phi, (x, 13)()))) and
let 2ÿ (PA")
8²
log pK (z, (PA())
0zz))²
0 log Vi, (u) 8u +
8 log p'N (u)
du Dag
Ou
0(z,)²
Ou
=0 d/ = =]
where we the fact that 85, pa' (X, PAA)), then oft =
8Pu 1,
Then it immediate follows that if (X(, PAA)
a log Do (2z Ju, f/ (PA))) )) if and only iX, PA7) pM (X, PAA ), where p² and pd are the probability density functions corresponding the probability measures Ipb and pp when X3 Zj.
Lemma 3. For any under Assumption B. da,
a log (x,- - (PA)))
18

---

Proof. Denote the associated density of (X, PAA) when X, Iy as PkA (z3 (PAA)) we proceed as follows:
①
log p (I, 1/ = a
dz, log kt' (x, 83 (PA))) logrk,
- logp's (2) - f3 (PAA)) = log pKx (2g 13 (PA )) + const Pkx, (2g (PA))) = pkx (x, f3 (PA) n 1) E
Pli, (x) -5 (PAA) connet / 1. [1g 14 (PA5) 1)dr,
JR
/ Pl, (z) - (25 -f) -f (PA()) = E const [ 1h (T, - 13 - S3 (PA)) ))
+1 const
consl = 0
Here, const a constant that is independent of zy. Integrating both sides with respect to æ, and using the fact that nd (x)dx we conclude that const = 0. Hence, we can establish the following:
log p', (z, PN, - (PA))) = dz, log ( - 13 (PA ))
- v&, (ag - (PA))) = 2\x (c, - 13" (PA (z (X, PAA) = 1pd (X, )
Proof of Theorem 1. Let us first expand the log density of the mixture distribution:
H
log (z) = llog
Σ
wkp (x))
h=1
Then, recall that s(x) log q(r), the j-entry reads:
whp"(x)
① log (og PPA]) + Σ Ô dz,
8z(x)= =
logr (z; PA!)))
h=1 k= (x))
SECHA
h
whp" (x)
0 logp't, () + Σ
i log (zi (o S4 (((At)
1= wap (x)
Ox
Or,
iECH3
Condition (i). First will prove condition (i). That is, given leaf node X, in all DAGs Gh, Xj is not a shifted node (i.e. an invariant node) if and only if Var( @e,(*), 0.
If xy is a leaf node in all the DAGs Gh,, then CH² = 0, vh E [M]. and we can write (9) as:
whp (x)
(y(x)= =
Ee weph(r) dx, log ph, (3
We use Den (x)(z) and Num Bay tively. Then we have:
(st(x) dry
to denote the denominator and numerator of (s1(z) Bny respec-
2
Den (o,(x)) = Ij
H
Σ
whp (x)
19

---

0 Omy
log PN, (z, - 1/ (PA))) + (x) (
a
log Logr'' (2, )(PA311)
Num
why (x)
J.j
Ony
H
H
x
WkP
(x)
Σ
wpp (x)
dz, log PN, (zy- - f3 (PA9))
k=1
Now, dividing Num (m) over Den( (x, du,
) we obtain:
dz,
(x)
Num (ar(ar)
H
wpp (x)
0²
log logph, (zy (PA)))
Ozj
Den
(xy
k=1 uxph (x) Oxy
wap"(x) (x)
8 Oz log logph (z) - 53 (PA))))
h= 1=1 (a)b(()
12
2
wrpp (x)
t 12 (2) E uxp (x) dx, log
(10)
Note that since 2j & PAG. the function f3 (PA)) is independent of Tj. Let On = wxp"(x) and let b Jæz RHS of (10) can be written as:
a log r/,, (x) - f3 (PA))). Then, the last two summands of the
H an
Σ
anbi
Σ an by
(11)
ΣHH an) )²
where the last inequality holds from Lemma Then, by Lemma 3, we have that by = by Ish (X, PAA) pph (X; PAA)) for all h, h1' E {H). Then if bn = by holds, by Lemma 2, we have c) = - = Cg and then the first term for (10) boils down to constant Cj. Finally, from Lemma 1, we have that equality in (11 holds and only if by by for all h, h [H]. Thus, we conclude that:
0sy(x) 0x.
If Xj is leaf node for all Gh, then X, a shifted node
Gj.
where (as(2) 8x, = ey is equivalent to Var,( (2 0z, = O.
(u) Condition (ii). We now prove that if Var,( 02, > O, then only one of the following two cases
holds: Case 1) Xj is a leaf node for all Th and shifted node. Case 2) Xj is not a leaf node in at least one DAG Gh.
Case follows immediately from the proof of condition (i) above.
For Case 2, we study whether there exists a non-leaf node X, with Var,, partial derivative of (I) in (9) w.r.t. Fj, we have:
( (8,(()))
O. Taking the
(x)(x) dr,
(x) 1= wayk (xr)
0² On?
0 log logr's (z. 1 (PA!))
log logr', (2) (PAR))
Σ
On?
iECH
2
1 (x)
0 Ozy
θ dzy
log ph; (2, - (PA)))
Σ SECH
log logrks (2: (2:- (PAA))
A= (x)
20

---

log Dog.N. (2; ( - ( |
(x) k= (xp
a doy
log logr2; (ay 13 (PAA)) iECHy
Σ
a dx,
By Assumptions B, we have 0e duy
logr' (5g - (P(A)) c/. For simplicity, let an a log p\,, (0; (PAA)) Then, we have:
Daj log Lgg,,, (xz - 13 (PA))) + Eiechh Je,
(3) whp" (x)
H
whp (z)
8
c2
Σ
Logp' (zt (PA\))
0z;
k=1 wk (z)
k=1 Wkp (x)
dzy
iECHş
H
whp (x)
whp wa) ah
o
Σ h=1
(12)
-= Wpp (2)
h=1 k=1
2
(z; We prove that L (x; (PA!)) is not constant under any circumstance, by contradiction. Let Gh be an environment's DAG where Xj is not a leaf, and let X, E CH² such that X.. & UsECHY PA Note that X,, always exist since X, is not a leaf, and it suffices to pick a child X appearing at the latest position in the topological order of Gh,, Now suppose that d 0ny
log Logri, (Tt
1h (PA = a, where a constant. Then we have:
i log ph, (Iu - (PAA)) = az, + g(x-5),
0z, M (PA))
0
log (mhe) az, g(x
On,
By deriving both sides w.r.t Tu. we obtain:
0² logpl. (nw) 0g(x=;)
dz, ft (PA
ong
Oxa
0g(x--))
f (PA c/ 0zy
Oer
ash Since the RHS does not depend on Tj, then dz, cannot depend on æÿ neither, implying that fl is linear in Iy. thus contradicting the non con-linearity assumption (Assumption A). Consequently, it becomes evident that term cannot be a constant, regardless of whether the node X, has undergone a shift not.
Now let us take a look to term (1.). We have:
2
wap' (x) H Ee wep* (x)
whp" (x) k=
Σ
Σ
ah
20
(x)
where the inequality follows by Jensen's 'ss inequality. Thus we conclude that X, is non-leaf node, we have Var,( (0, )> 0.
B.2 Proof of Theorem 2
To proof the theorem we will need the following assumptions:
Assumption C. Let MB,, denote the Markov Blanket of node X, under environment h. then assume
21

---

There are non-negative real number 3 and C' such that for any subset Xs C Pre(X)) of size S 1/6 + 2, any I, '' R1Xas and any 1.
P 2t
2 Xs = =x)) 1s) C(1 ||>|| + |<'1| - x'll
where (Xs] is the size of the set Xs.
There are positive rumbers C₁ and C2 such that for any Xs of size S 1/8 + 2 and any P(||Xs| 2 t) 2 Cie-Cst
For any subset Xs S (r)((X)) such that Xs &. MB there exists X; with X; E MBJ \ss. such that for any Xj with X, & MB
6/4 6/4 Q(xs) = - / Var(P(x) 21 (xsss)ddp(t
Q(xsu - Q (Xx))
where 5 the largest number such that for any subset X from Pre(X)), there some X₃ $ Xs such that Q( {X+}) 2 Q(Xs) + 6.
Proof. Under Assumption C, from Theorem 3.1 in [4] we have
P(MB = MB3)
Cym-Cam
where C3 and C4 constants that depend only the data generation process, and m is the number of samples. Since the estimated topological order से is assumed to valid for all environments, we can conclude that node Xj leaf node in the input data {Pre(X²), X for all h. As result, have MB; PA based on the Markov blanket definition. Therefore, the output of Algorithm is equal the true parent set PAA with high probability.
B.3 Proof of Theorem 3 in Appendix D
To prove Theorem we will make use of the following lemmas.
Lemma 4. Suppose X is an kn dimension matrix, Y is an x ky dimension matrix. Let the columns of the concatenated matrix Z (X,Y) be linearly independent. Consider Ba, de-dimensional vector, and By and By, both ky di-dmmnniional vectors. If B. + B, Y Sy. then follows that By = By and 8g = 0.
Proof.
3x
(ks Bg
B.
X 3, + YB, YB, = (X,Y)
X
Y)
2
0
B,
By
Since Z has full column rank, then the null space of Z is 0, which implics Br = 0, B, = By- Lemma 5. For any [11], if
Σ
V
kEPre(Xx)
kEPA
then 3/h B/, if PAJ, and B/, if & PAA.
Proof. Rearrange the set Pre(X,) so that Pre(X,) {Xk, X ... X Xkm+1..... X where (X X Pre (Xx) PAj, and {X km Xx,.) PAy. Then let X (TT ....... V kmm Y kmmm41:-- IK,,), and Z (X By the linear independence prop- erty of the basis functions of Pre(X;). we have that the columns of Z are linearly independent. Also. let
3h
2h Pyhm+z
B =
By =
By =
kmm
jikk
22

---

Then we must have:
Σ
Σ
X + Y = BB.
k)))
kEPAy
Then by Lemma 4, we have 3h = B
if E PAY. D 0ffk
Proof of Theorem 3. We know that Pre(Xx) contains the ancestors of Xj. Then, in environment h, we have:
Ext(X Pre(XX]] = 1ff (PA)) Pre(X;)] + B Pre(X,)]
= (PA)))
where the last equality follows since the first conditional expectation is equal to (PAy), due to PA, Pre(Xy). Moreover, in the second conditional expectation term, we have that N₁ is marginally independent of Pr)(X,) by the d-separation criterion. Thus the conditional expectation of N₃ equals to the marginal expectation of Nj. which is O. Finally,
Σ
V
2h Ski
Rag [X, Pre( X,) (PA))
Σ
V
kEPro(X,)
REPA
By Lemma 5, we have 3h
S if E By =0 if PA
C Additional Experiments
This section provides a thorough evaluation of the pipeline of our method. We begin by assessing the performance of our method in detecting the shifted nodes. Subsequently, extend the evaluation to include the recovery of the structurally shifted edges.
C.1 Experiments on detecting shifted nodes
Graph models. We ran experiments by generating adjacency matrices using the Erdōs Erdös-Rényi (ER) and Scale free (SF) graph models. For a given number of variables d, ERk and SFk indicate an average number of edges equal kd.
Data generation process. We first sampled a Directed Acyelic Graph (DAG) according to either the Erdös-Rényí (ER) model or the Scale-Free (SF) model for environment &
For environment &g, we used the same DAG structure in environment ε,, ensuring a direct comparison between the two environments. To introduce artificial shifted nodes, we randomly selected 0.2 d nodes from the non-root nodes, where d represents the total number of nodes in the DAG. These selected nodes were considered as the "shifted nodes," denoted as S, with [S] = 0.2d.
The functional relationship between node X, and its parents in environment & defined as follows:
X, Σ sin(X?) Nj,
iEPA,
while for environment E2, we defined the functional relationships between each node and its parents by:
X, EiEra, sin(X? + Nj, EiEPA, coos(2X? - 3X.) + N5, X, E.
if X, S
Experiment detail. In each simulation, we generated 500 data points, with the variances of the noise set to We conducted 30 simulations for each combination graph type, noise type, and number of nodes. The running time was recorded by executing the experiments on an Intel Xeon Gold 6248R Processor with 8 cores. For our method, we used the hyperparameters et G 0.005. eta H 0.005, and threshold f 2 (see Algorithm 3).
23

---

Evaluation. We conducted comparative analysis to evaluate the performance of our method in detecting shifted nodes compared DCI. The evaluation was based on F1 score, precision, and recall as the evaluation metrics. Furthermore, we examined the robustness of our method by conducting tests using Gumbel and Laplace as noise distributions.
Figures 6, 7, and illustrate our method's performance across varying numbers of nodes and sparsity levels of the graphs. Our method consistently outperformed DCI in terms of Fl score, precision, and recall.
Gaussian
Gumbel
Leplace
Gaussian
Gumbel
Laplace
100
75
75
F1 Score (%)
50-
ER2
Precision (%)
50
25-
25
0-
0
100
100
75
75
50
50
3
25
25
0-
0-
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
10
20
30
40
20 d (number nodes)
50
10
30
40
50
10
20
30
40
50
(number of nodes)
Gausslan
Gumbel
Laplace
Gaussian
Gumbel
Lapiace
10
Runtime (seconds)
233
Rocall (%)
0-
1001
10
75
2:8
50
25
0-
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
10
20
30
40
60
10
20
30
40
50
(number nodes)
(number of nodes)
Methods
DCI
ISCAN (ours)
Figure Shifted nodes detection in ER2 and SF2 graphs. For each point, we conducted 30 simulations as described in Section C.1. The points indicate the average values obtained from these simulations, while the error bars depict the standard errors. For each simulation, 500 samples were generated. Our method iSCAN (green) consistently outperformed DCI (red) in terms of Fl score, precision, and recall.
24

---

Gaussian
Gumbel
Laplace
Gaussian
Gumbel
Laplace
F1 Score (%)
Precision (%)
50-
25
25
100
100
75
75
50
4:35
50
7:8
25
25
0
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
10
20
30
40
50
20 (number
10
30
40
50
10
20
30
40
50
(number
nodes)
nodes)
Gaussian
Gumbai
Laplace
Gaussian
Gumbel
Lapiace
100
20
75
Runtime (seconds)
15
10
Recall (%)
25
100
20
F-S
15
PJS
10
30 10
40
10 d (number of nodes)
50
10
30
50
30
40
50
30
50
10
30
50
40
(number of nodes)
Methods
DCI
ISCAN (ours)
Figure Shifted nodes detection in ER4 and SF4 graphs. For each point, we conducted 30 simulations described in Section C. The points indicate the average values obtained from these simulations, while the error bars depict the standard errors. For each simulation, 500 samples were generated. Our method iSCAN (green) consistently outperformed DCI (red) in terms of F1 score, precision, and recall.
Gaussian
Gumbel
Laplace
Gaussian
Gumbel
Laplace
75
F1 Score (%)
50-
Precision (%)
50-
ER6
25
25
0-
100
75
75
50-
9:55
60-
9F5
25
25
0-
0
20
30
40
50
20
30
40
50
30
40
50
20
30
40
50
20
30
40
50
20
30
40
50
(number nodes)
(number nodes)
Laplace
Laplace
100
75
40
EW6
Runtime (seconds)
Recall (%)
20
40
9-5
9-55
20
20
40
50
20
40
50
40
60
30
40
50
20
60 20
30
40
50
d (number of nodes)
(number of nodes)
Methods
DCI
ISCAN (ours)
Figure Shifted nodes detection in ER6 and SF6 graphs. For each point, we conducted 30 simulations described in Section C. The points indicate the average values obtained from these simulations, while the error bars depict the standard errors. For each simulation, 500 samples were generated. Our method iSCAN (green) consistently outperformed DCI (red) in terms of Fl score, precision, and recall.
25

---

C.1.1 Experiments in detecting shifted nodes from Gaussian process
Data generation process. We first sampled DAG according to the ER or SF model. In our experiment, we considered two environments, & and Eg, with the same DAG structure. Each node in the graph had a functional relationship with its parents defined as = is an independent standard Gaussian variable. Recall that the superscript h denotes the function for environment E
+ Nj, where N;
To introduce shifted nodes, we randomly selected 0.2 - d nodes from the non-root nodes, denoted as S, be the shifted nodes. In other words, |S| 0.2d. For the non-shifted nodes Xj (i.e j & S), we set fj fj. However, for each shifted node X, in S, changed its functional relationship with its parents to X, = 2- 2- ()(()) + Nyz-
To test our method in a more general setting involving nonlinear functions, we followed the approach in [65, 39]. Specifically, for non-shifted nodes, we generated the link functions f] by sampling Gaussian processes with a half unit bandwidth RBF kernel, and we set f3 = S}. For shifted nodes, X, E S. we generated the link functions fj and f by sampling Gaussian processes with a half unit bandwidth RBF kernel independently. This allowed us simulate different functional relationships for the shifted nodes across the two environments.
Experiment detail. In each simulation, we generated 1000 data points, with the variances of the noise set to 1. We conducted 30 simulations for each combination of graph type, noise type, and number of nodes. The running time was recorded by executing the experiments on an Intel Xeon Gold 6248R Processor with 8 cores, For our method, we used the hyperparameters et G = 0.005, eta_ H 0.005, and elbow True (see Remark 6).
Evaluation. We conducted a comparative performance analysis between our proposed Algorithm 1 (iSCAN, green) and the DCI (red) method. The results for ER2 and SF2 graphs under Gaussian, Gumbel, and Laplace noise distributions are shown in Figure 9. In certain cases, our method may underperform DCI in terms of precision, resulting in lower F1 score. However, it is important to note that our method consistently outperforms DCI in terms of recall score.
Furthermore, Figure 10 and Figure 11 present the results for ER4/SF4, and ER6/SF6 graphs. In terms of precision, our method exhibits competitive performance and, in many cases, outperforms DCL. Notably, iSCAN consistently surpasses DCI in terms of recall score and Fl score.
These findings emphasize the strengths of our proposed method in accurately detecting shifted nodes and edges, particularly in terms of recall and overall performance. In denser graphs, our method demonstrates superior ability to recover shifted nodes compared to DCI. This suggests that our method is well-suited for scenarios where the graph structure is more complex and contains larger number of nodes and edges. The improved performance of our method in such settings further highlights its potential in practical applications and its ability to handle more challenging tasks.
Top-k precision. We have observed that in some cases, the precision of our method underperformed DCI. We attribute this to the elbow method rather than statsr. To further investigate this, we conducted an analysis using only stats,, and measured the precision based on different criteria. Specifically, we identified nodes shifted if their statsr, ranked first, first two, or within the top k, denoted as top-1 precision, top-2 precision, and top-k precision, respectively, where A: = [S].
Figure 12 presents the results of precision for top-1, top-2, and top-k criteria under various graph models and noise combinations. In most cases, the precision exceeds 80% and even approaches 100%. These results indicate that when using stats alone, our method still provides accurate information about shifted nodes. The findings suggest that the lower precision observed Figure can be attributed to the elbow strategy rather than the effectiveness of stats Overall, this analysis strengthens the reliability and usefulness of stats, S r. in accurately identifying shifted nodes in our method.
26

---

Caussian
Gurrbel
Laplace
Gaussibn
Gumbel
Laplace
90-
80
60
2333
5.orr (19
50-
Precision (%0
40%
90-
80
70
2.55
50
50-
42-
30
10
30
40
10
30
40
50
10
30
50
10
20
30
50
10
30
30
50
20
30
40
50
(number of nodes)
(number nodes]
Geussian
Gumbel
Laplace
Gaussan
Gumbel
Laplece
50
200
Rurning tire (Seconds)
50
Recall (%)
40
50
0:4
1.00
30
10-
50
10
30
50
20
40
10
30
30
20
50
10
10
20
30
50
6 [number of nodes)
[number of nodes)
Methods
Figure 9: Experiments detection of shifted nodes in ER2/SF2 graphs using Gaussian processes. Details described in Appendix C.1.1. The error bars represent the standard errors.
Caussian
Gurrbel
Laplace
Gaussian
Gumbel
Laplece
90-
70
50-
60
Score (%)
Precision (%0
50-
40
80
60-
80
5/5
40
20
10
30
40
50
10
30
40
50
10
20
30
50
10
20
30
40
50
10
30
30
40
50
10
20
30
40
50
(number of nodes)
[number of nodes)
Gaussien
Gumbel
Lablace
Gaussian
Gumbel
Laplace
50
200
Rurning time (Seconds)
V112
Recall (1))
50
100
50
20
15
30
50
20
40
10
30
30
20
30
30
10
20
08
10
30
50
d (number of nodes)
[number of nodes)
Methods
ISCAN
Figure 10: Experiments detection of shifted nodes in ER4/SF4 graphs using Gaussian processes.
Details described in Appendix C.I.I. The error bars represent the standard errors.
27

---

Gaussaan
Gurbbil
Lapace
Gaussion
Gunbel
Laplace
80
80-
60
60-
42-
Score (%
Precision (%6)
40-
20-
80
80-
60-
60-
9.15
20
20g
30
20
30
60
50
10
36
40
30
30
20
30
30
30
20
30
40
50
10 20
30
40
50
20
30
40
30
(number of nodes)
(number of nodes!
Laplace
100
200
Running time (Seonds)
1.50
51123
40
1.00
Recal (%)
250
200
60
9.255
6 (number of nodes)
(number of nodes)
Methods
ISCAN
Figure 11: Experiments on detection of shifted nodes in ER6/SF6 graphs using Gaussian processes.
Details described in Appendix C.I.I. The error bars represent the standard errors.
300
300
300
80
60
80
50
60
60
40-
40-
40-
100
RD3
40-
Top Pro
21/3
Procision
Top Pracistin
150
200
2.00
60-
60
60
60-
40
40
100
100
200-
10
10
80:
912
60
60
40
10 2 30
13 20 30 10 20 40.40050 jiuumberr modes)
10
50
20 40 50 (nuumbr modes)
35 50
10
12 20 50 (nuumbe nodess)
50
M Wethods
(SCAN (ound)
Figure 12: Top and K performance of iSCAN where functionals sampled from Gaussian processes. Details described in Appendix C.1.1. The error bars represent the standard errors.
28

---

C.2 Experiments on estimating structural shifts
Data generation. We first sampled DAG, G', of nodes according either the ER or SF model for env. &,. For env. Eg, we initialized its DAG structure from env. & and produced structural changes by randomly selecting 0.2 nodes from the non-root nodes. This set of selected nodes S. with cardinality |S| 0.2d, correspond to the set of "shifted nodes". In env. &g, for each shifted node X, E S. we uniformly at random deleted at most three of incoming edges, and use D2 to denote the parents whose edges X3 were deleted; thus, the DAG G2 is subgraph of G Then, in E,, each X, was defined as follows:
Xj
Σ IEPA]ID,
sin(X?)+ + Σ cos - 3X,) + N
iED;
In &g, each X, was defined as follows:
X,=
Σ REPA3
sin(X?) + N₁
Experiment details. For each simulation, we generated 500 data points per environment, i.c., mi 500, my 500 and m 1000. The noise variances were set to We conducted 30 simulations for each combination of graph type (ER SF), noise type (Gaussian, Gumbel, and Laplace), and number of nodes (d E {10, 20, 30, 50)). The running time was recorded by executing the experiments on an Intel Xeon Gold 6248R Processor with 8 cores. For our method, used (see Alg. 3).
7 = 0.05 for (q. and (7) and a threshold
In the case of the method introduced by Budhathoki et al. [11], we employed Kernel Conditional Independence (KCI) tests [86] for conducting conditional independence tests. As for CITE, KCD, UT-IGSP, and SCORE, we used their respective default parameter settings provided within their packages. Additionally, for SCORE, we employed estimate the DAGs independently for different environments and then compared the recovered DAGs to identify the shifted nodes. Given that Budhathoki's 's and KCD methods require information about the parents PA, for each node Xj. we employed the SCORE method to find the parent sets PAj.
Evaluation. In this experiment, we assessed the performance of our method in two aspects: detecting shifted nodes and recovering the strucutral changes (difference DAG). For the evaluation of shifted node detection, measured Fl score, recall, and precision. In the evaluation of difference DAG recovery, we compared the estimated difference DAG with the ground truth difference DAG using Fl score. Additionally, we considered the running time of the methods as another evaluation criterion.
Figures 13, 14, and 15 illustrate our method's 's performance in detecting shifted nodes across varying numbers of nodes and sparsity levels of the graphs. Our method consistently outperformed baselines in terms of F1 score, precision, and recall.
Figures 16 showcase the performance of our method in recovering difference DAG across different noise distribution, different numbers of nodes and sparsity levels the graphs. Our method achieves higher F1 score in recovering the difference DAG compared with DCI.
29

---

Laplace
Gaussian
Laplace
100
40
Recall (%)
300
2.0
80
00-
555
30
20
30
40
20
20
30
40
50
10
20
30
30
30
20
30
40
50
30
20
30
40
50
30
20
30
40
[number of nodes)
Gaussan
Gumbel
Laplace
Gauto an
Gumbel
Laplate
300
800
Prectsion (%)
200
40
Runtime
000
L200
0
400-
200
10
20
30
40
50
10
20
30
40
50
15
20
30
50
20
10
50
20
40
(number of nodes)
Inumber of nodes)
Methods
Buchathoki
CITE
DCI
KCO
ISCAN (ouns)
SCORE
UT-IGSP
Figure 13: Shifted nodes detection performance in ER2/SF2. See App. C.2 for experimental details.
iSCAN (light blue) consistently outperformed baselines in terms of F1 score, precision, and recall.
Laplace
Laplace
80
80
PHB
60
Recall (%)
core
20
2.0
80
60-
P-5
30
20
20
40
30
20
30
40
50
10
20
30
30
20
20
30
40
30
20
30
40
30
30
20
30
40
30
Inumber of nodes)
Gaussinn
Gumbel
Laplace
Gau55 an
Gumbel
Laplate
800
60
200
Predsion
40
Rurtime
000
4000
3000
2000
40
1000
30
20
30
40
10
20
30
40
50
10
20
30
25
10
30
40
50
20
40
(number of nodes)
Inumber of nodes)
Methods
Buchathoki
CITE
DCI
KCO
ISCM (ours)
SCORE
UT-6SP
Figure 14: Shifted nodes detection performance in ER4/SF4. See App. C.2 for experimental details.
iSCAN (light blue) consistently outperformed baselines in terms of Fl score, precision, and recall.
30

---

Lasbee
100
80
Recall (%)
40-
100
100
90
80
80-
70
0.5
80
60
40-
10
30
20
30
20
20
30
20
40
10
33
40
50
[number of nades))
100
800
Lapace
Runtime 65)
200
Precisoon
2004
5000
3000
2000
40
1000
10
30
50
10
30
50
20
10
30
50
10
30
30
Inumber nodes)
(number nodes]
Methods
BudhathoM
D21
ISCAN (uur))
SCORE
UT-IGSP
Figure 15: Shifted nodes detection performance in ER6/SF6. See App. C.2 for experimental details.
iSCAN consistently outperformed baselines in terms Fl score, precision, and recall.
Geussian
Gumbel
Laplace
Gaussian
Gumbel
Laplace
Running Liee (Secords)
30
(96
20-
10-
50-
10 20
20
30
40
50
10 30 40
30-
200
350
100
10-
50-
10
20
30
40
30
10 20
30
20
30
20
30
40
30
30 40 30 (number nades)
10 20
30 20 30
(number of nodes)
Gaussaan
Gumbel
Laplace
Gaussion
Gunbel
Laplaca
50-
Rurning ( (Secordts)
20-
200-
30-
10
40
30
40
30
20
60
30
10
20
30 40
50-
250
40-
200
0.25
20-
10
10
20
30
40
30
20
20
30
40
40
50
50 (number nodes)
10
20
30 40
10
20
30
40
50
(number of nodes)
Laplace
Rurning time (Seconds)
300-
30
953
20-
200
10
20 30
40
50
10
30
30
40
50
10
20
30
40
10
20
30
40
50
10
20
30
40
50
10
70
30
40
50
F35
30-
40-
300
9.55
25
(number nodes)
(number nodes)
Methods
0C1
SCAN-FOUI
Figure 16: Difference DAG recovery performance in all different graphs. iSCAN-FOCI (green) consistently outperformed DCI (red) in terms of Fl score.
31

---

C.3 Performance of Alg. 3 using the elbow method
In this section, we aim to understand the performance of our method when using the elbow approach discussed in Remark 6, random functions for shifted nodes, and different noise variances per variable within an environment.
Data generation process. We first sampled a Directed Acyclic Graph (DAG) according to either the Erdös-Rényí (ER) model or the Scale-Free (SF) model for environment &
For environment Eg, we used the same DAG structure as in environment ε,, ensuring a direct comparison between the two environments. To introduce artificial shifted edges, we randomly selected 0.2 d nodes from the non-root nodes, where represents the total number of nodes the DAG. These selected nodes correspond to shifted nodes, denoted as S, with |S] 0.2d. For each shifted node X, E S, we uniformly and randomly deleted 3 edges originating from its parents for environment E.. The parent nodes whose edges to X, were deleted denoted as Dj.
The functional relationship between shifted node X, and its parents Dj in environment & was defined as follows:
Xj
Σ i
sin(X?) + Σ cij 1u(-2X7 + +3X? + 4X,) + N IED;
where Cty ~ Uniform -001,,,-2-111 2 5. and fiy is function from (sinc(-), cos(-)) chosen uniformly random. For environment Eg, where the adjacency matrix has undergone deletions, defined the functional relationship between each node and its parents as follows:
X, Σ PEPA,
sin( + N₂
Experiment detail. In each simulation, we generated {500, 1000} data points, with the variances of the noises set uniformly at random in [0.25, 0.5]. We tested three types of noise distributions, namely, the Normal, Laplace, and Gumbel distributions. We conducted a 100 simulations for each combination of graph type, noise type, and number of nodes. The running time was recorded by executing the experiments on an Intel Xeon Gold 6248R Processor with 8 cores. For our method, used the hyperparameter 1) = 0.001. Different from the hard threshold of = 2 used in previous experiments, now used the elbow approach to determine the set of shifted nodes. To automatically select the elbow we made use of the Python package Kneed', with hyperparameters curve- 'convex', di direction= 'decreasing' online=online, interp_method= method= 'interpd''.
Evaluation. In this experiment, we assessed the performance of our method in two aspects: detecting shifted nodes and recovering the difference DAG. For the evaluation of shifted node detection, we measured F1 score, recall, and precision. the evaluation of difference DAG recovery, we compared the estimated difference DAG with the ground truth difference DAG using F1 score.
In Figures 17 and 18 we present the performances when using the elbow approach discussed in Remark 6. In Figure 17, we note that iSCAN performs similarly for number of samples 500 and 1000. We also show the top-1, top-2, and top-k precision of iSCAN when choosing the first, first and first k variables of stats s (see Algorithm 3) after sorting in decresing order, respectively. We remark that the superbly performance of iSCAN in top-1 or top-2 precision suggests that in situations that is difficult to choose threshold for Algorithm 3, the practioner can consider that the first or first two variables of stats are more likely to be shifted nodes. Finally, in Figure 18 we show that iSCAN outperforms DCI in recovering the underlying structural difference.
We used the latest version found
32

---

Laploss
Gumbel
Gauesien
Oumoel
Laplsce
75
50
50-
25
25
0
75
75-
50
60-
50-
100
100
100
75
75
75-
F1 Soore (%)
Recail (86
Precision
100
100
26-
0
75
50-
25-
100
100
75
76
50
50-
50-
d0 10 20 30 do 20 40 (number nodes)
10
do 50 10 30 do 50 20 30 (number nodes)
20 30
50 50 30 40
(number nodies)
Gaussian
Lagiace
Laplace
100
80
80-
80
80-
P6
80
60-
100
1-2
100
1.00
Top Precision
80
Top Precision
663
Top Precision
60
80
60
100
100-
100
2:50
80-
5:55
60
60
80
80-
80
9d5
80
BJ66
20 30 40 50 20 30 50 20 30 40 50 (number of nodes)
10 20 40 1010 20 30 40 50 10 20 30:30044
10 20 30 50 10 20 30 5010 20 40.50
d (number of nodes)
(number of nodes)
Methods
(ours)
Figure 17: Shifted nodes detection performance in ERk and SFK for & {2, 4, 6). For each point each subplot, we conducted 100 simulations as described the average values obtained from these simulations. The error bars depict the standard errors. Our method iSCAN (green) consistently outperformed DCI (red) in terms F1 score, precision, and recall.
Section C.3. The points indicate
33

---

Gaussian
Gumbel
Laplace
Goussin
Gurbel
Laplece
Gaussian
Gumbel
Laplace
50
30
50
25
25
25
04
50
50
PM
50-
25
28
25
04
50
25
Recail (%
25
Soore
0-
Precision
0-
60-
25
25-
90
80
26
25-
01
0
50
26
0
d0 60 do 20 so (number nodes)
30 20 40 sori 20 50 50
20 40 30 10 30 (number nodes)
(numbee nodes)
Methods
ISCAN + FOCI
500
1000
Figure 18: Difference DAG recovery performance in all different graphs. For each point in each subplot, we conducted 100 simulations as described in Section C.3. The points indicate the average values obtained from these simulations. The error bars depict the standard errors. Our method iSCAN with FOCI (green) consistently outperformed DCI (red) in terms of F1 score.
34

---

D Additional Discussion on Shifted Edges
In Section 4, focused estimating structural changes across the environments (Definition 4). However, in some situations it might be of interest to determine whether the functional relation- ship between two variables has changed across the environments. The latter could have multiple interpretations, in this section, we elaborate on particular type of functional change via partial derivatives.
Definition 5 (functionally shifted edge). Given environments & = (X, fth PN) for h E [H), an edge (X; Xy) is called functionally shifted edge if there exists h N [H] such that:
0 Oz.-
(PA)))
2 Oz,"
1y(PA))
Without further assumptions about the functional form of ff,, certain ill-posed situations may arise under Definition 5. Let consider the following example.
Example 1. Let A and Ê be two environments, each consisting of three nodes. Let the structural equations for node X3 be: =3= exp (Xi + Xx²) Na, and X exp (2-X + x + N3. In this scenario, one could consider that the causal relationship X2 + X3 has not changed. However,
we note that 0f1 3ug auy. thus, testing for changes the partial derivative would yield a false
+ 0f2
discovery for the non-shifted edge X2 X3-
Ill- Ill-posed situations such as the above example can be avoided by additional assumptions on the functional mechanisms. We next discuss a sufficient condition where the partial derivative test for functional changes is well-defined.
Assumption D (Additive Models). Let S be the set of shifted nodes across all the H environments. Then, for all j E S, E [H]:
(PA)) =
Σ (ax(),)
REPAS
where a is a constant. fyk is a nonlinear function, where ffx()) (-) lies in some space of function class F.
Remark Assumption D amouuts modelling each variable as a generalized linear model [27]. is widely used in nonparametries and causal discovery [12, 43, 80). Moreover, it not only provides practical framework but also makes the definition of shifted edges (as per Definition 5) well-defined and reasonable.
Remark 8. Note that Assumption D makes assumptions only on the set of shifted nodes. This is because the set of invariant nodes can be identified regardless of the their type of structural equation, and it is also clear these nodes cannot have any type of shift.
Now consider a function class F, which incorporates the use of basis functions to model the additive components f1x- Specifically, we express ffe (Ixx)) V ( where feature mapping Vi is a x r matrix whose columns represent the basis functions and 5hh is an r-dimensional vector containing the corresponding coefficients. Moreover assume that the functions fJk...,, share same feature mapping V (x =,.... ... = (-) but can have different coefficients 8h across the H environments. The latter has been assumed in prior work, e.g., [44]. The approach of using basis function approximation is widely adopted in nonparametric analysis, and it has been successfully employed in various domains such graph- methods [80], and the popular CAM framework [12]. Then, under Assumption D and Definition 5, we present the following proposition:
Proposition 2. Under Assumption D, an edge (X; Xj) is a functionally shifted edge, as in Definition 5, if and only if the basis coefficients are different. That is, Ofh''
0fh + Ox, dz,
= B/ +
3/
Proof. We have,
0fh (a.) d(w d(H))(z))))))
ddd ( g1. dr;
=
dx,
dx; dz;
35

---

Then.
05 05 Or;
dV (x()) (0)
+
3
- - +
dr.
dr;
The last < relation is due to the linear independence of the basis functions v je, then the null space of // can only be the zero vector
Note that the output of Algorithm 1 also estimates topological order f. However, the exact parents of a node Xj across the environments are not known, and they are possibly different. To estimate the coefficients without knowledge of the exact parents, we can consider the set Pr)(X)), which consists of nodes located before Xj the topological order By regressing Xj on Pre(X;) for each environment, can obtain coefficient estimations, which are the same coefficients obtained by regressing X, on its exact parents, in large samples.
Theorem 3. In large samples, let { (REPTe(X))) be the coefficients obtained by regressing X₃ the feature mapping of Pre((;)). and let (3) I the coefficients obtained by regressing X; on the feature mapping of PAg. Then, 8/ S/ if PAJ, and 8 if Pro((;) \PA
Proof. Proof can be found in Appendix B.3.
Motivated by Theorem 3, and given estimated { [ )ePre(X,)) could conduct hypothesis testing as follows:
Ho 33 = ........ 3/
(13)
If the null hypothesis Ho is rejected, it indicates that there is evidence of functionally shifted edge between nodes X₂ and Xj across the environments. In this paper we leave the hypothesis test unspecified to allow for any procedure that can test (q.(1
Algorithm Functionally shifted edges detection
Input: Sample data X ..... X shifted nodes set S, topological order स,, significance level a. Output: Set of functionally shifted edges E
1: for j E S do
2: Estimate 5i1 for all & Pre(X,) and h E [H]
3: for E Pre(X,) do
4: Conduct hypothesis testing Ho (equation 13) under significant level a. If Ho rejected, add edge (Xg Xj) to E.
5:
36


Reviews:
[{"review_id": "zMy1K3zGsW", "reviewer": "Reviewer_m244", "summary": "The paper proposes an efficient method to identify causal mechanism shifts across various interventional environments. The method's objective is to identify changes not only in the mechanism itself but also in the causal structure across the environments, which has been evaluated in artificial and real-world data sets.", "questions": "Overall, the paper is well written and easy to follow. The examples are very helpful. Here are some comments and questions:\n\n\n* While the overview of related work in the introduction is excellent, the actual related work section seems a bit thin. Consider adding literature such as “Conditional Distributional Treatment Effect with Kernel Conditional Mean Embeddings and U-Statistic Regression” by Park et al., and “Why did the distribution change?” by Budhathoki et al. Given that you build on previous work for inferring the graph structure, you could also frame this problem as comparing it node-wise, as suggested in these papers.\n* Line 112 seems to suggest that you do not support hard interventions. However, based on Eq. (3), it should be fine. Perhaps, you could clarify this more explicitly.\n* Assumption A seems to be a mix of the ‘independence of mechanism’ assumption and the ‘modularity/autonomy’ assumption. The former states that the input distribution has no information about the causal mechanism, i.e., the mechanism is independent of its input (see e.g. “Information-geometric approach to inferring causal directions” by Janzing et al. for a more intuitive explanation). The ‘modularity/autonomy’ assumption states that intervening on one mechanism does not change other mechanisms in the system. I assume you refer to the latter one.\n* In Remark 1, note that \\tilde PA_j = PA_j from related literature could also include a node that became irrelevant post-interventions. As you mentioned below, an intervention can change the functional form. For instance, in a linear model, the coefficient can become 0 of a parent, i.e., the parent can remain in the set, it is just not a ‘minimal’ set anymore. Maybe you can emphasize that you aim for ‘causal minimality’ and then having only a subset of parents after the intervention would make sense.\n* Assumption B is a quite strong, especially seeing that two environments might only change because of changes in external factors (the noise changes or becomes anomalous). Here, you would exclude these cases. In that sense, you restrict the shifts in a causal mechanism to changes in the functional form, which is significantly less general. Would be helpful to add a brief discussion about this and also emphasize the type of shift you are interested in more clearly.\n* You introduce the probability mass w (line 164), but wouldn't we know with a probability of 1 from which environment an observation comes from?\n* Stating that you do not rely on the knowledge about the causal structure seems slightly misleading as you still utilize Rolland et al.'s method. Hence, one could also use other methods that depend on knowing the DAG after applying the same algorithm. Maybe add a brief discussion regarding this.\n* Theorem 1 is a good result. Adding some intuitive explanation as to why the variance would be greater than 0 if the node shifted would be helpful.\n* The variable selection step could be more explicitly referred to as an edge pruning problem.\n* For the variable selection, would one use all data points from all environments?\n* Algorithm 2 seems slightly out of place. It should suffice to refer to that paper.\n* The experiment section is quite brief, especially the comparison with related methods. A possible approach could be to use Rolland et al.'s method to infer the graph and then apply methods as in the previously mentioned works by Park et al. and Budhathoki et al.\n\n--Update after rebuttal-- I have read the rebuttal and further discussed with the authors. Most of my concerns have been addressed and I raised the score by 1.", "limitations": "No concerns about societal impact were raised. The authors have reasonably discussed the technical limitations.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* Excellent introduction and illustrative examples to familiarize readers without strong backgrounds in this field with the problem setting.\n* Notation is concisely introduced, and the assumptions are detailed thoroughly.\n* The work is a theoretically valuable extension of existing work.\n* The method is easily applicable, and the details are sufficiently outlined for implementation.\n* The complexities are discussed in the paper.\n\nPlease see the questions section for additional points.", "weaknesses": "* The paper primarily builds on existing work, but it extends them in a sufficiently novel manner. One of my concern is the claim of not requiring knowledge about the causal structure, which is solely due to the utilization of Rolland et al.‘s existing work that aims at identifying the structure. Thus, one still relies on this; the discovery is merely part of the initial step.\n* The experiments seem quite limited, particularly with too few baseline methods for comparison (see the questions section for some suggestions).\n* The assumption about invariant noise distributions across environments is rather strong, although the other assumptions seem reasonably fair.\n* The overall complexity is relatively high.\n* The paper lacks a conclusion section.\n\nPlease see the questions section for additional points.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable feedback! We next respond to your questions, and please do not hesitate to ask follow-up questions.\n\n### Response to comments from Weaknesses section\n\nNext we only respond to the two points that are not repeated in the Questions section.\n\n> \"The overall complexity is relatively high.\"\n\nWe would like to remind two points: \n\n1. The computational complexity depends on the choice of score estimator. In this work, we used a kernel-based estimator and, as a result, the time complexity of Algorithm 1 is in the same order as that of Rolland et al. (2022). Here, any progress on score-estimation techniques will immediately improve our algorithm.\n\n2. Our approach does not make any assumption on the structural properties of the individual DAGs. Here, note that algorithms that aim to estimate only the structure of each individual DAG can have high computational costs without structural assumptions (e.g., sparsity, or small Markov boundary). For example, the PC algorithm runs in exponential time for dense graphs (Kalisch and Bühlman, 2007); for linear SEMs, Ghoshal & Honorio, (2018) developed an efficient peeling algorithm that runs on $\\mathcal{O}(d^5)$ for dense graphs; in recent gradient-based methods, the iteration complexity is $\\mathcal{O}(d^3)$ (Zheng et al., 2018) and the total time complexity is unknown. All these methods become prohibited for a large number of nodes. \n\n\n> \"The paper lacks a conclusion section.\"\n\nThank you for bringing this up. We will add a conclusion section for the camera ready.\n\n### Response to comments from Questions section\n\n> \"While the overview of related work ...\"\n\nThank you for suggesting these references. We will add and briefly discuss them in the revision. Regarding empirical comparisons against these approaches, we shall note that the work of Budhathoki et al. does not explicitly state which method they use for estimating the conditional distributions nor how they test for differences, moreover, they do not provide a source code. Similarly, there is no code available for the work of Park et al. and it is unclear how to reproduce them. We will be happy to provide numerical comparisons if more pointers are provided.\n\n\n> \"Line 112 seems to suggest that...\"\n\nThank you for noting this! To be precise, our results also hold for stochastic hard interventions. That is, an intervention can completely remove the parent set of a node (i.e. remove all incoming edges), but we still need to consider noise to have a well-defined ANM. We will add this explicitly.\n\n\n> \"Assumption A seems to be a mix of...\"\n\nGreat point! To be crystal clear, we assume that an environment admits a factorization of the form of eq.(3). As you correctly point out, we only need to consider modularity. We will fix this in the revision.\n\n> \"In Remark 1, note that...\"\n\nNote that Assumption C, which states that each functional $f_j^h(\\mathrm{PA}_j^h)$ is nonlinear in every component, implicitly ensures that no irrelevant parents are included in the set $\\mathrm{PA}_j^h$. Nonetheless, we will make a short footnote to discuss this for further clarity.\n\n> \"Assumption B is a quite strong, ...\"\n\nSee our response to all reviewers.\n\n> \"You introduce the probability mass w...\"\n\nCorrect, we know the environment membership of each sample. However, the introduction of $w_h$ is primarily to illustrate that the pooled data comes from a mixture distribution. As seen from Theorem 1, our identifiability result does not depend on $w_h$ and is helpful only for proving the theorem.\n\n> \"Stating that you do not rely on the knowledge about the causal structure...\" \n\nThroughout the text we carefully used the term “entire DAG structure” or “full structure” to explicitly indicate that our method does not require estimation of directed edges in the individual DAGs, which we argue is not misleading. As stated in our paper, we only use Rolland et al.’s leaf-identifiability result (Proposition 1 in our paper, or Lemma 1, item (i), in Rolland et al.), which has two important consequences:\n\n1. Leaf identifiability only requires estimation of the *diagonal of the score’s Jacobian*, whereas estimation of the edges in the individual DAGs would require estimation of the whole score’s Jacobian matrix. Thus, our approach is much more efficient.\n\n2. Even when a topological order is known, in finite samples, *accurate* estimation of all directed edges is very challenging for nonparametric models. For instance, in Table 3 of Rolland et al.’s paper, we can observe that SCORE, CAM and GraN-DAG, **all output a DAG with SHD of at least 130!** This means 130 edges are either missing, incorrectly oriented, or should not exist. Moreover, this high SHD value was obtained for Erdos-Renyi graphs with 50 nodes and 200 edges in expectation (an edge density of about 16%, i.e., *the graph is not dense at all*). One of the strengths of our algorithm is precisely that we do not need to estimate the edges on the individual DAGs in order to identify shifted nodes.\n\n> \"Theorem 1 is a good result...\" \n\nThank you for this suggestion! We will make sure to add a short high-level explanation of the theorem for the camera ready.\n\n\n> \"For the variable selection, would one use all data points from all environments?\"\n\nFor estimating structural changes across environments, we use data from each environment separately to estimate the parents of the *shifted nodes* and, thus, observe changes in the parent sets of the shifted nodes. We should note that estimating structural changes of shifted nodes is an optional step, it is possible that a scientist/practitioner is interested only in detecting shifted nodes, i.e., variables that have undergone a mechanism shift .\n\n\n> \"Algorithm 2 seems slightly out of place...\"\n\nFair point. For completeness, we will keep Algorithm 2 and move it to the appendix \n\n\n> \"The experiment section is quite brief...\"\n\nWe have run additional experiments including other baselines. Please see our rebuttal to all reviewers."}, {"review_id": "HSjUyCtu7f", "reviewer": "Reviewer_YGJy", "summary": "This work addresses the problem of estimating differences in the causal models in two (or more) environments, over a single set of variables with the same causal ordering. The authors are interested in discovering those variables that undergo a functional or structural mechanism shift. As their main insight, they show that under additive noise models, the Jacobian of the score function of the data distribution provides information about such shifts. With this, they follow up on a line of research that uses this Jacobian to identify leaf nodes in a causal DAG (Rolland et al. 22, Montagna et al. 23). \n\nThe authors propose the MSG algorithm, which given data from two environments, first discovers nodes which undergo a mechanism shifts, and in a second stage discovers the causal parents of shifted nodes as well as the edge differences between the causal DAGs in each environment.  Empirically, they evaluate MSG on  discovering nodes with mechanism shifts, and show that MSG outperforms the DCI approach (Wang et al. 18)  while scaling linearly in the number of variables.", "questions": "- What is the best way to extend MSG to more than two environments? For example, is it still possible to estimate the score from the mixture distribution over multiple contexts, and if so, how does this depend on the number of mechanism shifts over all contexts?  \n- Definition 3 seems to admit the case where the causal parent set is empty in a given environment, that is, where a hard intervention occurs. Given that this case is not considered in the paper (footnote 3), what problems does it cause for MSD? Regardless, it would be interesting to consider hard intervention in the experiments.\n- The authors mention the approach by Ghoshal et al. 21, which is more recent than and slightly outperforms DCI. Is there a reason that the experiments do not include it?  \n- Why is the evaluation result for 10 nodes sometimes missing in Appendix D (Fig. 10, 11)? I would also be interested in the case of 10 nodes in the main experiment in  Fig. 3, since it seems that MSD has lower precision on smaller graphs (e.g. 20 nodes).\n- Fig. 3 seems to suggest that shifted node discovery improves with increasing number of nodes d, whereas in Fig. 1 it degrades, why is this? \n- In the experiments on difference DAG discovery in Fig. 11,  MSG has F1 scores of at most 0.5. What most likely causes this? For example, is it related to the estimation of the causal parent sets with Alg. 2? While Theorem 2 establishes the consistency of discovering the correct parent sets, is this also true in practice? \n- Can the hypothesis tests in Appendix E be extended to also accomodate changes in noise distributions?  \n- *smaller comments*: Citation style in Section 2.2 (numbers as subjects), Precision and Recall missing in Fig. 11 (while included in Fig. 13 which has the same setting).\n\n", "limitations": "- *Types of Mechanism Shifts*: The setting is motivated in the context of both functional and structural differences in SCMs (Def. 3, 4). While Alg. 1 for discovering shifted nodes can deal with both of these shifts, the evaluation focuses only on a specific type of *functional* shifts. Alg. 3 for discovering the causal parents, on the other hand, is only evaluated for *structural* shifts (Appendix D). The authors  discuss the theoretical aspects of using Alg. 3 under functional shifts in Appendix E, but it would be interesting to also see an evaluation.  \n- *Experiments*: As elaborated under weaknesses, the evaluation is limited because the authors only consider the linear competitor DCI, as well as do not evaluate parts of the proposed algorithm (Alg. 2, Alg. 5) or only in specific settings (Alg. 1, Alg. 3).\n", "rating": 5, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "- *Presentation*: The authors motivate the problem well and place the approach appropriately into the context of existing work. The paper has a clear structure and is easy to follow.\n- *Soundness*: The authors state the theoretical results clearly and provide detailed proofs for each claim.\n- *Scalability*: As the authors point out, their approach can identify mechanism shifts while avoiding a computationally costly estimation of the complete causal DAG. In addition, the proposed algorithm scales linearly in the number of DAG nodes, contrary to most existing work. Lastly, they propose to estimate the score function using the mixture distribution over multiple environments, suggesting an improved sample efficiency compared to approaches that use per-environment distributions.\n- *Assumptions*: MSG does not make any structural assumptions on the DAGs, nor restrictive  functional modelling assumptions.\n- *Contribution*: Existing work on discovering difference DAGs assumes linear  functional models (Wang et al. 18, Ghoshal et al. 21). Meanwhile, nonlinear approaches that can indirectly solve this problem by discovering a causal DAG over multiple contexts and identifying mechanism changes for each node (e.g. Mooij et al. 16, Huang et al. 19) do not scale well in practice. Hence, this work is a valuable addition to the literature. \n", "weaknesses": "- *Novelty*: This work follows a recent line of research which also uses the Jacobian of the score of the data distribution to identify leaves in a causal DAG  (SCORE, Rolland et al. 22, Montagna et al. 23). The proposed approach applies similar techniques and algorithms, making the approach incremental in nature. \n- *Comparison to SCORE*: To demonstrate the practical relevance of their insights and the novelty compared to the SCORE framework, it would be helpful to see a comparison to SCORE in the experiments. For example, is there a noticeable benefit of MSD compared to applying SCORE to each environment and comparing the resulting DAGs for structural differences?  \n- *Comparison to DCI*: Given that DCI makes a linearity assumption, the superior performance of MSD in the experiments might only be due to DCI's model misspecification in the nonlinear setting. It might be interesting to include the linear case, where MSD's functional model is misspecified, or to consider additional competitors.\n- *Competitors:* There are many approaches that solve this problem indirectly, by (a) discovering a joint DAG over multiple environments and functional shifts for each node, or (b) discovering a DAG per environment, which we can compare for  structural differences. For discovering shifted nodes, MSG has a computational advantage of not needing a direct estimation of DAGs, but it would be interesting to see whether it also has an advantage in terms of accuracy.  For discovering structural differences, MSG has an F1 score of around 0.5 (Fig. 11), thus it would be informative to see whether direct DAG estimation improves upon this.\n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your valuable input. We have new experiments on additional settings, and other baselines, as per your suggestions. Please see response to all reviewers.\n\n### Response to comments from Weaknesses section\n\n> \"Novelty: This work follows ...\"\n\nWe firmly disagree with the statement that “our approach is incremental in nature”. Here we provide some points to take into account.\n\n1. *The problem itself is very different to that of Rolland et al.*, and coming up with a more efficient algorithm that is different from the naive approach (i.e., estimating individual DAGs, fitting functionals, and testing for differences) is far from trivial. The effectiveness of our approach is demonstrated through empirical experiments, and it is grounded on solid theory.\n\n2. The proof of Theorem 1 is very different and more convoluted than the proof of Proposition 1 from Rolland et al.\n\n3. There is a consensus from the reviewers that the identifiability result of Thm 1 is easy to digest. The latter does not necessarily imply a lack of novelty or creativity; on the contrary, it often means that the method is more accessible for others to understand and apply.\n\nWe sincerely hope that you will reconsider our method's novelty. We believe a simple yet effective approach should be celebrated, not criticized. \n\n\n> On additional comparisons\n\nWe tested additional settings such as functions drawn from Gaussian processes, and other baselines based on indirect methods. Please see response to all reviewers.\n\n### Response to comments from Questions section\n\n> \"What is the best way to extend MSD to more than two environments?...\" \n\n* The problem setting, as stated in L141, considers data from $H$ environments. All the theoretical statements, including proofs, are also given in full generality, i.e., considering $H$ environments. For the experiments, we mainly consider $H=2$ because DCI only works for two datasets.\n\n* Regarding estimation, in L5 of Algorithm 1 (or L9 of Alg 4), we estimate the diagonal of the score’s Jacobian of the mixture distribution, which consists of data gathered from all the environments. This estimation process does not depend on the number of mechanism shifts. \n\n\n> \"Def 3 seems to admit ...\"\n\nGood catch! We will assume that you meant Def 2, as Def 3 is related to shifted nodes. Correct, our results also support the case where the parents of a variable are completely removed. Formally, our identifiability result also holds for stochastic hard interventions, i.e., the parents are removed but we keep the noise to have a well-defined ANM for each environment.\n\n\n> \"...authors mention the approach by Ghoshal et al...\"\n\nThe reason was primarily because DCI's code was accessible, and Ghoshal et al.'s was not. However, we should point out that, as illustrated in Fig 3 of Ghoshal et al., the difference in performance is not significant. Given that both methods study linear models, we believe that DCI was sufficient as a baseline. \n\n> \"Why is the evaluation result for...\"\n\nGood catch! The absence is due to not being able to generate graphs with the given specification, in this case, ER6/SF6 graphs. For clarity, at $d=10$, an ER6/SF6 graph corresponds to a graph with 60 edges in expectation, which is impossible to generate since there are at most 45 edges in a DAG of 10 nodes. The reason there exist a point for $d=10$ in SF6 graphs (e.g. in Figs 10, 11) is because the package we used to generate SF graphs automatically sets to the maximum number of possible edges if larger numbers are given; whereas the package to generate ER graphs, does not. We will remove all points for ER6/SF6 graphs at $d=10$.\n\n> \"Fig. 3 seems to suggest...\"\n\nThere should be a typo in your question since Fig 1 does not show any empirical results, thus, we will assume that you refer to Fig 6, not Fig 1.  A high-level explanation can be given by looking at the data generation process in L262-L268. In this setting, the structures in the two environments are exactly the same, and only the functional mechanisms of the shifted nodes change. Then, for a given number of nodes $d$, the mechanisms of shifted nodes in ER6/SF6 graphs will observe higher “fluctuations” since the expected number of parents is higher than that of ER2/SF2 graphs, which our method successfully captures.\n\n\n> \"In the experiments on difference DAG...\"\n\nGood question. To estimate a difference DAG (Alg 3) we use: Alg 1 to detect the shifted nodes, and Alg 2 for estimating local parents. While both Alg 1 and Alg 2 affect the performance of estimating the difference DAG, we attribute your observation primarily to Alg 2. While Thm 2 provides the consistency result of Alg 2, it may not always translate to strong performance in finite samples. Finally, it is clear that incorrect estimates of shifted nodes from Alg 1 will also affect the final difference DAG accuracy, but we believe Alg 1 is not the main source of errors as it performs quite well as seen in the experiments. \n\n\n> \"Can the hypothesis tests...\"\n\nSee response to all reviewers. We can drop Assumption B, i.e., our method will also identify changes in noise distributions.\n\n### Response to comments from Limitations section\n\n> \"Types of Mechanism Shifts...\"\n\nOur paper's primary goal is to identify mechanism shifts without specific functional or structural assumptions, as stated in our contributions (L63-L70). While our method can identify structural and functional shifts in nodes, estimating shifted edges may require an additional set of assumptions. In App E we explored a potential semi-parametric assumption to identify functional shifted edges, but this is not central to our work, and was provided mainly to motivate future directions.\n\n> \"Experiments...\" \n\n* See new experiments in the response to all reviewers.\n* We executed Alg 3 (and thus Alg 2) in the process of recovering the difference DAG, as detailed in Section D.2. Given that App E is not central to the goal of this paper, we decided not to evaluate Alg 5.\n"}, {"review_id": "eUsXLx55MF", "reviewer": "Reviewer_PS2y", "summary": "The authors propose a two-step algorithm for identifying causal mechanism shifts. The algorithm is motivated by a characterization of the nodes with mechanism shifts and the procedure is based on several existing tools. The proposed algorithm outperforms the baseline.", "questions": "Regarding the statistics defined in (6) from Appendix, if the leaf node L is not shifted, the statistics may not be small since both the numerator and denominator are approximately zeros. Did I miss something here?\n\n\n\n\n", "limitations": "No potential negative societal impact.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. Theorem 1 provides a new perspective for the characterization of nodes with shifted mechanisms.\n\n2. The combination of existing tools in Algorithm 2 is interesting. \n\n3. The computation complexity of the algorithm is reported. ", "weaknesses": "1. Assumption B (invariant noise distribution) makes the setting quite restrictive. There is no discussion or experiments about how sensitive the proposed method is w.r.t. shifted noise distributions, or whether the proposed idea can be generalized to handle shifted noise distributions.\n\n2. I suspect that the statistic is not well-defined in some cases (see Questions below). \n\n3. The threshold $t$ is a key hyper-parameter for the method, but the suggested elbow strategy is not a principled approach. \n\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments! We next respond to your concerns, and hope our answers will make you feel more positive about our work. Please do not hesitate to ask any follow-up questions.\n\n> Assumption B (invariant noise distribution) makes the setting quite restrictive. There is no discussion or experiments about how sensitive the proposed method is w.r.t. shifted noise distributions, or whether the proposed idea can be generalized to handle shifted noise distributions.\n\nWe have excellent news here, we can drop Assumption B and still have our method to identify distribution shifts without major changes to our proof of Theorem 1! Please see our response to all reviewers. Given that your main concern is related to this assumption, we hope you will feel more positive about our contributions.\n> The threshold t is a key hyper-parameter for the method, but the suggested elbow strategy is not a principled approach. \n\nWe agree but note that we never claimed that the elbow method is a principled approach, in fact, it is well known to be a classical heuristic for estimating the number of clusters in a dataset. In this case, as the experiments in Appendix D.3 suggest, it can lead to good performances.  Finally, we also believe that a principled approach would be nice to have, which makes this an exciting future direction!\n\n\n> Regarding the statistics defined in (6) from Appendix, if the leaf node L is not shifted, the statistics may not be small since both the numerator and denominator are approximately zeros. Did I miss something here? \n\nGreat question! If the score estimator is well behaved (meaning that as the number of samples increases, the estimated value gets closer to its population value), then, for a non-shifted leaf node $L$, the numerator would converge faster towards zero. This is due to the numerator using all the pooled data for the estimation of the score; on the other hand, the denominator uses only data from a single environment. Thus, when the leaf node $L$ is not shifted, we can expect the statistic $\\mathtt{stats_L}$ to be lower than 1.\n"}, {"review_id": "JmNnwQ0sdG", "reviewer": "Reviewer_appw", "summary": "This paper focuses on identifying mechanism shifts in related structural causal models without learning the structures of the individual SCMs. It builds on the recent advances in score-matching based causal discovery and consider nonlinear additive noise models, unlike the prior work that assumed linear Gaussian models. The proposed algorithm first identifies the shifted variables and then estimates the structural changes associated with each shifted variable.", "questions": "- After finding the shifted nodes, finding the parents should not be a big deal besides the computational burden. One angle can be your efficiency of Algorithm 2. Can you elaborate on its efficiency and which parts are your contribution?\n- Invariant noise (Assumption B): How do you deal with the root nodes? Do you simply assume that root nodes are never intervened? I see that this is how simulations are conducted. It should be clearly stated in the paper as well, as part of the assumptions perhaps. \n- If I recall correctly, DCI algorithm of Wang et al. [77] does not require invariant noise assumption. In contrast, the algorithm of Ghoshal et al. (2019) [23] relies on invariant noise and offers some advantages over DCI. I believe the algorithm of Ghoshal et al. can be an important comparison too.\n- Speaking of related work and comparisons, I believe the present work misses an important reference, “Scalable intervention target estimation in linear models” Varici el al. (2021). Although it’s focused on just identifying the intervened (shifted) nodes, their algorithm also finds some edges. For the comparison of Algorithm 1, along with Ghoshal et al. (2019), Varici et al. (2021) should also be a valid comparison in addition to the already used DCI. The formulation in Varici et al. (2021) differs by assuming that the noise variables necessarily change, but I believe a change in the linear weights is enough to relax the “varying noise” assumption. (you can check their follow-up Varici et al. (2022) in which they discuss the same problem for causally insufficient models with multiple environments)\n- Finally, all these three papers (Wang et al., Ghoshal et al., Varici et al.) are based on linear models whereas the present paper requires non-linear models. More complex algorithms that can work with non-parametric causal models can be used. For instance, UT-IGSP algorithm of Squires et al. (2020), JCI algorithm of Mooij et al. (2020), and $\\psi$-FCI algorithm of Jaber et al. (2020) all can return intervention targets along with the causal graph. The latter two can be computationally prohibitive, but I think UT-IGSP can work on fairly large graphs as well for a fair comparison of identifying shifted nodes.\n- Including code in supplementary could have been helpful.\n\n Varici, B., Shanmugam, K., Sattigeri, P., & Tajer, A. (2021). Scalable intervention target estimation in linear models. Advances in Neural Information Processing Systems, 34, 1494-1505.\n\nVarici, B., Shanmugam, K., Sattigeri, P., & Tajer, A. (2022, August). Intervention target estimation in the presence of latent variables. In Uncertainty in Artificial Intelligence (pp. 2013-2023). PMLR.\n\nSquires, C., Wang, Y., & Uhler, C. (2020, August). Permutation-based causal structure learning with unknown intervention targets. In Conference on Uncertainty in Artificial Intelligence (pp. 1039-1048). PMLR.\n\nMooij, J. M., Magliacane, S., & Claassen, T. (2020). Joint causal inference from multiple contexts. The Journal of Machine Learning Research, 21(1), 3919-4026.\n\nJaber, A., Kocaoglu, M., Shanmugam, K., & Bareinboim, E. (2020). Causal discovery from soft interventions with unknown targets: Characterization and learning. Advances in neural information processing systems, 33, 9551-9561.\n\n**Minor notes**:\n\n- Montagna et al. “Scalable Causal Discovery with Score Matching”: This is cited in the appendix. I believe it deserves more discussion in the main paper as it has some improvements over Rolland et al. (2020) especially for parent selection and efficiency of their algorithm. I see that the recent version became public at arXiv just shortly before the submission of the present paper, though an earlier version was available (NeurIPS 2022 Workshop on Score-Based Methods).\n- *ICM assumption*: Do you really need ICM assumption? It seems to me that you only need modularity principle of Pearl, that is each causal mechanism (conditional prob. dist.) is independent of the others. I see that you use it to say that “there are only a few shifted nodes” but if I am not mistaken, your analysis does not depend on such sparsity assumption.\n-  [62] (Saeed et al. 2020): To my knowledge, their result for identifying shifted variables comes under poset compatibility assumption that can be quite restrictive for the structural changes.\n", "limitations": "The main limitations of the present work would be their assumptions. They are not discussed in a separate paragraph or pointed out somewhere in the paper (if I am not missing). But they are clearly written and can be understood at where they appear in the main text. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- I find the main problem — identifying the changes without learning the graphs in their entirety — important.\n- Score-based (gradient of log likelihood) are gaining momentum in causal inference and this can be a timely paper to expand the existing recent work. Also, I think soft interventions are also relatively understudied but are important for modeling different environments with sparse changes.\n- Paper is written clearly, the flow is natural and easy to follow.", "weaknesses": "- *Invariant noise assumption (Assumption B) and shifted node definition*: I see that invariant noise assumption is crucial for your results, especially it is used in Lemma 2 (and subsequently for the proof of Theorem 1). However, I find this a somewhat unrealistic assumption. For instance, your current definition of shifted nodes excludes root nodes — they can be intervened only upon a change in the noise distribution. I believe Assumption B needs a convincing discussion. I add more comments on questions below.\n- The main contribution (Theorem 1) is rather a straightforward extension of the existing result (Proposition 1). Sure, there may be technical challenges, but I don’t see a fundamental difficulty. If I am missing anything, please let me know. \n- Although the paper provides extensive experiments in numerous settings, it misses comparisons with some important related work.\n\n========= AFTER AUTHORS' REBUTTAL ==========\nA major concern, the invariant noise assumption, is addressed by the reviewers in their rebuttal. Hence, I updated my score accordingly. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable feedback! We next respond to your questions, and hope our answers will make you feel more positive about our work. Please do not hesitate to ask any follow-up questions.\n\n### Response to comments from Weaknesses section\n\n> \"Invariant noise assumption (Assumption B) and shifted node definition...\"\n\nWe have great news regarding Assumption B, we can relax it and allow for changes in the noise distribution! Our method is now capable of identifying mechanism shifts of a node, stemming from either changes in the functional relationship w.r.t. its parents or shifts in the noise distribution. Please see our rebuttal to all reviewers for more details.\n\n> \"The main contribution (Theorem 1) is rather a straightforward extension of the existing result (Proposition 1). Sure, there may be technical challenges, but I don’t see a fundamental difficulty. If I am missing anything, please let me know.\"\n\nWe respectfully disagree with your comment, but we earnestly hope to clarify our perspective and appeal to your understanding.\nOur work's novelty can be seen from two angles: the question we address and the method we propose.\n\n* The Question: As you kindly acknowledge in the strengths, the problem we tackle is both interesting and important. To the best of our knowledge, this paper is the first to consider this specific question, a natural extension motivated by previous work such as DCI. We genuinely believe that coming up with a more efficient algorithm that is different from the naive approach (i.e., estimating individual DAGs, fitting functionals, and testing for differences) is already a novel contribution. \n\n* The Method: First, our proposed method is innovative in that it bypasses any parametric form assumption on the functional mechanisms, and does not rely on any structural assumption in the individual DAGs, as pointed out in our contributions (L63-L70). The effectiveness of our approach is demonstrated through empirical experiments, and it is grounded on solid theory. Second, the proof of Theorem 1 is very different and more convoluted than the proof of Proposition 1 from Rolland et al. Finally, Theorem 1 is easy to digest and our algorithm is simple to implement, making our paper more accessible for others to understand and apply. We sincerely hope you can see that the simplicity of our approach should not be mistaken for a lack of innovation.\n\n\n### Response to comments from Questions section\n\n> \"After finding the shifted nodes, finding the parents should not be a big deal besides the computational burden. One angle can be your efficiency of Algorithm 2. Can you elaborate on its efficiency and which parts are your contribution?\"\n\nWe argue that finding the parents of a node (a.k.a. Local causal discovery) is also a very challenging task, not just computationally but also methodologically; see for instance Gao & Aragam (2021) “Efficient Bayesian network structure learning via local Markov boundary search”. Our approach addresses this task by using an estimated topological ordering returned by Alg 1, and adapting recent methods like FOCI, which efficiently measures the level of dependence between random variables. We've shown FOCI's efficiency empirically and its consistency for local recovery in Theorem 2. \n\n\n> \"If I recall correctly, DCI algorithm of Wang et al. [77] does not require invariant noise assumption...\"\n\nTo be precise *both DCI and Ghoshal et al.’s method assume invariant noises in order to identify all edge orientations.* Without such assumption both methods can only output a partially oriented difference DAG. See Theorem 4.4 in Wang et al.’s DCI paper, and Theorem 1 in Ghoshal et al.’s paper.\nRegarding a comparison against Ghoshal et al.'s method. We should point out that, as illustrated in Fig 3 of their paper, the difference in performance against DCI is not significant. In addition, their code is not available. Given that both DCI and Ghoshal et al.'s study linear models, we believe that DCI is sufficient as a baseline. \n\n\n> \"Speaking of related work and comparisons, I believe the present work misses an important reference...\"\n\nThank you for suggesting these references! We will make sure to include them in the revised version. Additionally, we have incorporated some additional baselines in our experiments, as detailed in our rebuttal to all reviewers. \nIt is also worth noting that *our problem setting is not the same as identifying intervened nodes.* Please see our response to all reviewers.\n\n\n> \"Finally, all these three papers (Wang et al., Ghoshal et al., Varici et al.) are based on...\"\n\nThank you again for sharing additional related work. We have included UT-IGSP as a baseline as well, see our response to all reviewers.\n\n\n> Including code in supplementary could have been helpful.\n\nWe will make sure to release code for the camera ready.\n\n\n> \"Montagna et al. “Scalable Causal Discovery with Score Matching”: This is cited in the appendix...\"\n\nWe will certainly discuss this work in our revised version. However, from our understanding, Montagna et al.’s work does not improve on the efficiency, nor in parent selection. It cannot be more efficient as they also perform a score estimation step as in Rolland et al.’s, and the edge pruning step (or parent selection) is based on CAM’s pruning procedure. On the other hand, Montagna et al.’s work does support arbitrary noise distributions, which makes for intriguing future work to investigate how to use some of their ideas for our problem setting.\n\n\n> \"ICM assumption...\"\n\nThanks for pointing this out. Indeed, we only need the modularity principle, we will correct this in the revision.\n"}]

Author rebuttal:
We thank all reviewers for their time and effort put into our work. We truly appreciate your feedback and we are happy to see that your comments will lead to a stronger version of our manuscript. Overall, there is a positive feeling about the significance of the problem, our contributions, and presentation of our work. We detected two major concerns from most of the reviewers that apparently lead them to assign a borderline score: One concern is on Assumption B, and the other on the lack of additional baselines in the experiments. We next tackle both, and additionally highlight important aspects of the problem setting and our contribution that may have been overlooked. In light of the rebuttal process, we hope our answers will successfully address your concerns, and we welcome any follow-up or additional questions you may have.

### 1. On Assumption B.

TL;DR: we can relax Assumption B and allow for changes in the noise distributions, without many changes in the proof. The definition of shifted node will be updated to changes in the mechanisms (conditional distribution), and not just *functional changes*, as the mechanism shift can now stem from changes in the functional or in the noise distribution.

Nevertheless, we would like to point out that Assumption B was not unique to our work and was actually considered in the two most closely related works to ours, namely, by Wang et al. and by Ghoshal et al. The consistency result in Wang et al.’s Thm 4.4 requires invariant noises in order to orient all edges in the difference DAG. Likewise, Ghoshal et al. also assume invariant noises in their Thm 1. Motivated by these works, we almost by default assumed invariant noises and focused on identifying functional shifts. 
Thanks to your criticisms on this assumption, we were challenged to think about ways to relax it, and ended up proving that without Assumption B, our method can more generally identify distribution shifts (not only functional shifts). We briefly mention the updates to be made.

Updates: Without Assumption B, Thm 1 (line 187) can be modified as follows:

* if $j$ is a leaf in all DAGs, then the conditional distribution $p^h(X_j^h \mid PA_j^h)$ is shifted if and only if $\mathrm{Var}_X[\frac{\partial s_j(x)}{\partial x_j}] > 0$.

The proof updates start at eq.(8) (L589), where the last two summation terms can only be canceled out if and only if:
 $ \frac{\partial}{\partial x_j} \log p^h_{N_j}(x_j-f_j^h(PA_j^h))=\frac{\partial}{\partial x_j}\log p^{h'}_{N_j}(x_j-f_j^{h'}(PA_j^{h'}))$ for all $h,h'\in [H].$

For brevity, consider Gaussian noises, then the above equation boils down to:
$ \frac{x_j-f_j^h(PA_j^h)}{(\sigma_j^h)^2}=\frac{x_j-f_j^{h'}(PA_j^{h'})}{(\sigma_j^{h'})^2}$, for all $x_j, PA_j^h, PA_j^{h'}.$ 
With some algebra, we can prove that this equation holds if and only if $\sigma_j^h=\sigma_j^{h'}, PA_j^h = PA_j^{h'}$, and $f_j^h(PA_j^h) = f_j^{h'}(PA_j^{h'}).$
Since the conditional distribution is $p^h(X_j\mid PA_j^h)\sim N(f_j^h(PA_j^h),(\sigma_j^h)^2)$, if both the mean and variance parts are the same, there are no conditional distribution shifts, and the last two summation terms in eq.(8) will be canceled out, resulting in zero variance. Otherwise, any difference in functional relationships or noise variance across environments will prevent the cancellation of the last two summation terms in eq.(8), inducing positive variance. Lastly, Lemma 2 will also be updated.

### 2. Comparison to other baselines

We include plots of our new experiments in the attached PDF. Due to space constraints, we only show F1 scores and Runtimes, but will provide all details in the revision.

**New baselines:** As per your suggestions, we ran more experiments, similar to the setting described in App. D.2, by adding three more baselines: SCORE (Rolland et al.), UT-IGSP (Permutation-Based Causal Structure Learning with Unknown Intervention Targets), and CITE (Scalable Intervention Target Estimation in Linear Models). Figures R.1 and R.2 correspond to results on shifted nodes and structural differences, respectively.

**Experiments using Gaussian Processes (GPs):** We also include experiments where the functionals are drawn from GPs. These results are in Figure R.3. We will include all experimental details in the revision, but will be happy to clarify any questions in the discussion period. 

As seen in the figures, our method consistently outperforms the other baselines in terms of the F1 score, further affirming the robustness and efficiency of our approach.

### 3. Our problem setting is not the same as identifying intervention targets

While some reviewers have suggested that methods for estimating intervention targets could be applied to our problem, we would like to clarify that there is a fundamental difference between "estimating mechanism shifts" and "estimating intervention targets." Briefly, **all shifted nodes are intervened nodes, but not all intervened nodes are shifted nodes.** This is because we are interested in *changes across environments*. For example, it is possible that a variable has the same intervention across the environments, and thus, it would not be considered a shifted node. Finally, we do not assume access to observational data, which is a typical assumption for methods designed to detect intervention targets.

### 4. Important assumptions that we do not make 

With the exception of Reviewer YGJy, we note that two important aspects of our approach have been largely overlooked, and we would like to reiterate their importance: (1) We do not make parametric assumptions on the functional mechanisms; and (2) We do not make any structural assumptions on the individual DAGs, such as sparsity or small Markov boundaries. Without these assumptions, the naive approach of learning the individual DAGs, fitting functionals, and then testing for changes, is doomed to fail. Not only this naive approach becomes computationally intractable, but also is statistically inefficient.


Meta-review: The paper provides a clear to read framework, makes validated algorithmic contributions and is a pleasure to read. All reviewers are in favor of acceptance. While there were considerable (positive) discussions with the reviewers and it is highly recommended and assumed that the authors include the take-aways and improvements as well as the additional results into the final version of the paper.
Decision: Accept (poster)