# Answer to Reviewer Vx9e

We answer here the reviewer's follow-up questions.

**Regarding the I-NLL metric**. The I-NLL is the importance-weighted ELBO (Burda et al., 2016) of each model, evaluated for a set of hold out data points with different interventions. The formula for the iwELBO appears in Appendix C. Because all models have similar building blocks (e.g., same neural network architectures, and same observation likelihood), we benchmark the impact of the modelling choice aschoice of how the prior $p(z | a)$ changes with respect to the action $a$. Also, we note that, as expected, those hold-out interventions are the same for all baselines. With that in mind, the I-NLL measures how well each model fits the data after the transfer learning procedure. We will add a few sentences to clarify this in the main text.

**Regarding the variability assumption**. We thank the reviewer for clarifying the question. In the original sVAE manuscript, there are two kind of assumptions related to the sparse mechanism shift hypothesis. The first one relates to time-measurements, and sparsity in interactions between latent units through time (sparse $G_z$). The second one relates to sparsity between interventions and latent variables (sparse $G_a$). The main theorem of Lachapelle et al. (Theorem 5) encapsulates both assumptions, and therefore also has two types of sufficient variability, one for $G_z$, and one for $G_a$. In our work, we *do not* consider time varying measurements, so only consider sufficient variability with respect to Ga, which cannot be relaxed, neither for sVAE nor for sVAE+ (as we claimed in our previous response). However, the sufficient variability for $G_z$ can be relaxed, and this special theoretical analysis is already conducted in Lachapelle et al. (Theorem 22). We will make this clear in the manuscript to avoid confusion.

# Answer to Reviewer Vx9e

We thank the reviewer for their feedback and are glad they appreciated the potential impact for the biological community.

**Regarding the results on simulated data**. We agree with the reviewer that the simulated data resembles the generative model from sVAE+. However, we would like to point out that the generative model of sVAE shares a similar level of similarity to the simulated data, because both models assume sparsity of the Bernoulli probabilities $\pi$, as explained in paragraph “Connection to prior work” of Section 3.2. Because, sVAE+ still provides a substantial improvement over sVAE across several of our benchmark datasets, we believe this is strong evidence that the method improves over sVAE.

**About the high recall for iVAE**. We thank the reviewer for pointing this out. Unfortunately, the iVAE method does not directly identify a causal relationship, because it simply assumes a Gaussian location family in the prior of $z$ for the interventional data. Consequently, we need a heuristic for constructing a graph, using the inferred parameters $\mu_a$. In this work, we used an outlier detection method, but it seems to call too many edges (as noted by the high recall score compared to other methods). We have investigated a variant that instead selects the top 2 affected components of $\mu_a$ (in terms of absolute value), and here the Precision, Recall and F1 scores are: 0.16, 0.18, 0.17 for the experiments in Table 1. Those results are still lower than sVAE+.

**About the significance of results in Section 5**. We thank the reviewer for pointing this out. We are currently rerunning the results for Section 5, and will update them for the camera ready version. We may not be able to report the entire updated results before this rebuttal period, because each run requires around 20 min of computing time, and Table 3 only comprises 5 seeds x 8 data splits x (5 methods + hyperparameters) runs, totalling around a week of wall clock time on a single GPU.

However, we managed to obtain a subset of the results in a short amount of time. More precisely, we report the mean and standard deviation of the results for the first (Exosome and mRNA turnover) and second column (Spliceosome) of Table 3, for the top-2 performing models (iVAE and sVAE+), using 5 random seeds. The standard deviations are very small, and the difference in mean are significant under an unpaired t-test (p < 0.0001).

**Exosome and mRNA turnover**
| I-NLL | iVAE | sVAE+ |
| -------- | -------- | -------- |
| Mean | 1784.118 | 1782.831 |
| SD | 0.140|0.129|
|SEM| 0.062 | 0.058|

**Spliceosome**
| I-NLL |iVAE| sVAE+
| ------- | ----- | ----- |
| Mean |1779.235 | 1777.817
| SD | 0.142 | 0.159
| SEM |0.063 |0.071

We have proceeded similarly for Table 2, but comparing the results of the top-2 performing methods (sVAE and sVAE+). Here as well as have significant differences (p < 0.0001)

**Norman et al.**
| I-NLL | sVAE | sVAE+|
| ------- | ----- | ----- |
|Mean |776.271 | 775.441|
| SD | 0.084 | 0.102|
| SEM | 0.037 |0.046

In the light of these results, we believe that sVAE+ **significantly** outperforms other methods on real data sets.

**About the conditions for identifiability between sVAE and sVAE+**. As far as we understand, the generative models of all of sVAE+, sVAE, and iVAE, as inspired by the ICA framework, assume conditional independence of the latent components $z$ given the auxiliary covariates $a$. The main difference between both of sVAE+ and sVAE compared to iVAE is the sparsity assumption, which is one of the sufficient conditions for identifying the latent variables up to a permutation.
However, to the best of our knowledge, treating the proportions from the Bernoulli distribution from a Bernoulli-Beta instead of a Bernoulli with a Sigmoid-Laplace prior, is such a minimal change to the model that it may not warrant relaxation of the sufficient variability assumption.

