2.8 Could you explain more about "The clustering accuracy and performance metrics show high standard deviation when K is set to the true number of mixture components K* = 10." which is stated in the appendix and explain more about the "buffers" when K>K*?

We have added a new section Appendix B.3 with additional experiments for the $D=10, K*=10$ synthetic dataset with $K=10$ and $K=20$ to illustrate our assertion. We plot both the number of unique graphs and the clustering accuracy as a function of number of training steps. We observe that when $K=20$, as training progresses, the algorithm groups together points from different causal graphs until they more or less converge to the “true” number of causal graphs and clustering accuracy converges to (approximately) 100%; however, when $K=10$, we observe that for some runs, the number of unique graphs can sometimes fall below $K^\ast=10$, and thus the clustering accuracy can suffer.

3. Other questions

3.1 Why maximizing ELBO is equivalent to maximizing or minimizing each term separately in the log-likelihood in terms of $\theta$, $\phi$, and $\psi$ as $\phi$ and $\psi$ appear together in the second term?

We did not mean that each term in the log-likelihood is maximized or minimized separately; indeed we mean the joint optimization of the ELBO in terms of $\theta, \phi$ and $\psi$ is equivalent to optimizing the log-likelihood. We have updated Section 4 in the draft to clarify this.

3.2 Is it better to mention the Rhino algorithm in the related work section and briefly explain how it works as the proposed algorithm is built on top of Rhino?

Thank you for the suggestion! we have already included a self-sufficient description of the Rhino algorithm in Appendix C due to page limit.

3.3 Is it essential to verify whether the data satisfies condition (*) to ensure the reliability of the results? If so, how to verify this?

Yes, it is essential to ensure the identifiability of the underlying mixture model in order to ensure that the problem of causal discovery is a statistically feasible problem. Condition (*) gives one sufficient condition to verify this. While condition (*) can be difficult to verify exactly in most scenarios since the log-likelihood is rarely explicitly known, in practice one can “approximately” check it given good enough approximations to the log-likelihood functions of the underlying causal models, by explicitly checking if any of the samples satisfy the condition.

3.4 The paper mentions that people can also infer functional equations through the proposed model. Could you help me locate the related outputs in the experiment section? Does this mean the causal effect?

We do not benchmark the causal effect, since our focus in this paper is on causal discovery. We will work on benchmarking the causal effect in a future iteration of the paper.

3.5 Has any code been provided?

We will release our code after the peer review process.

References: [1] Gong, Wenbo, et al. "Rhino: Deep causal temporal relationship learning with history-dependent noise." arXiv preprint arXiv:2210.14706 (2022).

We thank the reviewer for their useful feedback and questions. Below we address the main concerns:

Weaknesses:

3. The baselines used in the paper are not algorithms designed for multiple DAGs and mixture SCMs. Even though the authors mentioned other algorithms designed for multiple DAGs, none has been applied in the comparison.

We are not aware of any baselines dealing with time-series data that are designed for multiple DAGs. All the algorithms mentioned in the Related Work section are designed for independent data, and hence not directly applicable to our setting. We would appreciate any concrete pointers to comparable algorithms.

Questions

1. Notations and explanations

1.1 The SCM described in equation (1) differs from equation (6) and also the equation below equation (6); which one is correct? Do you assume additive noise?

We presented equation 1 for preliminaries, however we assume additive noise in our work due to its structural identifiability. We have clarified this in the updated draft, in Section 3.

1.2 In Theorem 1 equation (*), what is the definition of $a_i$? Does it mean one sample from $i$th SCM?

$a_i$ refers to any point in the support of the mixture distribution such that the condition (*) is satisfied (as specified in mathematical notation in the theorem statement). It is not a ‘sample’ from the ith SCM in the sense that it is not randomly sampled from the SCM, but it can be potentially explicitly chosen for the condition to hold. We have added this clarification to the draft in Appendix A.3.

1.3 In Theorem 1, what are $g_1$ and $g_2$? Should they be $h_1$ and $h_2$?

Thank you for pointing out the typo. Yes, they should be $h_1$ and $h_2$. We have corrected this in the new draft in the theorem statement.

1.4 As mentioned in the main paper, "This highlights the idea that learning a mixture model is only beneficial when the underlying SCMs differ from one another significantly. " could you clarify what the significant difference is here?

By “significant difference”, we mean that the causal graphs in the mixture distribution differ from one another significantly in terms of their edges and resulting data likelihood. This is also exemplified by the following table of pair-wise SHDs between the constituent graphs of the mixture models.

We note a notably higher clustering accuracy and improved performance in causal discovery when evaluating the Netsim-permuted dataset. This is in contrast to the Netsim dataset, where the pairwise Structural Hamming Distance (SHD) is smaller, indicating that the mixture graphs in this dataset are more similar to each other.

