Distinguishing Cause from Effect with Causal Velocity Models

We thank the reviewer for their critical reading of our paper. The reviewer raises a few concerns and points of confusion that we believe can be addressed.

Reliance on score/density estimate: Please see response to reviewer 9rbC.

Questions for Authors:

1. Do flows represent ANMs? Yes. We can regard the velocity/flow parameterization as a new way to express SCMs. ANMs, PNLs, and LSNMs are special cases (see Table 1). The results in Tables 2-3 refer to the velocity parametrization/estimation of ANM/LSNM. We can also express new SCMs easily via their velocities (see Table 1, Fig. 2, and Sec. D.2 for examples). In fact, Thm. 3.2 tells us that bijective SCMs and velocity models are in one-to-one correspondence. One of the advantages of our method is that entirely new classes of SCMs based on velocities can be fit and used for causal inference/discovery. ANMs are known or suspected to be mis-specified in many settings, in which case they may not be reliable for causal inference/discovery.
2. Performance of method with varying sample sizes: The method is sensitive to sample size mainly through the accuracy of score estimation. Given the exact score, Eq. 14 holds pointwise, so there is no statistical element remaining in minimizing Eq. 18. We designed an experiment to evaluate this empirically. Please see our response to reviewer 9rbC, Tables 9rbC.1 and 9rbC.2 for details.
3. Comparisons to additional methods: We have evaluated ANM + independence tests and other methods from the cdt package, as well as CGCI, on the Tuebingen dataset (after filtering out discrete/ordinal datatsets, for which our model is mis-specified; see our response to reviewer sNSP) as well as the three simulated datasets in the paper (Table H4Vs.1 below). We find that our method performs best overall, although IGCI works well due to the Gaussian noise when generating the data.
4. Difficulty of optimizing loss: Please see our response to reviewer xsfW about computational considerations and optimization.

Table H4Vs.1: Performance of our method compared to other methods

Questions on Experimental Design/Analysis

• LOCI not only SOTA method. See answer 3 above.
• Noise level: To our knowledge, methods such as NOTEARS have been shown to be sensitive to marginal variances, in particular where the variance in the effect is larger than that of the cause. Following recent convention in the field first established by [Reisach et al., 2021], we standardize the data prior to causal discovery. Note that we use a noise scale of 3 when generating data to reduce the signal-to-noise ratio and increase the difficulty of causal discovery, but this affects all methods equally after standardization.
• Sample size. See answer 2 above.
• Optimization of loss. See answer 4 above.

Questions on Claims and Evidence: Optimized loss in causal direction as sample size diverges. Please see Table 9rbC.1, which shows that when the exact score is known and used, the optimized loss is small even for relatively small $n$, and appears to tend towards 0 for the causal direction, while it stays an order of magnitude larger in the anti-causal direction. For estimated scores, Table 9rbC.2 shows a similar trend towards increasing discovery accuracy as score estimation improves with increasing $n$, which is supported theoretically by Thm 6.1.

On Other Comments/Suggestions/References: Thank you for the detailed suggestions. Due to space restrictions we cannot respond to each in detail. We will add the suggested clarifications and references in revisions. We would like to clarify that in the context of Sec. 2.3, $s$ can be viewed as a “starting time” and $t$ as a time variable at which we evaluate the flow that was started at time $s$. Sec. 2.3 is background on flows and ODEs. Starting with Sec. 3, we substitute a cause variable $x$ to play the mathematical role of time, as stated at the start of Sec. 3.

We thank the reviewer for their careful reading and analysis of our paper. The reviewer lists a few weaknesses and questions that we would like to address.

Extension to multivariate settings: Please see our response to reviewer xsfW.

Performance on cause-effect pair dataset collection: Indeed, there is room for improvement. There are a number of possible explanations.

• If datasets are generated by a process for which a specific model class (e.g., ANM or LSMN) is well specified, we would expect the specific model-based methods to perform better. Whether that is the case for datasets in the collection is not something we can test.
• Some of the constituent datasets in the Tuebingen collection have discrete or ordinal data, while our framework requires densities that are continuous on $\mathbb{R}$. We re-ran the experiments after discarding 29 total datasets that have discrete or ordinal data (e.g., “rings” in pair 5 is integer-valued). On this data, our best performing methods (B-QUAD and B-LIN with KDE), obtained 83% (B-QUAD) and 79% (B-LIN) accuracy (up from 59% and 69%, respectively) which is more competitive with other methods. We also re-ran experiments with competing methods in this setting and found that their performance generally did not improve after removal of discrete/ordinal datasets. Please see our response to reviewer H4Vs, specifically Table H4Vs.1 for full results.

Violation of assumptions:

• Non-invertible mechanisms are not an issue for our method. They would be an issue for making unit-level counterfactual inference, but our method only relies on the distributions and their representations via functional models (i.e., SCMs), so there is no loss of generality in assuming that the conditional distribution can be represented by a bijective SCM. This is noted after Eq. (1) in the paper, but we will be sure to emphasize the point more clearly in revisions.
• Confounders: We have not studied this setting and therefore cannot say anything specific, though we consider it important future work. We note that all methods that assume causal sufficiency would likely have issues in the presence of confounders.

