Bivariate Causal Discovery with Proxy Variables: Integral Solving and Beyond

We appreciate your efforts and valuable suggestions in reviewing our paper. We address your concerns below.

Q1. The method relies on a partially known prior structure, specifically that $U$ causes both $X$ and $Y$, and that $Z$ and $W$ serve as proxy variables for $U$. This assumption may be restrictive, as in practice, we may not have direct knowledge of the existence of $U$ or whether $Z$ and $W$ are valid proxies.

A. This is a standard setting in proximal causal inference. Indeed, since the observational causal inference/discovery typically involves unobserved measurements, it is essential to introduce substitution variables, such as instrumental variables and proxy variables, as considered here. In many scenarios, these assumptions naturally hold and such proxy variables are easy to obtain, for example, they can be noisy measurements of the confounders [1] or as descendants in time-series data [2].

Q2. The asymptotic properties of the proposed method depend on a set of assumptions, which may limit its applicability in real-world scenarios.

A. The most key assumption is completeness, which is standard in proximal causal inference and can easily hold in our scenario. Others are regularity conditions that are also standard for kernel regression.

Q3. The uniqueness of the solution (not only the existence) to the integral equation seems crucial for the hypothesis test. Could the authors provide an intuitive explanation or illustration of its role?

A. As claimed in lines 168-176 (right column), uniqueness is not required in our procedure, since our goal is to determine whether a solution exists. Among all solutions to the integral equation, our estimate converges to the least-norm solution (Lemma D.19 in Appendix D.5).

Q4. If $U$ is a vector rather than a single variable, does the proposed theoretical result still hold?

A. Our procedure remains valid even when $U$ is a high-dimensional variable, but requires the proxy variable $W$ to satisfy the completeness condition (Assumption 4.1). It implies that the dimension of $W$ is greater than that of $U$. When this holds, the completeness is easy to hold.

Q5. The theoretical contribution is intriguing! I am curious whether the proposed method can be viewed as a generalization of Tetrad constraints when considering $Z$ and $W$ as proxy variables.

A. Thank you for your insightful question. Our methods can be viewed as a generalization of Tetrad constraints. Specifically, the tetrad constraint was introduced to test whether $X, Y, Z, W$ are conditionally independent given $U$, under linear models. Specifically, the constraint is:

$$\mathrm{cov}( X,Y ) \mathrm{cov}( Z,W )= \mathrm{cov}( X,Z ) \mathrm{cov}( Y,W )= \mathrm{cov}( X,W ) \mathrm{cov}( Y,Z ).$$

When proxies $W,Z$ are assumed to satisfy $W \perp Y|U$ and $Z \perp X|U$, the constraint degenerates to:

$$\mathrm{cov}( X,Y ) \mathrm{cov}( Z,W )= \mathrm{cov}( X,W ) \mathrm{cov}( Y,Z ),$$

and it can be used to test $\mathbb{H}_0: X \perp Y|U$. In contrast, our constraint for testing $\mathbb{H}_0$ is formulated by integral equations (1) and (10). Under the linear Gaussian setting with standard normal exogenous noise, our constraint gives rise to the one in #. Specifically, suppose $U \sim \mathcal{N}(0,1)$, $W=\mu_UU +\varepsilon_W$, $Z=\beta_UU +\varepsilon_Z$, $X=\alpha_UU +\varepsilon_X$, and $Y=\gamma_UU+\varepsilon_Y$, where $\varepsilon_W,\varepsilon_Z, \varepsilon_X, \varepsilon_Y$ are standard normal. We have $h(W) = \frac{\gamma_U}{\mu_U}W$. Since $\mathrm{cov}( Z,W ) = \beta_U \mu_U$ and $\mathrm{cov}( Z,Y ) = \beta_U \gamma_U$, we have $\mathrm{cov}( X,Y ) \mathrm{cov}( Z,W )=\gamma_U \alpha_U \mu_U \beta_U = \mathrm{cov}( X,W ) \mathrm{cov}( Y,Z )$.

[1] Kuroki, M. and Pearl, J. Measurement bias and effect restoration in causal inference. Biometrika, 101(2):423-437, 2014.

[2] Liu, M., Sun, X., Hu, L., and Wang, Y. Causal discovery from subsampled time series with proxy variables. Advances in neural information processing systems, 2023.

Thank you for the positive assessment and valuable suggestions about our paper. We address your questions below. We will correct typos and make the number of lemmas consistent in the updated version.

Q1. How sensitive is PMCR to kernel and bandwidth choices?

A. We would like to clarify that PMCR requires kernels to be bounded, continuous, and integrally strictly positive definite (Assumptions C.2 and C.3). Typical types of kernels satisfying these conditions include Gaussian RBF and Laplacian kernels. We adopt the Gaussian kernel because it is commonly used in the literature. Its bandwidth is typically initialized via the median distance heuristic, ensuring robust performance.

Q2. What is the computational efficiency of the proposed approach? How well does it scale to larger datasets and/or high-dimensional variables?

A. PMCR has a time complexity of $O(pn^3)$, where $n$ is the sample size and $p$ is the dimension of $W$. The complexity is linearly proportional to $p$ because, for multivariate $W$, the kernel can be constructed as the product of scalar kernels applied to each input dimension $w$, i.e., $k(w, w') = \prod_{j=1}^{p} k(w[:, j], w[:, j]')$, where $w[:, j]$ denotes the $j$-th feature. The complexity with respect to $n$ can be reduced using standard Cholesky techniques when computing the inverse. We record the running time of a single trial as $n$ varies (when $W$ is univariate). As shown, the actual running time increases much more slowly than $O(n^3)$.