2. Details in the experiment

2.1 In the experiment section, "(per sample) indicates that the baseline predicts one graph per sample." how to apply the algorithm with only one sample? For example, PCMCI+ needs CI tests, and CI tests need a set of samples instead of only one sample.

Each sample is a time-series of length $T$, while PCMCI+ performs conditional independence testing over time-windows of length $L$, where $T > L$. Thus, CI tests are conducted on a set of time windows per time series.

2.2 In the experiment performance section, AUROC and F1 are used as metrics. However, it is not mentioned which kind of AUROC and F1 refer to adjacency or orientation AUROC/F1?

Thanks for asking for the clarification. All the reported metrics are orientation AUROC and F1. We have updated the Experiments section of the draft to clarify this as well.

2.3 As far as I know, the Rhino algorithm is not designed for heterogeneous time series data. In the experimental section, would it be feasible to compare the proposed algorithm with "Rhino (grouped)" as well? This would involve applying the Rhino algorithm to the grouped data based on the true underlying causal graph. Considering this, the results regarding Rhino in Figures 3 and 6 could be more comprehensive if grouped data were utilized.

You are correct that Rhino is not designed for heterogeneous time series data. Based on your suggestion, we have added Rhino (grouped) results on the synthetic datasets in Figures 3, 6 and 7. We observe that Rhino (grouped) on average performs at the same level or worse than MCD on the synthetic datasets (especially for higher number of graphs). This is due to the reduced number of samples per graph, which MCD is more robust to due to weight sharing (as implemented in Equation 10).

Thanks for your response. Below we answer your queries:

1.1. If the SCM assumes additive noise, perhaps mentioning it in the abstract or contribution section would be beneficial.

Thank you for the suggestion. We have clarified that we work with the additive noise assumption in the contributions section of the introduction.

2.1 Based on my understanding of PCMCI, the algorithm operates on a p-variate time series with length T, conducting CI tests across the entire set of time series samples rather than within a sample window. If my understanding is incorrect, could you direct me to the specific section in the PCMCI$^+$ paper discussing the time window?

Quoting from the PCMCI paper [1] "The nodes in a time series graph represent the variables at different lag times and a causal link $X^i_{t \rightarrow \tau} \rightarrow X_t^j$ exists if $X_{t-\tau}^i$ is not conditionally independent of $X_t^j$ given the past of all variables, formally defined by $X^i\_{t-\tau} \not\perp X\_t^j | (\mathbf{X}\_t^{-}- \\{ X\_{t-\tau}^i \\})$ with $\not\perp$ denoting the absence of a (conditional) independence, the vertical bar $|$ meaning 'conditional on', and $\mathbf{X}\_t^- - \\{ X_{t-\tau}^i\\}$ denoting the past of all $N$ variables up to a maximum time lag $\tau_\text{max}$ excluding $X^i_{t-\tau}$." Since the PCMCI algorithm assumes causal stationary, the CI tests are carried out for time windows of length $\tau_\text{max}$, which denote the maximum time length for which the causal relationships are inferred.

2.1' I have a follow-up question based on 2.1. Considering PCMCI$^+$ requires only one p-variate time series with length T (where the sample of time series needed is $n=1$) I'm curious about the input of PCMCI$^+$ in your experiments. Specifically, I had assumed that the input dimension for PCMCI$^+$ should be $p \times T$. How is a dataset sized $n \times p \times T$ fed into PCMCI$^+$, where $n$ signifies the number of time series samples?

PCMCI$^+$ as implemented in the Tigramite package, has two options based on whether there is a single time-series input of size $p \times T$ or multiple time-series of size $n \times p \times T$. When inferring a single graph per sample, we use the analysis_mode='single' setting. When inferring a graph per dataset we use the analysis_mode='multiple' option [2].

2.6 Do you mean that some settings (e.g., D=10) could be inherently easier than others (e.g., D=5) simply due to the stochasticity of the generation process? I understand your point; just want to ask whether it is a typo in your answer.

Yes, this is what we meant. Apologies for the typo in our response.

References: [1] PCMCI: J. Runge, P. Nowack, M. Kretschmer, S. Flaxman, D. Sejdinovic, Detecting and quantifying causal associations in large nonlinear time series datasets. Sci. Adv. 5, eaau4996 (2019). https://advances.sciencemag.org/content/5/11/eaau4996 [2] Tigramite package. https://jakobrunge.github.io/tigramite/#module-tigramite.data_processing

We express our gratitude to the reviewer for their valuable feedback and inquiries. In the following section, we will tackle the primary issues raised.

Weaknesses

If we turn to the experimental section of the paper, we get two examples: one for NetSim (which I’m very familiar with) and another for the DREAM3 gene network. Both of these have problems with respect to the goal of this paper, suggesting that the choice of experimental datasets could be improved.

