Thank you for your extensive response; we truly appreciate your time and constructive feedback.

Please consider our responses to your comments below.

Response to General Comments

Presentation

• In the related work, we will include a reference to time series methods such as RPCMCI. These methods use EM to discover regimes along the temporal dimension and, therefore, additionally rely on the temporal ordering of the data.
• We also value the suggestions for streamlining our formulation of the causal model; we agree with your viewpoint and will express $\beta$ to be a function of the respective mixing variable. We also note that you are correct that $Z$ are themselves causal parents, but we write $\text{Pa}_i$ to distinguish the observed from the unobserved parents.
• We appreciate the effort taken to point out our typos, we have corrected them in the current manuscript.

Nonlinear Models

While we consider focusing our analysis on linear models useful on its own, we would like to justify our claim on the "easy" extensibility of our framework to nonlinear models. We have therefore we conducted a small proof-of-concept study to support this claim.
Please consider our results, below.

In the following, we replace the linear-regression mixtures with nonlinear variants, here a natural spline. We again use the EM algorithm to infer the mixture components and the BIC criterion to pick their number. The following table shows how well we can reconstruct the mixture components (Fig. 4 in the manuscript), measured by Adjusted Mutual Information (AMI) averaged over each node in 10-node graphs.
We show the CMM with the MLR in as our main presentation (left) and the one with a mixture of natural splines (right) under different true generating functions (rows).

The results match our expectations, where the models perform very closely in the linear case, but the MLR degrades in performance under non-linearity. Replacing it with the mixture of splines allows covering the non-linear cases reasonably well.

This proof-of-concept experiment supports exploring our algorithm with more flexible regression mixtures, and we can include it as an ablation study to encourage this. Indeed, our guarantees for the latent-aware oracle BIC can be transferred to a nonlinear setting. However, this would require access to oracle mixture parameters, and the connection to estimates obtained using EM is less clear under non-linearity. This point would require a deeper analysis in future work, but we hope that the reviewers agree that the parametric case is of sufficient contribution for the scope of this work.

Implementation Details

As we cannot provide source code here, we summarize the changes we made to the implementation. The source code included in the supplement requires the following small changes,

1. using a different regression-mixture within our algorithm. This can be done in src/mixtures/mixing/mixing.py by changing the syntax of the regression formula passed to the R-script.
2. using a different data-generating function $f$ for each causal relationship in the randomly generated DAGs. Implementation support for this already exists under src/exp/gen/synthetic/data_gen_mixing.py where the linear function used in each cause-effect relationship can be replaced by other functional forms. To visualize these functional forms, one can refer to the notebook src/examples/demo_mixing.ipynb and modify the parameter F, such as params['F'] = FunType.CUB for cubic functions.

The experimental setting corresponds to the same base parameters that we used in main manuscript.

Response to Questions

Motivating Example

This is a good point. We can indeed extend our motivational example to motivate multiple latent mixing variables. For example, independent latent variables could correspond to (1) unknown batch effects due to different laboratories from which measurements come and (2) prior exposure related to demographic location. Following your suggestion, we improved the motivational premise in the updated manuscript as follows.

Take for example a country-wide study of anti-microbial resistance in hospitalised patients focused on a resistant pathogen such as Methycillin-Resistant Staphylococcus Aureus (MRSA) [4]. For each patient, phenotypic susceptibility to multiple antibiotics is measured; for each antibiotic, each medical center may choose one of the collaborating laboratories to perform each test, thereby introducing a latent variable corresponding to a batch effect, each independent from the others.
Further, each patient's medical history includes prior exposure to pathogens with known cross-resistance to MRSA, such as Enteroccocus [5]. The plasmid profile of the latter in each geographical region largely affects its susceptibility, say to Vancomycin [2], which in turn affects the cross-resistance of MRSA [1]. Despite being well documented, this variable is not routinely measured by medical centers, and is hence an independent latent variable that defines the mechanism under which the presence of Enteroccocus affects the MRSA cross-resistance [3].
Such settings, where multiple observed variables are each affected by one independent latent factor, are common in practice. Yet, this premise is not currently addressed by standard causal inference methods, which often rely on assumptions like causal sufficiency or no unmeasured confounders. Both of these assumptions are violated in this setting, often leading to the discovery of spurious dependencies, or to missing true causal relations, altogether [6, 7].