Q3. Are the non-identifiability results (Prop. 5.1) generalizable beyond linear models?

A. We believe that nonlinear models also suffer from the non-identifiability issue when $W$ has a strong effect on $Y$, but the theoretical analysis is challenging as there is no closed-form solution for non-linear models.

To empirically verify this point, we consider a non-linear setting: $U \sim \mathcal{N}(0,1)$, $W=U +\varepsilon_W$, $X=U +\varepsilon_X$, and $Y=X^2+\gamma_W W^2+U+\varepsilon_Y$, where $\varepsilon_W, \varepsilon_X, \varepsilon_Y$ are standard normal. We generate $n=800$ samples and record the average power over 20 times. We can see that the power decreases as $\gamma_W$ increases, which suggests the difficulty of identifying the causal relation when $W$ strongly affects $Y$.

We appreciate your efforts and suggestions in reviewing our paper. We address your concerns below.

Q1. About the type I error level in Liu's paper.

A. Liu's type-I error control requires a Lipschitz smooth function, but sqrt does not satisfy this. To make it hold, we conduct experiments with the sigmoid function. As shown, Liu's method approximates the type-I error as $n$ increases. It is also seen that our method has better sample efficiency than theirs.

Q2. No evidence to support sample efficiency.

A. The table above shows our method's superior sample efficiency in type-I error control. As claimed in lines 107-109, this is because Liu's method requires the bin number to go to infinity to make the discretization error vanish.

Q3. The power does not approach to one in Fig. 3b.

A. We further conduct experiments as $n$ increases. Combined with Fig. 3b, they show a clear trend of power approaching one. When $n=4000$, the power is 0.97.

Q4. About definitions and claims.

A. About $\hat{H}^\lambda(w,t)$, its relation to the paragraph before "equivalent speaking". As mentioned in lines 201-211, $\hat{H}^\lambda(w,t)$ is the minimizer of (5), which is designed to solve the least-norm solution as introduced in lines 176-181. From lines 182-208, we explain how the risk (5) is derived. We first introduce the risk in lines 192-193, followed by its equivalent form in lines 197-198. Indeed, the risk (5) is its empirical version, augmented with a regularization term. Additionally, we provide a more detailed explanation in Appx. C.4-C.5 due to space limitations. We will offer further clarifications in the updated version.

About the motivation of Tikhonov regularization. The regularization provides a way to solve ill-posedness problems and is commonly adopted in the literature of kernel regression. We have explained in more detail in Appx. C.5.

About the bridge function. The bridge function is the solution to the integral equation used in proximal causal inference, serving a similar role to $H$ in (4). It is a fundamental concept in proximal causal inference, and we have provided references in lines 167-168 due to space limits.

Q5. About power analysis after successive approximation.

A. Our power analysis focuses on alternatives in (4), which are defined based on characteristic restrictions. By utilizing the characteristic function, our approach demonstrates significantly better power compared to first-order moment methods in the literature (Fig. 4).

Q6. Does (4) hold if (1) does not hold under the assumptions in Thm 4.5?

A. We are sorry that we're unsure if we fully understand the question. Under assumptions in Thm. 4.5, both (1) and (4) should hold, as (1) is the conclusion of Thm. 4.5, and (4) naturally follows if (1) holds. We would appreciate your clarification if we've misunderstood.

Q7. Does assump. 6.1 say that all variability in U is captured by X? That seems to make a "confounder" immaterial.

A. Yes, it means variability in $U$ is captured by $X$. However, we want to clarify that this assumption is easy to hold and will not make $U$ immaterial.

First, the assumption applies to a wide range of models, as long as the dimension of $X$ is greater than that of $U$. Under this condition, [Andrews et al. 2017] demonstrated that completeness generically holds. Even if this condition holds, the confounding problem still matters as it also causes bias since completeness does not imply that $p(x)$ would determine (up to transformation) $p(x,u)$ (this is known as the ``equivalence condition" in Miao et. al. 2023) and hence the distribution of $p(u)$. In this regard, we still cannot determine whether $p(y|do(x)) \neq p(y|x)$ when $U$ is unobserved.

Q8. About "there exist" statement in $H^0(w,t)$" (line 234).

A. Local alternative refers to alternatives that are very close to the null hypothesis. That means, it adds a small deviation $r(X)/n^\alpha (\alpha > 0)$ to some $H^0$. To ensure that it does not degenerate to $\mathbb{H}_0$, we require $r(X)/n^\alpha$ cannot be written as $E(H-H^0|X)$ for any $H$, as claimed in lines 239-241.

Q9. Why need h(w,y) to be square-integrable?

A. If $h(w, y)$ is not square-integrable, it cannot be solved using regression methods, as regression approaches minimize the squared loss that involves the second-order moment of the solution. Besides, our method and analysis are built upon the completeness assumption, which is fundamental in proximal causal inference and is imposed on the square-integrable function class.

Q10. What does "feasibility of solutions" mean in Line 204?

A. It refers to the capability of regression methods to achieve consistency for integral solving (Lines 200-203).