# Answer to Reviewer Lstr

We appreciated the reviewer's feedback, and are glad that they appreciated the strength of the experiments.

**About the technical contributions in this manuscript**. We agree with the reviewer that this paper may not be adequate if we were claiming to introduce new methodology for causal representation learning. Instead, we introduced a new application to the theory of Lachapelle et al., with several technical components:
1. the reformulations of existing VAEs for single-cell RNA sequencing data
2. the development of sVAE+, a new VAE variant that outperforms sVAE
3. the extension of the benchmarking done in Lachapelle et al. with new metrics (OOD evaluation), new robustness tests, and open-source sandbox
4. a practical discussion of the theoretical assumptions in Lachapelle et al. paper, and whether they hold in a practical context.

Our modification of Lachapelle et al.’s approach to adjust for single-cell data, can also serve as a guiding procedure for practitioners that wish to apply sparse VAEs on other data modalities (such as medical images, sequences etc.), namely, how to change the generative model and the inference.

Considering those contributions in the context of an application paper, we hope that the reviewers will see this manuscript fit in the categories (i) Causal generative models for machine learning (ii) Benchmark for causal discovery and causal reasoning and (iii) Applications of any of the above to real-world problems.

# Answer to Reviewer orYU

We thank the reviewer for their review, and are glad they enjoyed the idea and the quality of the experiments.

**About the noise model for scRNA-seq**. Thanks to the reviewer for pointing this out. Briefly, gene expression data captured by single-cell RNA sequencing can be composed as technical noise, as well as biological signal. The technical noise is complex, but usually described by two components: (A) the number of RNA transcripts captured in a single-cell is (mostly) treated as an artifact of the assay, and must be factored out of the learned representation. (B) the data takes the form of counts, with overdispersion. The best noise distribution is (still) a topic of active research, but there is a general consensus that the noise model presented by Lopez et al. (2018), with a scaled negative binomial distribution provides satisfactory results (Grün et al. 2014). We will add a longer discussion of the appendix in order for the paper to be self-contained.

Grün et al. Validation of noise models for single-cell transcriptomics. Nature Methods (2014)

**About the presentation**. We agree with the reviewer's suggestions for improving the presentation, and will add such elements to the final version of the paper.

# Answer to Reviewer t3PU

We thank the reviewer for their feedback, and for appreciating the quality of the experiments, as well as the importance of the problem we considered.

**About the discussion in Appendix G**. We thank the reviewer for bringing this important point. We will add a brief summary of some of the points in Appendix G in the main paper. We also conducted a follow-up study in which we target a subset of the latent variables, and report the evolution of the Pearson MCC for sVAE+ with varying subset size. As expected, performance gets worse when a lower number of latent variables are targeted. We will add this results to the Appendix.

| || | | | |
| -------- | -------- | -------- | --| --| --|
| Number of targeted latent variables (out of 15) | 1 | 4 | 7 | 12 | 15|
|MCC | 0.5 | 0.59 | 0.71 | 0.85 | 0.88|

**About the identifiability of the prior in the transfer learning experiments**. In our transfer learning experiments, we fix the parameters of the conditional distribution $p(x | z)$ and learn the parameters for the prior $p(z | a)$ (that is, $\pi$ and $\gamma$), for new interventions a, from data $x_a$. In this setting, recovery of the latent variables $p(z | x_a, a)$ is not guaranteed. However, our intuition is as follows (informal). If (A) the interventions in the training data are enough to correctly estimate $p(x | z)$ during training (i.e., disentanglement) and (B) the variational posterior $q(z | x)$ closes the variational gap, then the problem reduces to maximum likelihood for the parameters of the prior $p(z | a)$. Because this prior is essentially a Gaussian location family, its parameters are identifiable (because we fixed $p(x | z)$). Consequently, the transfer learning experiment is set up in a way that we expect the interventions on the *training* set to identify the latent variables, and we evaluate how well those latent variables can be used to describe the data points from held-out interventions.

The (informal) sufficient conditions described above are not satisfied in practice for our experiments in Table 3 of Section 5, because we are holding out interventions targeting entire pathways at a time. In this case, this presumably guarantees only partial identification of the latent units. However, we believe that the transfer learning procedure is still a reasonable way of quantifying the adequacy of the learned latent variables for describing hold-out interventions.
From the results in Appendix D, we also see how the performance decreases when the hold out interventions have undergone some shifts and they do not correspond to the distribution of the interventions seen in training.

We will incorporate those nuances in the appropriate sections of the papers, and add more details about the transfer learning experiments in Appendix.