* [1] Arredondo-Alonso S, Top J, McNally A, Puranen S, Pesonen M, Pensar J, Marttinen P, Braat JC, Rogers MRC, van Schaik W, Kaski S, Willems RJL, Corander J, Schürch AC. Plasmids Shaped the Recent Emergence of the Major Nosocomial Pathogen Enterococcus faecium. mBio. 2020 Feb 11;11(1):e03284-19. doi: 10.1128/mBio.03284-19. PMID: 32047136; PMCID: PMC7018651.

* [2] Boumasmoud M, Dengler Haunreiter V, Schweizer TA, Meyer L, Chakrakodi B, Schreiber PW, Seidl K, Kühnert D, Kouyos RD, Zinkernagel AS. Genomic Surveillance of Vancomycin-Resistant Enterococcus faecium Reveals Spread of a Linear Plasmid Conferring a Nutrient Utilization Advantage. mBio. 2022 Apr 26;13(2):e0377121. doi: 10.1128/mbio.03771-21. Epub 2022 Mar 28. PMID: 35343787; PMCID: PMC9040824.

* [3] Cong Y, Yang S, Rao X. Vancomycin resistant Staphylococcus aureus infections: A review of case updating and clinical features. J Adv Res. 2019 Oct 12;21:169-176. doi: 10.1016/j.jare.2019.10.005. PMID: 32071785; PMCID: PMC7015472.

* [4] Hasanpour AH, Sepidarkish M, Mollalo A, Ardekani A, Almukhtar M, Mechaal A, Hosseini SR, Bayani M, Javanian M, Rostami A. The global prevalence of methicillin-resistant Staphylococcus aureus colonization in residents of elderly care centers: a systematic review and meta-analysis. Antimicrob Resist Infect Control. 2023 Jan 29;12(1):4. doi: 10.1186/s13756-023-01210-6. PMID: 36709300; PMCID: PMC9884412.

* [5] Li G, Walker MJ, De Oliveira DMP. Vancomycin Resistance in Enterococcus and Staphylococcus aureus. Microorganisms. 2022 Dec 21;11(1):24. doi: 10.3390/microorganisms11010024. PMID: 36677316; PMCID: PMC9866002.

* [6] Huang B, Zhang K, Zhang J, Ramsey J, Sanchez‑Romero R, Glymour C, Scholkopf B. Causal Discovery from Heterogeneous/Nonstationary Data: Skeleton Estimation and Orientation Determination. Journal of Machine Learning Research. 2020 May;21(89):1–53.

* [7]  Mazaheri B, Squires C, Uhler C. Synthetic Potential Outcomes and Causal Mixture Identifiability. In: Proceedings of the 28th International Conference on Artificial Intelligence and Statistics (AISTATS 2025); 2025 May; Mai Khao, Thailand. PMLR, vol. 258.

Algorithm Notation

In the description of the algorithm, we will refer to $\text{Pa(Y)}$ instead of $\text{Pa}_j$, thank you for pointing out this inconsistency. To clarify, we refer indeed to the edge $X \rightarrow Y$ that we consider including in the current parent set $\text{pa}(Y)$ for $Y$, for which we fit mixture models under these given parent sets.

We are happy to address any additional points that may arise.

Dear reviewer, thank you very much for your response.

We propose the following mild changes to Section 2, and welcome any suggestions.

To formalise this, we also consider a set of discrete, unobserved random variables $\boldsymbol{Z} =${$Z_1,\ldots,Z_m$}, with $m \leq n$, each following a categorical distribution $Z_i \sim \mathrm{Categorical}(\boldsymbol{\gamma}^i)$ with $Z_i \in \{1, \ldots, K_i\}$. That is, each $\boldsymbol{\gamma}^i$ lies on a $K_i$-dimensional probability simplex $\boldsymbol{\gamma}^i = (\gamma^i_1, \ldots, \gamma^i_{K_i})$ with $\sum_{k=1}^{K_i} \gamma^i_k = 1$, so that $\mathbb{P}(Z_i = k) = \gamma^i_k$. We call these random variables mixing variables.

