We thank you for your thoughtful and constructive feedback. We are especially grateful for the recognition of the clarity of our writing, the novelty of our theoretical contribution, and the empirical improvements in efficiency without sacrificing accuracy. We hope the following responses properly address your concerns.

Q1. In your definition for a maximal ancestral graph, is that the same as an acyclic directed mixed graph (ADMG)? Further, what is an “almost directed cycle” in lines 90-91?

A1. Thank you for the question.

To clarify: maximal ancestral graphs (MAGs) are a subclass of acyclic directed mixed graphs (ADMGs). MAGs are used to represent the structure of data-generating processes in the presence of latent variables [Richardson, 2003; Spirtes et al., 1997] (See the first sentence of the third paragraph of [Richardson, 2003]), and they satisfy additional constraints that ADMGs do not. For example, in MAGs, every missing edge corresponds to a conditional independence relation, making them maximal in that sense.

Regarding the term ``almost directed cycle'', it refers to a scenario where a node $X$ is both a spouse and an ancestor of another node $Y$. Specifically, $X$ is considered a spouse of $Y$ if there exists a bi-directed edge $X \leftrightarrow Y$, and an ancestor if there is a directed path from $X$ to $Y$. In an ancestral graph, the presence of a bi-directed edge $X \leftrightarrow Y$ implies that neither $X$ is an ancestor of $Y$ nor $Y$ of $X$. Therefore, allowing both a bi-directed edge and a directed path between two nodes would violate the ancestral property and create an "almost directed cycle," which is inconsistent with the semantics of MAGs [Richardson et al., 2002].

Reference

• Richardson T. Markov properties for acyclic directed mixed graphs. Scandinavian Journal of Statistics, 2003.
• Spirtes P, Richardson T, Meek C. The dimensionality of mixed ancestral graphs. Technical Report CMU-PHIL-83, 1997.
• Richardson T, Spirtes P. Ancestral graph Markov models. The Annals of Statistics, 2002.

Q2. Theorems 2,3, and 4 are nearly exactly the same as the inference rules and Theorems 2 and 3 presented in Entner et al. (2013). Is there a specific reason why the authors did not cite Entner et al. (2013) when presenting these three theorems? It seems to me that the proofs of Theorems 2, 3, and 4 should follow in a fairly straightforward manner from the proofs in Entner et al. (2013) once Theorem 1 is established. I believe that it is important to cite Entner et al. (2013) when presenting these three theorems as they are, at the very least, heavily inspired by the work of EHS.

A2. Thank you for pointing this out.

Theorems 2, 3, and 4 are indeed heavily inspired by the inderence rules proposed in Entner et al. (2013). In the revised version, we have added appropriate citations to Entner et al. (2013) when presenting these theorems.

We would also like to highlight a key distinction: our formulation is fully local, relying only on the Markov Blankets of the treatment and outcome variables, rather than performing an exhaustive search over all variable combinations, as in the original inference rules. Importantly, our local approach is not merely a heuristic: it is theoretically guaranteed to be both sound and complete under the stated assumptions. This localization significantly improves computational efficiency and enhances the method's applicability in high-dimensional settings. We believe this constitutes a novel and meaningful extension beyond the original work.

Q3. Is faithfulness implicit in the theorems of this paper? Line 257 is the first mention of the faithfulness assumption in the main paper. If the authors are assuming faithfulness throughout, then they should define faithfulness and state the assumption somewhere in Section 2.

A3. Yes, the faithfulness assumption is implicit in our theoretical results. Following your suggestion, we have moved its definition and formal statement from Appendix B to Section 2 of the main text for greater clarity.

We would also like to emphasize that faithfulness is a standard assumption in causal discovery methods. In fact, many classical algorithms—such as PC, FCI, and RFCI—fundamentally rely on faithfulness [Spirtes and Zhang, 2016]. Moreover, existing covariate selection methods, including EHS [Entner et al., 2013], CEELS [Cheng et al., 2022], and LDP [Maasch et al., 2024], also make use of this assumption to ensure the correctness of their independence-based inference procedures.

Further, earlier work shows that reducing unnecessary conditional independence tests can help mitigate the impact of statistically weak violations of faithfulness [Isozaki T, 2014]. While faithfulness remains a theoretical requirement, our method’s localized design naturally reduces the number and complexity of conditional independence tests, which empirically improves robustness to potential violations of this assumption.

Reference