We found it quite challenging to find large-scale real-world experimental datasets for time series causal discovery. If you have any suggestions or pointers for datasets, we would be glad to run experiments on them and incorporate them into our draft.

The problem with the NetSim data is that it’s not an extremely convincing time series, as the records in the simulation are spaced far enough apart in time to render the data nearly i.i.d. In fact, analyzing it as i.i.d. often yields better results than analyzing it as time series, frustratingly, as in this paper:

Ramsey, J. D., Hanson, S. J., & Glymour, C. (2011). Multi-subject search correctly identifies causal connections and most causal directions in the DCM models of the Smith et al. simulation study. NeuroImage, 58(3), 838-848. This paper also treats the distributions as a mixture, though doesn't assume non-i.i.d.

We respectfully disagree. We note that several contemporary papers on time series causal discovery use this dataset to benchmark their methods [1, 2, 3]. Furthermore, we examined the paper cited by the reviewer and discovered that the paper demonstrated that the IMaGES algorithm performed better than other contemporary i.i.d causal discovery algorithms; however, they did not make a claim that analyzing it as i.i.d yields better results than analyzing it as time series. Indeed, they did not compare against any temporal causal discovery baselines.

As a result, it seems that any study proposing a time series analysis of this data should do a comparison of this result to one obtained by treating the data as i.i.d. instead, since this is a known phenomenon for this particular dataset. This is an issue because the proposed method is specifically designed to deal with time series.

Following your suggestion, we benchmarked well-known i.i.d causal discovery algorithms against simulation 3 data from the NetSim dataset. We were unable to run the IMaGES algorithm from [4] since we could not find a publicly available implementation. We compare the performance against reported numbers from temporal causal discovery methods in literature. The results are summarized below:

We observe that the temporal causal discovery algorithms handily outperform the i.i.d causal discovery baselines.

The problem with the DREAM3 examples is that none of the reported results have any lift with AUC; they’re all stuck around 0.5, which is more or less random. The differences between the methods are slight. This result doesn’t give a case where the proposed method is particularly helpful.

DREAM3 is a challenging dataset due to the low number of samples (46 samples per graph) in the dataset. Nevertheless, we quite significantly outperform the baselines. More importantly, the clustering accuracy is high (94.83%) as reported in Section 5.3, which means that despite being unable to infer the causal relationships themselves, our method can distinguish among data coming from different causal models despite the limited samples. We believe that this is most definitely a significant outcome.

Questions:

Can you find better examples to show the usefulness of the theory?

Identifiability is a cornerstone theory in causal discovery. Without it, any empirical case studies would be meaningless. That is, if there is no guarantee that there exists an underlying true set of graphs that our algorithm can find, then even strong empirical performance would be in vain. We would like to point out that the theory results presented in the paper pertain to the question of structural identifiability of the mixture model, which is a general statistical question rather than pertaining to our model specifically. Identifiability results specify conditions under which the constituent causal graphs can be uniquely determined from the mixture distribution.

Where can the theory really shine so far as the empirical application is concerned? Where can it be expected to do significantly better than alternative methods?

Although not particularly pertaining to the theoretical results presented in the paper, we would expect our method to work better than alternative methods when the underlying causal graphs are sufficiently diverse and different from each other.

References: [1] Gong, Wenbo, et al. "Rhino: Deep causal temporal relationship learning with history-dependent noise." arXiv preprint arXiv:2210.14706 (2022). [2] Cheng, Yuxiao, et al. "CUTS: Neural Causal Discovery from Irregular Time-Series Data." arXiv preprint arXiv:2302.07458 (2023). [3] Löwe, Sindy, et al. "Amortized causal discovery: Learning to infer causal graphs from time-series data." Conference on Causal Learning and Reasoning. PMLR, 2022. [4] Ramsey, J. D., Hanson, S. J., & Glymour, C. (2011). Multi-subject search correctly identifies causal connections and most causal directions in the DCM models of the Smith et al. simulation study. NeuroImage, 58(3), 838-848.

We appreciate the reviewer for their valuable feedback and inquiries. In the subsequent section, we will address the key concerns that have been raised.

Weaknesses:

The proposed objective function can be seen as a straightforward extension of the standard variational inference optimisation framework. Therefore, the overall novelty and significance of the work may be somewhat limited in this regard.

We politely disagree. The problem of inferring multiple causal graphs is quite challenging and has received scant attention in the literature despite its widespread applicability. Our method makes significant contributions with (1) designing an effective deep latent variable model and deriving an ELBO objective to infer the relevant parameters (2) proving the identifiability conditions and (3) solid empirical results. Even though the objective function is an instance of variational inference, it is not the novelty that we claim in this paper. Despite its relative “simplicity”, it achieves competitive empirical results.