The causal mechanism of each observed random variable in $\boldsymbol{X}$ depends on a set of observed causal parents as well as exactly one of the unobserved, latent $\boldsymbol{Z}$, as determined by the surjective map $\mathrm{La} : \boldsymbol{X} \to \boldsymbol{Z} : X_j \mapsto Z_i$, in which case we simply write $\mathrm{La}_j = Z_i$. The mixing variable directly affects the parameters of the causal mechanism, which we therefore express as a function $b_j: [K_i] \to \mathbb R^{|\mathbf{Pa}_j|}: k \mapsto \beta^{jk}$ mapping each value $k$ of $Z_i$ to a linear coefficient vector $\beta^{jk} \in \mathbb R ^{|\mathbf{Pa}_j|}$. Hence, the parameters of the functional dependency $f$ are the collection of vectors $B_j=(\beta^{j1}, ..., \beta^{j K_i}$); this consists of one coefficient vector $\beta^{jk}$ for each mixing coefficient $1\leq k \leq K_i$, and each such vector has a dimension equal to the number of parents $\mathbf{Pa}_j$ of $X_j$.

Summarizing, we can now model each random variable $X_j$ as generated from its observed causes $\mathbf{Pa}_j \subseteq \boldsymbol X$ by the causal function $f$ and the coefficients $b_j$ that depend on the latent $Z_i = \mathrm{La}_j \in \boldsymbol{Z}$, where we recall that $\boldsymbol Z \cap \boldsymbol X = \emptyset.$ Then, we have

$X_j = f(\mathbf{Pa}_j, b_j)+ N_j \qquad \text{ with }\qquad f( x , b_j (k)) =$ $\beta^{jk \top}$ $x \qquad \qquad (1),$

where $N_j \perp\kern-5pt \perp \mathbf{Pa}_j$ is additive Gaussian noise, $N_j \sim \mathcal{N}(0, \sigma^2)$.

The main changes follow your suggestions, where we

• explicitly write $\mathbf{Pa}_j \subseteq \boldsymbol X$ to refer to the observed causal parents while $\mathrm{La}_j$ is a single unobserved parent, and
• express the linear coefficients as a function $b_j$ of $\mathrm{La}_j$, so that the dependence on the latents enters directly (and only) through $b_j$ instead of appearing as a general argument of $f$; see Eq (1).

---

Note: There seems to be an issue with LaTeX code display where subscripts of two or more letters, as in \beta_{jk}, result in display errors when used in longer equations. Since the same works fine for superscripts, \beta^{jk}: $\beta^{jk}$, we temporarily moved the indices for all $\beta$ to the superscript in our response above. In the manuscript, we will keep using subscripts as in the submitted version.

We thank you for your time and thought-provoking comments.

Please consider our response to your questions below.

Response to Questions

1. Experimental States in Sachs et al.

Yes, there is a known correspondence of the mixing variables to experimental conditions for this data.

The dataset arises from single-cell protein-signalling data of the human immune response system, wherein each of the $11$ observed variables corresponds to the activity of one compound of interest: either a protein or a phospholipid. The dataset contains (among others) 5 experimental conditions, in each of which a particular molecular modifier has been applied to the cells, such that the activity of exactly one of the $11$ compounds of interest is affected. This results in a known change of the measured activity for the corresponding compound.

In our experiments, we combine the data from all experimental conditions, and we do not make their origin available to the algorithm. Hence, in the pooled data, each variable is affected by a single latent one. The usefulness and relative prevalence of this dataset as a real-world case study among causal inference methods arises from both

1. the well-studied underlying causal mechanisms, and, in our case, also
2. from the known ground-truth values of the resulting latent variables.