• Peter Spirtes and Kun Zhang. Causal discovery and inference: concepts and recent methodological advances. In Applied informatics, 2016.
• Isozaki T. A robust causal discovery algorithm against faithfulness violation. Information and Media Technologies, 2014.

We sincerely thank you for the thoughtful comments and for recognizing that our paper is well written, with clear remarks and discussion. Below, we address each point raised:

W1. The paper is notationally heavy. I am not sure there is much the authors can do.

A1. We acknowledge that the notation may be dense at times. This is largely due to the technical nature of the identification problem under latent confounding. To assist readers, we have included a notation table (Table 1) in Appendix B. In the final version (thanks to the availability of an additional page), we will move this table to the main text to improve accessibility and readability.

W2. The work realies on stronger assumption about ordering of covariates and knowledge of MAG/DAG. I am not sure if it is realistic.

A2. We would like to clarify that our method does not assume access to the true MAG or DAG. Instead, it is a local, data-driven approach that operates without relying on global structural knowledge. All experimental results are based solely on observed data, without any access to the underlying MAG or DAG.

Regarding the pretreatment assumption, it realistically reflects how data are collected in many fields such as economics and epidemiology [Hill, 2011, Imbens and Rubin, 2015, Wager and Athey,2018], where all covariates in $\mathbf{O}$ are measured before treatment and outcome. This assumption is particularly common in observational studies [Cheng et al., 2022, Entner et al., 2013, De Luna et al., 2011, Vander Weele and Shpitser, 2011, Wu et al., 2022]. For example:

• In medical studies, covariates such as age, gender, and past medical history are collected before administering a treatment and the outcome;
• In economic policy evaluation, covariates like income, education level, and demographic data are collected before a policy is implemented or the downstream outcome of interest.

We adopt this assumption to focus on covariate selection without requiring causal sufficiency, which is a strong and often unrealistic assumption in latent variable settings.

Moreover, LSAS outperforms other methods even when the pretreatment assumption may not be satisfied, as demonstrated by results on benchmark graphs (Figures 4 and 11).

Q3. Can authors provide sensitivity analysis to violation of correct specification of MAG/DAG?

A3. We would like to clarify that our method does not assume access to the true MAG or DAG. LSAS is a fully local and data-driven procedure that identifies valid adjustment sets without requiring knowledge of the global causal structure. As such, there is no dependency on the correct specification of the full MAG/DAG, and hence no direct source of sensitivity to its misspecification.

Q4. In Algorithm 1, what if R1 or R2 is satisfied for multiple different combinations of S and Z? What if they result in different estimates?

A4. Thank you for this valuable question.

When multiple $(S, \mathbf{Z})$ pairs satisfy condition $\mathcal{R}1$, we follow the strategy of [Cheng et al., 2022] and average the resulting causal effect estimates across these valid adjustment sets. This approach helps improve stability and reduce variance.

In contrast, if any $(S, \mathbf{Z})$ pair satisfies condition $\mathcal{R}2$, this implies that there is no causal effect from $X$ to $Y$. In such cases, LSAS returns a null estimate, consistent with the interpretation of R2 as a test for the absence of a causal effect.

Reference

• Debo Cheng, Jiuyong Li, Lin Liu, Jiji Zhang, Jixue Liu, and Thuc Duy Le. Local search for efficient causal effect estimation. IEEE Transactions on Knowledge and Data Engineering, 2022.

Q5. Typically, MB(Y) is always much larger than MB(X). In that case, their algorithm thats exponentially more time to run. How are authors claiming scalability?

A5. We would like to clarify that MB(Y) is not necessarily much larger than MB(X) in general. The relative sizes of Markov Blankets depend heavily on the structure of the underlying graph.

To demonstrate the empirical scalability of LSAS, we conducted additional experiments on denser random graphs, reporting both runtime and average conditioning set sizes. Specifically, we conducted additional experiments on random graphs with increasing average degrees (3, 5, 7, and 9), under the same settings as in the main paper (40 nodes, sample size = 5K). These denser graphs result in larger MBs, which present a greater computational challenge.

As shown in the table below, we find the follwing observations:

1. LSAS consistently outperforms all baseline methods in terms of RE and nTests, even as graph density increases. One exception is the (40, 9) network, where the local method LDP has lower nTest. While performance degrades with increasing density—as expected for all methods—LSAS maintains notable advantages, particularly in RE accuracy.

2. The runtime of LSAS is slightly higher than that of LDP in some denser settings; however, its RE is substantially lower, reflecting its completeness in identifying VAS. This is because we search for all adjustment sets within the Markov Blanket of $Y$. Limiting the number of valid adjustment sets found before stopping would help speed up the running time.