A major drawback of the work is that the reported results are based on training data, as the proposed method depends on learnable sample-specific parameters. It raises questions about why an encoder producing a K-way categorical random variable given a sample was not considered by the authors.

There seems to be a misunderstanding of the problem. Causal discovery is an “unsupervised” problem, in the sense that the ground truth causal graph(s) are assumed to be unknown. Reporting causal discovery results on training data is a standard practice in the literature. We point out that many similar works on optimization-based causal discovery methods work on similar principles, e.g. [1,2,3] Our proposed method reframes the task to an optimization problem, which we solve by training neural networks. In other words, our inference procedure itself involves training on the input data $X$. We have clarified this in the Problem Setting in Section 3 of the paper.

We learn sample-specific parameters, as the task involves inferring membership information solely for the samples within the input data $X$. There is no need for generalization to unseen points. Utilizing an "encoder producing a K-way categorical random variable" would essentially achieve the same objective of learning sample-specific parameters. Using a K-way encoder is a good suggestion that would work especially well for a large sample size. However, since all our experiments are run in the low-sample regime, sample-specific weights work sufficiently well with low memory footprint.

The lack of reported results on generalisation limits insight into the proposed method's ability to perform beyond the training data.

As noted earlier, the causal discovery problem only involves inferring causal graphs for the training data $X$, and hence the notion of “generalization beyond the training data” is not applicable in this context.

Although the method claims flexibility in terms of core causal structure learning algorithms, performance is only demonstrated based on one causal discovery method.

We focus on one “core” causal structure learning algorithm (i.e. Rhino [1]) in order to thoroughly evaluate its capabilities and performance. While other algorithms could potentially be more suitable for specific applications, we find that MCD with Rhino seems to work reasonably well for the tasks we selected and to demonstrate the utility of our framework.

Questions:

Can you please clarify how the AUROC metric is computed for the different methods being compared?

For Rhino and MCD, which parameterize the causal graphs as Bernoulli distributions over each edge, we use the inferred edge probability matrix as the “score”, and evaluate the AUROC metric between the score matrix and the true adjacency matrix. For DYNOTEARS, we use the absolute value of the output scores and evaluate the AUROC. For PCMCI+ and VARLiNGAM, since they only output adjacency matrices, we directly evaluate the AUROC between the predicted and true adjacency matrices. We have updated Appendix D.1 with this clarification.

Regarding the experiment on Netsim-permuted data, wouldn’t one expect the results of the proposed method to be inherently favorable compared to non-mixture alternatives? This is because by permuting the variables, we intentionally manipulate the data to explicitly align with the underlying assumption of a mixture model.

We performed random permutation of the nodes to demonstrate that even though methods that predict a single graph can perform reasonably well when the different causal models are similar, they are susceptible to poor performance under very mild “perturbations” i.e. permutation of the nodes. Through node permutation, we generate a mixture distribution where the constituent causal graphs differ significantly. This process validates our hypothesis that the proposed method (MCD) should be effective in such scenarios.

We would like to thank the reviewer for their valuable feedback and thoughtful inquiries. Below, we respond to the main concerns that have been raised.

Weaknesses:

my main skepticism is due to the fact that all the results are for the training set, which I am surprised (and had a stronger positive impression until that point). Trivial enough: For an unquestionable positive score, I would rather need result on test data, on unseen graphs.

As pointed out in the reply to Reviewer vYSG, we reiterate that this is standard practice in causal discovery, where the causal graph is inferred without any supervision. We do not assume access to a training corpus of ground-truth graphs or functional relationships; as such, all results are on “unseen graphs”. We have clarified this in the Problem Setting in Section 3 of the paper.

How different the causal graphs in these mixture models (and also the data i.e., average SHD within the cluster is missing. This is really important to really understand what is going on behaviour of the method. (hence, my question). No surprise to see its effect on highly-imbalanced one.

Within the cluster, the average SHD would be 0, since we predict the same graph for each cluster. We plot the average pair-wise SHD (over all pairs of causal graphs) for each dataset here:

on minor point, an illustrative toy example is lacking would be very useful.

Thank you for the suggestion. We have added an illustrative toy example with two causal graphs in Appendix E to explain the importance of using a mixture model for heterogeneous causal discovery.

Questions

What is the effect of the distance between causal graphs to the performance?

We believe that our model is better able to cluster points when the SHD between causal graphs is larger. This is especially prominently seen with the Netsim and Netsim-permuted datasets.

g_1 and g_2 are not defined in Theorem 1, is it a typo: should be h? or a particular graph?

Yes, it is a typo. They should be $h_1$ and $h_2$. This has been fixed in the new draft.