For example, considering the node "Akt", there are two mixture components: one is the experimental condition where the so-called Akt inhibitor was applied, inhibiting the activation of Akt directly; the other mixture component comprises the remaining samples where Akt is in its baseline condition. Similar interventions have been applied to each of the remaining $4$ variables shown in Table 1 (PKC, PIP2, Mek, and PIP3) which have been intervened upon in a different experimental condition.

We hope that this explanation makes the setting more clear; we will expand appropriately the description of the Sachs dataset in the improved manuscript.

2. Interpretability

We agree that interpretability is a valuable property of any causal model. Like many statistical models, the full causal structural model may contain information that can overwhelm a human researcher, more so when it includes latent variables. Below we give some pointers on how to interpret their meaning.

• Interpretation of mixing variables:
  The categorical nature of the mixing variables has the benefit that a domain expert can directly inspect the components and their number. For example, someone interested in a specific target node $Y$ with cause $X$ and affected by a latent $Z$ can readily visualize the $k$ components of $Z$ as "clusters" and inspect how the causal mechanism $X\rightarrow Y$ differs in the respective groups, similar to in Fig. 1 of the manuscript.
  When $Y$ has more than one cause, a pair-wise analysis or a dimensionality reduction technique could allow for a similar visual analysis. It is also interesting to see which nodes $Z$ jointly influence a subset of nodes, as this is a strong indication of latent subgroups influencing a part of the system. When sufficient background domain knowledge is available, one could match the $k$ components to known genetic subgroups or biological states. If no such expert knowledge is available, one might gain novel insights into the existence of possible causal mechanisms that were previously unknown.

• Follow-up work:
  Your question inspires some ideas for future work where we further interpret the meaning of a given $Z$, especially when they point to previously unknown states. To this end, one could explore explainability approaches such as LIME or Shapley values. For this, one could see the latent mixing variables $Z$ as classification targets, and then test whether any known or background knowledge features are important predictors of $Z$ using LIME or Shap. For example, one might find whether the latent components arise from a population group of a certain age, demographic, or other background.

Please let us know if we can address any remaining questions or concerns.

We appreciate your time and thank you for the detailed feedback.

Please consider our responses to your points below.

Response to Questions

1. Source Variables

This is a good question; we will clarify the meaning of edges such as $Z_1 \rightarrow X_1$ (Fig. 2), whenever $X_1$ is a source node of the causal graph.
Note that the GMM used in the sources can be seen as a special case of the MLR model if we use an implicit input variable with a constant value of 1. Indeed, re-visiting Eq. (4) we would then have $\mathbf Y= \{\mathbb 1\}$, where $\mathbb 1$ is a deterministic input variable $\mathbb 1 = 1$.
Then, we would obtain the GMM as $p_{X|\{\mathbb 1\}}^{\mathrm{MLR}}(x;\mathbf\beta,\mathbf\gamma,\sigma^2)=\sum_{k=1}^K\gamma_k p^{\mathcal N}_{X}(x;\beta_k\cdot 1,\sigma^2),$ which the reader can easily verify to be the p.d.f. of the GMM model.

This also gives a seamless approach: We use EM and the BIC to infer latent mixing for both sources and non-sources, where we fit a GMM in the case of sources and an MLR otherwise.

We will include this explanation after Eq. (4) to be clear on its meaning for source nodes.

2. Equal Noise Variances and Identification

Indeed, methods exist to identify linear additive Gaussian noise models with equal variance (under further conditions, e.g., non-perfect dependencies), and what is more, also when other constraints are added (e.g., that the coefficients are all less than the unit). Our identification results do not rely on such assumptions. Instead, we base identification on the effects of the latent mixing variables and the asymmetries they introduce in the causal directions.
Thus, our results can be seen as orthogonal to the approach of restricting the functional model with further assumptions, and can therefore address scenarios where these perhaps unrealistic assumptions do not hold.

3. Formulation of Theorem 3.4

• (a) We thank you for spotting this typesetting error. Here the correct distribution is $X|\mathbf Y$, in alignment with the notation of Eq. (4).
  We have made the correction $\mathbf X\to X$ and $Y\to \mathbf Y$ in the improved manuscript.
  The intended meaning is as you correctly assumed.