3. As graph density increases, the Markov Blanket of each node grows, leading to increased computational complexity for LSAS. This is reflected in the observed runtimes and conditioning set sizes. Importantly, this trend affects all methods, not just LSAS.

In summary, while LSAS shares the worst-case exponential complexity inherent to constraint-based methods, its localization strategy yields strong empirical efficiency, especially on sparse-to-moderate graphs. We acknowledge that performance degrades with graph density, but LSAS continues to offer better RE accuracy and competitive efficiency compared to all baselines, even in denser settings.

Q6-Limitations. The authors discuss limitations in appendix but I think it is worth moving to maintext. Further, it will be useful to discuss the reliance on assumptions and what happens when they are violated.

A6. Thank you for the suggestion. Due to page constraints, we originally placed the discussion of limitations in the appendix, but we will move it to the main text in the final version as recommended.

Regarding assumptions, our method relies on the pretreatment assumption, which, as discussed in Response A2, is widely adopted and considered realistic in many observational studies across fields such as medicine and economics.

When the pretreatment assumption is violated, some observed covariates may be descendants of the treatment, potentially lying on causal paths from $X$ to $Y$. In such cases, even if $(S, \mathbf{Z})$ satisfies condition $R1$, adjusting for variables in $\mathbf{Z}$ may introduce post-treatment bias and invalidate the causal estimate. For example, in the graph:

if we choose $S = F$ and $\mathbf{Z} = \{B, D\}$, then $(S, \mathbf{Z})$ satisfies condition $R1$, but $\mathbf{Z}$ violates the generalized adjustment criterion since $B$ is a descendant of $X$.

We will incorporate this discussion into the revised version to more clearly acknowledge the method’s assumptions and limitations.

We sincerely thank you for your positive evaluation and thoughtful questions. We greatly appreciate your recognition of the novelty, theoretical foundation, clarity of presentation, and comprehensive experimental results in our work. We hope that the following responses adequately address your concerns.

Q1. I still cannot understand why the pretreatment assumption is realistic from the example "every variable within the set is measured prior to the implementation of the treatment and before the outcome is observed". Could you please provide more specific explanations?

A1. Thank you for pointing this out. As discussed in Remark 1 of the main paper, the pretreatment assumption is realistic as it reflects how data are collected in many fields such as economics and epidemiology [Hill, 2011, Imbens and Rubin, 2015, Wager and Athey,2018], where all covariates in $\mathbf{O}$ are measured before treatment and outcome. This assumption is particularly common in observational studies [Cheng et al., 2022; Entner et al., 2013; De Luna et al., 2011; Vander Weele and Shpitser, 2011; Wu et al., 2022]. For example:

• In medical studies, covariates such as age, gender, and past medical history are collected before administering a treatment and the outcome;
• In economic policy evaluation, covariates like income, education level, and demographic data are collected before a policy is implemented or the downstream outcome of interest.

In such a setting, covariates are guaranteed not to be affected by the treatment or outcome, justifying the assumption. We will clarify this with concrete examples in Section 2.3 of the revised manuscript.

Reference

• Jennifer L Hill. Bayesian nonparametric modeling for causal inference. Journal of Computational and Graphical Statistics, 20(1):217–240, 2011.
• Guido W Imbens and Donald B Rubin. Causal inference for statistics, social, and biomedical sciences: An introduction. Cambridge University Press, 2015.
• Stefan Wager and Susan Athey. Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523):1228–1242, 2018.
• Debo Cheng, Jiuyong Li, Lin Liu, Jiji Zhang, Jixue Liu, and Thuc Duy Le. Local search for efficient causal effect estimation. IEEE Transactions on Knowledge and Data Engineering, 2022.
• Doris Entner, Patrik Hoyer, and Peter Spirtes. Data-driven covariate selection for nonparametric estimation of causal effects. In Artificial intelligence and statistics, pages 256–264. PMLR, 2013.
• Xavier De Luna, Ingeborg Waernbaum, and Thomas S Richardson. Covariate selection for the nonparametric estimation of an average treatment effect. Biometrika, 98(4):861–875, 2011.
• Tyler J Vander Weele and Ilya Shpitser. A new criterion for confounder selection. Biometrics, 67(4): 1406–1413, 2011.
• Anpeng Wu, Junkun Yuan, Kun Kuang, Bo Li, Runze Wu, Qiang Zhu, Yueting Zhuang, and Fei Wu. Learning decomposed representations for treatment effect estimation. IEEE Transactions on Knowledge and Data Engineering, 35(5):4989–5001, 2022.