Computational complexity: We thank the reviewer for pointing this out. We will add a discussion to the appendix. Please see our response to reviewer xsfW about computational considerations.

Performance of B-QUAD: We believe that the B-QUAD model is generally flexible enough to fit most data, while still being simple enough (it only has 9 scalar parameters, see Eqn 67 + 68) to distinguish the causal direction. See Figure 6 of the Appendix for an example of fitting B-QUAD causal curves to real data. Developing a better understanding of this trade-off is an interesting direction for future work.

Using the velocity parameterization for SCMs: A great point–we will add a discussion in revision. In short: In the bivariate scalar setting, the velocity parameterization could be used to define novel classes of functional causal models for causal modelling. Since they generate bijective SCMs, they would inherit the counterfactual identifiability results of BCMs [Nasr-Esfahany et al., 2023]. In particular, they can be used to compute queries at all levels of Pearl’s ladder of causation. We can use numerical integration to evaluate the causal curve (Figure 1). Note also that the abduction and prediction steps are especially simple in BCMs, only involving inverting and forward evaluation of the model. In particular for velocity models, this just corresponds to integrating the velocity from an observed condition $x$ to a query condition $x'$, see the discussion above Eqn 2.

We thank the reviewer for their careful reading and analysis of our paper. The reviewer seems to have appreciated the main benefits of the proposed method (bivariate causal discovery with minimal model assumptions), with some reservations about its reliance on nonparametric score estimation. Our rationale for the reliance on a nonparametric estimate is as follows:

• We view this as an instance of the fundamental statistical trade-off between strength of assumptions and required data/computation. In causality, ANMs sit near the strong assumptions/small sample end of this spectrum. LSNMs already require more data/computation (i.e., two-stage regression) in general. Our method sits further into the weak assumption/more data regime. So do others, such as doubly-robust estimation methods. We would argue that all of these methods are interesting; each is the right tool for some class of problems, with weak assumptions being preferred in settings with limited a priori evidence for favoring specific models.
• Our approach is to reduce causal discovery to score estimation. We do not claim to introduce a tool that works well on every data set. Our goal is to design a tool that turns reasonable score estimates, where available, into conclusions about cause and effect. This is conceptually similar to methods that reduce causal discovery to independence tests.
• We find this conceptually appealing since it cleanly separates the causal aspect of the problem (computed once the score estimate is obtained) from the statistical ones. All statistical aspects are, in a sense, encapsulated in the score estimate—the dependence on sample size, for example, only enters through the score estimate.
• Since score estimation is itself an active research area, we also point out that our method is agnostic to how the score estimate is obtained; it only depends on the quality of the estimate. As new or improved score estimators for a given problem become available, these can simply be plugged into our method.
• We also note that recent work on causal discovery in ANMs (e.g., NoGAM [Montagna et al., 2023]) also relies on score estimation as a first stage of the discovery method.
• For computational complexity, please see our response to reviewer xfsW.

Discovery with known score

To highlight the separation of causal/statistical components, we performed the following experiment. Using synthetic data generated from an ANM (Gaussian noise; mean function is a 3-layer MLP with random weights and tanh activation) such that the score could be computed analytically/numerically in both directions, we compute the ground truth scores as input to our method and find that sample size plays no role (beyond $n>10$, which is required to evaluate the criterion at a sufficient number of approximation points), and causal direction can be inferred with certainty. We also report the GoF statistic from Eq. 18, showing that it is an order of magnitude lower on average in the causal direction.

Table 9rbC.1: Velocity-based causal discovery with known score (results averaged over 100 replications)

Discovery with estimated score

We repeated the experiment using the same synthetic data as above, using varying sample sizes to estimate the score and subsequently perform velocity-based causal discovery (as in the experiments reported in the paper). We find that as score estimation improves, the success rate of causal discovery also improves.

Table 9rbC.2: Velocity-based causal discovery with score estimated by Stein score/Gaussian kernel (results averaged over 100 replications)

We thank the reviewer for their comments.

Extension to multivariate settings: The multivariate case is not a straightforward extension of the bivariate setting. To extend the velocity interpretation we would need to carefully define a multivariate time. We agree that this extension is interesting, but we have not identified a best approach yet and consider this future work.

Scalability: The entire procedure, even for $n=5000$, takes less than 3 seconds for KDE estimators, and less than 20 seconds for the Stein estimator on a 2020 M1 Macbook air. The bottleneck is the score estimation, which is asymptotically $O(n^3)$. However, KDE only requires a single matrix multiplication, while Stein only requires several matrix multiplications and inversions, ahead of optimization. One can also use nearest-neighbour based subsampling to reduce the complexity to $O(n \log n)$, which enables scalability to datasets with $n \gg 10000$. Once the score is estimated, optimization of Eq. 18 is a cheap $O(n)$ regression task given the score estimate.

To give this discussion some detail, we report the average computation times/steps for the experiment reported in Table 9rbC.2 of our response to reviewer 9rbC.

Table xfsW.1: Computational performance of the proposed method (averaged over 100 replications; velocity optimization convergence tolerance is 1e-6)