• (b) We thank you for and agree with this recommendation. Following your suggestion, we refrain from using 'Markov Equivalence Class' here to avoid confusion.

• (c) Indeed, the two are the same. We have corrected the theorem numbering in the updated appendix.
  As for Lemma 3.3, we agree and will make the appropriate change.

4. Specifying the Number of Mixtures

Exactly, the BIC in lines 210-212 accounts for an unknown number of mixtures. That is, we rely on the (theoretical) latent-aware BIC as in Eq. (6), where the likelihood takes the form of Eq. (4).
To compute the BIC in practice, we need to consider each 1\leq K\leq K_\max and use the EM algorithm to obtain estimates for the MLR parameters; finally, we then use the BIC score of the best such $K$ in our algorithm (and the corresponding estimates for the rest of the model parameters). We updated our manuscript text to make this explicit.

We are happy to dispel any remaining doubts and clarify any persisting source of confusion.

I thank the authors for their response. I would like to maintain my current evaluation of the paper. I have two minor comments:

1. Regarding Question 1: I recommend that the authors separate the bias term from the vector product in $f$—for example, by introducing an explicit bias term (e.g., $\beta_{j,k}^{(0)}$) instead of setting $\mathbf{Y}=1$.

   My understanding is that $\mathbf{Pa}_j$ is a subset of $\mathbf{X}$, and for $X_1$, we have $\mathbf{Y}=\mathbf{Pa}_1=\emptyset$. Please correct me if I am wrong.

2. Regarding Question 3(a): I suggest replacing $\mathbf{Y}$ with $\mathbf{Pa}$, or explicitly referring to Equation (3) or (4) to avoid confusion.

We thank you for your response and appreciate your time in the review process.

Please consider our response to your comments below.

Response to Questions

• Consistency of Alg. 1.
  Indeed, studying the consistency of the estimates derived by the EM algorithm is a difficult problem, and the literature is limited to rather strict settings (see the references in our manuscript). Therefore, it is correct that we cannot easily provide strict guarantees for Alg. 1.
  Having mentioned this, we can further note that the experimental evidence supports the connection we make in our conjecture. There is also a reasonable mechanism that could theoretically explain the observed behaviour (asymptotically).
  Simply put, when the data truly follow the MLR model, it is more likely for our EM implementation to align with the oracle; when this does not happen, it is likely due to a faulty model, in which case our downstream process would anyways (asymptotically) reject the specific MLR hypothesis in favour of the correct alternative.
  In summary, we confirm your observation that Alg. 1 is only guaranteed to be consistent under appropriate conditions, such as a good-enough EM-initialization regime, and should otherwise be viewed as a reasonable approximation to the problem.

• Choice of $K_\text{max}$.
  This is a fair point; indeed, we assume that a reasonably high maximum number of $K_\text{max}$ is given that caps the number of mixture components.
  Here, the following comments may be made.

  • First, even when $2 \leq K_\text{max}<K*$, for $K*$ the true number of clusters, the causal structure would still be (likely) correctly inferred, as the capped MLR would still offer a higher likelihood versus any alternative model (under appropriate assumptions, that is, non-zero additive noise and that the total back-door noise being higher than the remaining variance due to the mixture collapse.)
  • Further, if the MLR model ever selects $K < K_\text{max}$, it is (asymptotically) sure (under the oracle assumption, assuming the latter to be capped, too) that this estimate of $K$ is the true number of mixture coefficients. When $K=K_\text{max}$, one could keep increasing the bound and repeat. Hence, we can consider this bound as more of an implementation issue than a truly limiting factor.
  • For practitioners, we consider small values of $K_\text{max}$ (in the order of 10 or less) sufficient, with motivational examples in mind where $K$ counts distinct biological states, subpopulations, or conditions with distinct causal mechanisms, so that very large $K$ are most likely not meaningful. We recommend that users inspect the BIC values across different $K$ and increase the maximum if necessary.

Feel free to let us know whether we sufficiently addressed the questions and if any concerns remain.