Q2. While the paper mentions experimental effectiveness even when the pretreatment assumption is violated, could the authors elaborate on the theoretical implications or provide a more detailed empirical analysis of LSAS's performance in such scenarios?

A2. Thank you for the insightful question. When the pretreatment assumption is violated, some observed covariates may be descendants of the treatment, potentially lying on causal paths from $X$ to $Y$. In such cases, even if $(S, \mathbf{Z})$ satisfies condition $R1$, adjusting for variables in $\mathbf{Z}$ may introduce post-treatment bias and invalidate the causal estimate. For example, in the graph:

if we choose $S = F$ and $\mathbf{Z} = \{B, D\}$, then $(S, \mathbf{Z})$ satisfies condition $R1$, but $\mathbf{Z}$ violates the generalized adjustment criterion since $B$ is a descendant of $X$. A potential workaround is to identify descendants of $X$ first, then apply $R1$—but this lacks completeness, and may fail to recover adjustment sets even when they exist. To the best of our knowledge, no existing method can soundly and completely identify descendants of a treatment variable locally, especially in the presence of latent variables and without recovering the full graph. Addressing this remains an open and non-trivial challenge.

Q3. The first step of LSAS is learning the Markov Blankets. What specific MB discovery algorithms are used? What are the computational complexities and known practical limitations of these MB discovery methods that might impact the overall efficiency of LSAS?

A3. Thank you for your question. In our implementation, we use the TC (Total Conditioning) algorithm [Pellet and Elisseeff, 2008b] to discover the MB; see Appendix C.2 for further details. The computational complexity of this method is $O(n)$, where $n$ is the number of observed variables, making it efficient for use as a first-stage procedure in LSAS.

It is well known that the accuracy of MB discovery directly affects the overall performance of MB-based causal inference methods. In our setting, if relevant variables are omitted from the estimated MB, LSAS may fail to identify a VAS. To mitigate this risk, one can adopt a conservative approach: relaxing the inclusion criteria to ensure that as many true MB variables as possible are retained [Wang et al., 2014; Liu et al., 2020]. This increases the chance of including extraneous variables beyond the true MB, which in turn may increase the computational burden in the subsequent search phase. Nonetheless, when runtime is not a primary concern, this trade-off improves robustness and helps preserve the completeness of the adjustment set search.

Reference

• Jean-Philippe Pellet and André Elisseeff. Using Markov blankets for causal structure learning. Journal of Machine Learning Research, 9(7), 2008b.
• Wang, C., Zhou, Y., Zhao, Q., & Geng, Z. Discovering and orienting the edges connected to a target variable in a DAG via a sequential local learning approach. Computational Statistics & Data Analysis, 2014.
• Liu, Y., Fang, Z., He, Y., Geng, Z., & Liu, C. Local causal network learning for finding pairs of total and direct effects. Journal of Machine Learning Research, 2020.

We thank you for the thoughtful feedback and appreciate the positive remarks on the soundness and originality of our work. We address your concerns below.

W1&Q1. the pretreatment assumption...generalization to non-pretreatment settings...fail without this assumption...beyond m-structures LSAS aims to exclude.

A1. We would like to clarify the motivation and scope of our work as follows.

• First, regarding the pretreatment assumption, it realistically reflects how data are collected in many fields such as economics and epidemiology [Hill, 2011, Imbens and Rubin, 2015, Wager and Athey,2018]. We adopt this assumption to focus on covariate selection without requiring causal sufficiency. While the full set $\mathbf{O}$ forms a VAS in the absence of latent confounders, conditioning on all covariates—especially in high-dimensional settings—can lead to poor performance [Luna et al., 2011; Abadie and Imbens, 2006]. Our goal is to identify the local sufficient VAS that enables unbiased and efficient causal effect estimation, even in the presence of latent variables.

• Second, even under the pretreatment assumption, selecting a VAS remains non-trivial in the presence of latent confounders. Not all pretreatment covariates are harmless to condition on—for example, variables involved in m-structures or other types of collider paths (e.g., $X \leftrightarrow A \leftrightarrow B \leftrightarrow C \leftrightarrow Y$) may introduce bias by opening backdoor paths. Our method identifies which covariates to include or exclude based on their roles in the local causal structure, which is learned directly from data without requiring global graph recovery.

• Third, existing methods-such as EHS, CEELS, and ours-do not directly extend to settings where the pretreatment assumption is violated, as they may select invalid adjustment sets that include descendants of the treatment, leading to biased estimates. For example, in the graph:

if we choose $S = F$ and $\mathbf{Z} = \{B, D\}$, then $(S, \mathbf{Z})$ satisfies condition $R1$, but $\mathbf{Z}$ violates the generalized adjustment criterion since $B$ is a descendant of $X$. A potential workaround is to identify descendants of $X$ first, then apply $R1$—but this lacks completeness, and may fail to recover adjustment sets even when they exist. To the best of our knowledge, no existing method can soundly and completely identify descendants of a treatment variable locally, especially in the presence of latent variables and without recovering the full graph. Addressing this remains an open and non-trivial challenge.

• Finally, as discussed in Section 1 and Appendix A, existing global adjustment set methods are often inefficient, while existing local approaches suffer from incompleteness—even under the pretreatment assumption. This highlights the need for a method that is both sound and complete, operates locally, and remains valid in the presence of latent variables, to reliably select the VAS for specific causal queries.

1. Abadie A, Imbens G W. Large sample properties of matching estimators for average treatment effects. Econometrica, 2006.

W2&Q2. ...report runtime metrics alongside nTests... Analyze conditioning sets...theoretical or empirical analysis of how conditioning set sizes scale with problem dimensions?

A2. Thank you for the thoughtful comments.

We agree that nTests alone may not fully reflect runtime performance, especially when conditioning set sizes differ across methods. As suggested, we have now included both wall-clock runtimes and average conditioning set sizes for all benchmark networks. Due to high computational cost, EHS did not finish within 2 hours; we report its runtime and conditioning set sizes using a 200-second early stopping threshold for fair comparison.

Our key observations are as follows:

1. EHS, as a global search method, exhibits the highest runtime and largest conditioning sets across all settings, as expected.
2. LSAS consistently outperforms baselines in runtime across most graphs and sample sizes. One exception is the WIN95PTS network at large sample sizes (15K), where the local method LDP is faster; however, LSAS achieves substantially higher RE accuracy (see Figures 4 and 11), due to LDP’s incompleteness in identifying valid adjustment sets.
3. Among the three local search methods, LSAS uses slightly larger conditioning sets on average. This is because LSAS is designed to search for all valid adjustment sets within $MB(Y)$, rather than stopping after identifying the first one. Despite this overhead, LSAS achieves consistently lower RE than CEELS and LDP (see Figures 4 and 11), highlighting the advantage of its completeness in VAS identification.

In summary, these results suggest that LSAS offers a favorable trade-off between runtime and estimation accuracy.

W3&Q3. ...the performance of LSAS in denser graphs...how LSAS complexity scales with graph density?...characterize the settings where LSAS provides advantages?

A3. Thank you for the valuable suggestion. We conducted additional experiments on random graphs with average degrees 3, 5, 7, and 9 (40 nodes, 5K samples), as in the main paper. Due to high computational cost, EHS did not finish within 2 hours; we report its runtime and conditioning set sizes using a 200-second early stopping threshold for fair comparison.

Our key observations are as follows:

1. As shown in the table below, LSAS consistently outperforms all baseline methods in terms of RE and nTests, even as graph density increases. One exception is the (40, 9) network, where the local method LDP has lower nTest. While performance degrades with increasing density—as expected for all methods—LSAS maintains notable advantages, particularly in RE accuracy.
2. The runtime of LSAS is slightly higher than that of LDP in some denser settings; however, its RE is substantially lower, reflecting its completeness in identifying VAS. This is because we search for all adjustment sets within the Markov Blanket of $Y$. Limiting the number of valid adjustment sets found before stopping would help speed up the running time.
3. As graph density increases, the Markov Blanket of each node grows, leading to increased computational complexity for LSAS. This is reflected in the observed runtimes and conditioning set sizes. Importantly, this trend affects all methods, not just LSAS.

In summary, while LSAS shares the worst-case exponential complexity of constraint-based methods, its localization strategy yields strong empirical efficiency, especially on sparse-to-moderate graphs. While performance degrades with increasing density—as with all methods—LSAS consistently delivers better RE accuracy and competitive efficiency across all settings.

I thank the authors for their thoughtful response and for addressing the points I raised. I will increase my score to weak accept (4), and am inclined to raise it further if the authors agree to incorporate the additional experimental results present in answers A2 and A3, as well as the discussion in A1, into the camera ready version of this paper. Such results are crucial to contextualize the strong performance of LSAS, relative to baselines, across different contexts.