A Conditional Independence Test in the Presence of Discretization

## Abstract

Abstract Testing conditional independence (CI) has many important applications, such as Bayesian network learning and causal discovery. Although several approaches have been developed for learning CI structures for observed variables, those existing methods generally fail to work when the variables of interest can not be directly observed and only discretized values of those variables are available. For example, if X 1 subscript 𝑋 1 X_{1} italic_X start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , X ~ 2 subscript ~ 𝑋 2 \tilde{X}_{2} over~ start_ARG italic_X end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and X 3 subscript 𝑋 3 X_{3} italic_X start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT are the observed variables, where X ~ 2 subscript ~ 𝑋 2 \tilde{X}_{2} over~ start_ARG italic_X end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT is a discretization of the latent variable X 2 subscript 𝑋 2 X_{2} italic_X start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , applying the existing methods to the observations of X 1 subscript 𝑋 1 X_{1} italic_X start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , X ~ 2 subscript ~ 𝑋 2 \tilde{X}_{2} over~ start_ARG italic_X end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and X 3 subscript 𝑋 3 X_{3} italic_X start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT would lead to a false conclusion about the underlying CI of variables X 1 subscript 𝑋 1 X_{1} italic_X start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , X 2 subscript 𝑋 2 X_{2} italic_X start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT and X 3 subscript 𝑋 3 X_{3} italic_X start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT .
Motivated by this, we propose a CI test specifically designed to accommodate the presence of discretization. To achieve this, a bridge equation and nodewise regression are used to recover the precision coefficients reflecting the conditional dependence of the latent continuous variables under the nonparanormal model. An appropriate test statistic has been proposed, and its asymptotic distribution under the null hypothesis of CI has been derived.
Theoretical analysis, along with empirical validation on various datasets, rigorously demonstrates the effectiveness of our testing methods.

## 1 Introduction

Independence and conditional independence (CI) are fundamental concepts in statistics. They are leveraged for exploring queries in statistical inference, such as sufficiency, parameter identification, and ancillarity (Dawid, ). They also play a central role in emerging areas such as causal discovery (Koller & Friedman, ), graphical model learning, and feature selection (Xing et al., ). Tests for CI have attracted increasing attention from both theoretical and application sides.

Formally, the problem is to test the CI of two variables $X_{j_{1}}$ and $X_{j_{2}}$ given a random vector (a set of other variables) $\bm{Z}$. In statistical notation, the null hypothesis is written as $H_{0}:X_{j_{1}}\perp X_{j_{2}}\mid\bm{Z}$, where $\perp$ denotes “independent from.” The alternative hypothesis is written as $H_{1}:X_{j_{1}}\not\perp X_{j_{2}}\mid\bm{Z}$, where $\not\perp$ denotes “dependent with.” The null hypothesis implies that once $\bm{Z}$ is known, the values of $X_{j_{1}}$ provide no additional information about $X_{j_{2}}$, and vice versa. Different tests have been designed to handle different scenarios, including Gaussian variables with linear dependence (Yuan & Lin, ; Peterson et al., ; Mohan et al., ; Ren et al., ) and non-linear dependence (Fukumizu et al., ; Zhang et al., ; Strobl et al., ; Sen et al., ; Aliferis et al., ) (For detailed related work, please refer to App. D).

Given observations of $X_{j_{1}}$, $X_{j_{2}}$, and $\bm{Z}$, the CI can be effectively tested with existing methods.
However, in many scenarios, accurately measuring continuous variables of interest is challenging due to limitations in data collection. Sometimes the data obtained are approximations represented as discretized values.
For example, in finance, variables such as asset values cannot be measured and are binned into ranges for assessing investment risks (e.g., sell, hold, and strong buy) (Changsheng & Yongfeng, ; Damodaran, ). Similarly, in mental health, anxiety levels are often assessed using scales like the GAD-7, which categorizes responses into levels such as mild, moderate, or severe (Mossman et al., ; Johnson et al., ).
In the entertainment industry, the quality of movies is typically summarized through viewer ratings (Sparling & Sen, ; Dooms et al., ).

When discretization is present, existing CI tests can fail to determine the CI of underlying continuous variables. This issue arises because existing CI tests treat discretized observations as observations of continuous variables, leading to incorrect conclusions about their CI relationships. More precisely, the problem lies in the discretization process, which introduces new discrete variables.
Consequently, although the intent is to test the CI of the underlying continuous variables, what is actually being tested is the CI involving a mix of both continuous and newly introduced discrete variables. In general, this CI relationship is inconsistent with the one among the underlying continuous variables.

As illustrated in Fig. LABEL:fig:1, we show different data-generative processes using causal graphical models (Pearl, ) in the presence of discretization. A gray node indicates an observable variable, while a white node indicates a latent variable. Variables denoted by $X_{j}$ (without a tilde $\sim$) represent continuous variables, which may not be observed; while variables denoted by $\tilde{X}_{j}$ represent observed discretized variables derived from $X_{j}$ due to discretization.
In Fig. LABEL:fig:1a, $X_{2}$ is latent, and only its discrete counterpart $\tilde{X}_{2}$ is observed. In this case, rather than observing $X_{1}$, $X_{2}$, and $X_{3}$, we only observe $X_{1}$, $\tilde{X}_{2}$, and $X_{3}$. Existing CI methods use these observations to test whether $X_{1}\perp X_{3}\mid\{X_{2}\}$, but what is actually being tested is whether $X_{1}\perp X_{3}\mid\{\tilde{X}_{2}\}$.
In fact, according to the causal Markov condition (Spirtes et al., ), , it can be inferred from Fig. LABEL:fig:1a that $X_{1}\perp X_{3}\mid\{X_{2}\}$ and $X_{1}\not\perp X_{3}\mid\{\tilde{X}_{2}\}$. This mismatch leads to existing CI methods, that employ observations to check the CI relationships between $X_{1}$ and $X_{3}$ given $X_{2}$, to reach incorrect conclusions. Due to the same reason, checking the CI also fails in Fig LABEL:fig:1b and Fig LABEL:fig:1c.

In this paper, we design a CI test specifically for handling the presence of discretization. An appropriate test statistic for the CI of latent continuous variables, based solely on discretized observations, is derived. The key is to build connections between the discretized observations and the parameters needed for testing the CI of the latent continuous variables. To achieve this, we first develop bridge equations that allow us to estimate the covariance of the underlying continuous variables with discretized observations. Then, we leverage a node-wise regression (Callot et al., ) to derive appropriate test statistics for CI relationships from the estimated covariance. By assuming that the continuous variables follow a Gaussian distribution, we can derive the asymptotic distributions of the test statistics under the null hypothesis of CI. The major contributions of our paper include that

- •
We develop a CI test for ensuring accurate analysis in scenarios where data has been discretized, which are common due to limitations in data collection or measurement techniques, such as in financial analysis and healthcare.
- •
Our CI test can handle various scenarios including 1). Both variables $X_{j_{1}}$ and $X_{j_{2}}$ are discretized 2). Both variables $X_{j_{1}}$ and $X_{j_{2}}$ are continuous. 3). One of the variables $X_{j_{1}}$ or $X_{j_{2}}$ is discretized.
- •
We compare our test with the existing methods on both synthetic and real-world datasets, confirming that our method can effectively estimate the CI of the underlying continuous variables and outperform the existing tests applied on the discretized observations.

## 2 Problem Setting and necessity of correction

##### Problem Setting

Consider a set of independent and identically distributed (i.i.d.) $p$-dimensional random vectors, denoted as $\tilde{\bm{X}}=(X_{1},X_{2},\ldots,\tilde{X}_{j},\ldots,\tilde{X}_{p})^{T}$. In this set, some variables, indicated by a tilde ($\sim$), such as $\tilde{X}_{j}$, follow a discrete distribution. For each such variable, there exists a corresponding latent Gaussian random variable $X_{j}$. The transformation from $X_{j}$ to $\tilde{X}_{j}$ is governed by an unknown monotone nonlinear function $g_{j}$ and a thresholding function $f_{j}$. The function $f_{j}\circ g_{j}:\mathcal{X}\rightarrow\tilde{\mathcal{X}}$ maps the continuous domain of $X_{j}$ onto the discrete domain of $\tilde{\mathcal{X}}_{j}$.
Specifically, for each variable $X_{j}$, there exists a finite constant vector $\mathbf{d}_{j}=(d_{j,1},\dots,d_{j,M})$ characterized by strictly increasing elements such that

$$ $\tilde{X}_{j}=f_{j}(g_{j}(X_{j}))=\begin{cases}1&0<g_{j}(X_{j})<d_{j,1}\\ m&d_{j,m-1}<g_{j}(X_{j})<d_{j,m}\\ M&g_{j}(X_{j})>d_{j,m}\\ \end{cases}$ (1) $$

This model is also known as the nonparanormal model (Liu et al., ).
The cardinality of the domain after discretization is at least 2 and smaller than infinity. Our goal is to assess both conditional and unconditional independence among the variables of the vector $\bm{X}=(X_{1},X_{2},\ldots,X_{j},\ldots,X_{p})^{T}$. In our model, we assume $\bm{X}\sim N(0,\Sigma)$, $\Sigma$ only contain 1 among its diagonal, i.e., $\sigma_{jj}=1$ for all $j\in[1,\dots,p]$. One should note this assumption is without loss of generality. We provide a detailed discussion of our assumption in App. A.9.

##### Why the correction is essential?

We aim to propose a CI test that serves as a correction to infer the correct CI relationships among the latent continuous variables of interest. One question that arises is whether the discretized variables exhibit the same conditional independence as their original continuous counterparts, i.e., the correction is not needed. This concern becomes more significant when the level of discretization is high. To highlight the effect of discretization, we propose the following theorem. In essence, the discretization inevitably introduces distortions, leading to potentially false conclusions. The proof can be found in Appendix A.1.

###### Theorem 2.1 .

Let $X_{1},X_{2}$ and $X_{3}$ be jointly Gaussian random variables such that $X_{1}\perp\!\!\!\perp X_{3}|X_{2}$, $\tilde{X}_{2}=f_{j}(g_{j}(X_{2}))$ is the descretized observation as defined in equation 1. Then the conditional independence between $X_{1}$ and $X_{3}$ given $\tilde{X}_{2}$ doesn’t hold, i.e., $X_{1}\not\perp\!\!\!\perp X_{3}|\tilde{X}_{2}$.

###### Theorem 2.1 .

## 3 DCT: A CI Test in the Presence of Discretization

To develop an independence test, one needs to design a test statistic that can reflect the dependence relation and be calculated from observations. Next, it is essential to derive the underlying distribution of this statistic under the null hypothesis that the tested variables are conditionally (or unconditionally) independent. By calculating the value of the test statistic from observations and determining if this statistic is likely to be sampled from the derived distribution (i.e., calculating the p-value and comparing it with the significance level $\alpha$), we can decide if the null hypothesis should be rejected.

Our objective is to deduce the independence and CI relationships within the original multivariate Gaussian model, based on its discretized observations.
In the context of a multivariate Gaussian model, this challenge is directly equivalent to constructing statistical inferences for its covariance matrix $\mathbf{\Sigma}=(\sigma_{j_{1},j_{2}})$ and its precision matrix $\bm{\Omega}=(\omega_{j,k})=\mathbf{\Sigma}^{-1}$ (Baba et al., ). The covariance matrix $\mathbf{\Sigma}$ captures the pairwise covariances between variables, while the precision matrix $\bm{\Omega}$ provides information about the CI between variables. Specifically, the entry $\omega_{j,k}$ in the precision matrix is the partial correlation coefficient between variables $X_{j}$ and $X_{k}$, which can be used to test whether these variables are conditionally independent given some other variables.
Technically, we are interested in two things: (1) the calculation of the covariance $\hat{\sigma}_{j_{1},j_{2}}$ and the precision coefficient (or the partial correlation coefficient) $\hat{\omega}_{j,k}$, serving as the estimation of $\sigma_{j_{1},j_{2}}$ and $\omega_{j,k}$ respectively (in this paper, a variable with a hat indicates its estimation); and (2) the derivation of the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ and $\hat{\omega}_{j,k}-\omega_{j,k}$ under the null hypothesis of independence and CI.

In the remainder of section, 1). we first introduce bridge equations to address the estimation challenge of the covariance $\sigma_{j_{1},j_{2}}$; 2). we proceed to derive the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$, demonstrating it is asymptotically normal; 3). utilizing nodewise regression, we establish the relationship between the covariance matrix $\bm{\Sigma}$ and the precision matrix $\bm{\Omega}$, where the regression parameter $\beta_{j,k}$ acts as an effective surrogate for $\omega_{j,k}$. Leveraging the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$, we further illustrate that $\hat{\beta}_{j,k}-\beta_{j,k}$ is also asymptotically normal.

### 3.1 Design Bridge Equation for Test Statistics

##### Estimating Covariance with Bridge Equations

The bridge equation establishes a connection between the underlying covariance $\sigma_{j_{1},j_{2}}$ of two continuous variables $X_{j_{1}}$ and $X_{j_{2}}$ with the observations.
When in the presence of discretization, the discrete transformations make the sample covariance matrix based on $\tilde{\bm{X}}$ inconsistent with the covariance matrix of $\bm{X}$.
To obtain the estimation $\hat{\sigma}_{j_{1},j_{2}}$ of $\sigma_{j_{1},j_{2}}$, the bridge equation is leveraged.
In general, its form is as follows.

$$ $\hat{\tau}_{j_{1},j_{2}}=T(\sigma_{j_{1},j_{2}};\hat{\bm{\Lambda}}),$ (2) $$

where $\sigma_{j_{1},j_{2}}$ is the covariance needed to be estimated, $\hat{\tau}_{j_{1},j_{2}}$ is a statistic that can also be estimated from observations, and $\hat{\bm{\Lambda}}$ is a set of additional parameters required by the function $T(\cdot)$. The specific form of the function $T(\cdot)$ will be derived later. Both $\hat{\tau}_{j_{1},j_{2}}$ and $\hat{\bm{\Lambda}}$ should be able to be calculated purely relying on observations. Then, given the calculated $\hat{\tau}_{j_{1},j_{2}}$ and $\bm{\hat{\Lambda}}$, $\hat{\sigma}_{j_{1},j_{2}}$ can be obtained by solving the bridge equation $\hat{\tau}_{j_{1},j_{2}}=T(\sigma_{j_{1},j_{2}};\hat{\bm{\Lambda}})$.
As a result, the covariance matrix $\bm{\Sigma}$ of $\bm{X}$ can be estimated, which contains information about both unconditional independence and CI (which can be derived from its inverse).

To estimate the covariance of a latent multivariate Gaussian distribution, we need to design appropriate $\hat{\tau}_{j_{1},j_{2}}$, $\hat{\bm{\Lambda}}$, and $T(\cdot)$. Notably, bridge equations have to be designed to handle all three possible cases: C1. both observed variables are discretized; C2. one variable is continuous while the other is discretized; and C3. both variables remain continuous. We will show that cases C1 and C2 can be merged into a single form of bridge equation with different parameters and a binarization operation applied to the observations. Our bridge equations are presented in Def. 3.2, Def. 3.3, and Def. 3.4.

##### Bridge Equations for Discretized and Mixed Pairs

Let us first address the challenging cases where both observed variables are discretized or where one variable is continuous while the other is discretized. In general, different bridge equations would need to be designed to handle each case individually. However, in our analysis, we provide a unified bridge equation that is applicable to both cases. This is achieved by binarizing the observed variables, thereby unifying both cases into a binary case. As some information may be lost in the binarization process, this unification may require more examples compared to using tailored bridge functions for each specific case. Developing specific bridge equations for each case to improve sample efficiency is left in future work.

Intuitionally, for the original continuous variable $X_{j}$, binarization separates it into two parts based on a boundary $h_{j}$: the part for $X_{j}$ larger than $h_{j}$ and the part for $X_{j}$ smaller than $h_{j}$. In this case, we can estimate the boundary by calculating the proportion of $X_{j}$ that exceeds the boundary. In the scenario of two variables where the threshold $h_{j_{1}}$ and $h_{j_{2}}$ divide the space into four regions, the proportions of these areas are influenced by the covariance $\sigma_{j_{1},j_{2}}$, which connects the relation between the binarized variables with the latent covariance. This approach allows us to initially estimate the threshold $h_{j_{1}}$, $h_{j_{2}}$ of a pair of variables, followed by estimating the covariance $\sigma_{j_{1},j_{2}}$.

Let $\mathbb{P}_{n}Z$ denote the average of a random variable $Z$ given $n$ i.i.d. observation of $Z$ and $E[Z]$ as the true mean of $Z$, $\mathbb{P}$ as the probability and $\hat{P}$ as the empirical probability. We then define the boundary $h_{j}$ as follows: for any single discretized variable $\tilde{X}_{j}$, there exists a constant $c_{j}$ such that:

$$ $\mathbbm{1}\{\tilde{x}_{j}^{i}>E[\tilde{X}_{j}]\}=\mathbbm{1}\{g_{j}(x_{j}^{i} )>c_{j}\}=\mathbbm{1}\{x_{j}^{i}>h_{j}\},$ $$

where $h_{j}=g_{j}^{-1}(c_{j})$. Specifically, $h_{j}$ is the boundary in the original continuous domain to determine if the discretized observation $\tilde{X}_{k}$ is larger than its mean. When the continuous variable $X_{j}$ follows a normal distribution, there is a relation $\mathbb{P}(\tilde{X}_{j}>E[\tilde{X}_{j}])=1-\Phi(h_{j})$, where $\Phi$ is the cumulative distribution function (cdf) of a standard normal distribution.
We then provide the following definition:

###### Definition 3.1 .

The estimated boundary can be expressed as $\hat{h}_{j}=\Phi^{-1}(1-\hat{\tau}_{j})$, where $\hat{\tau}_{j}=\sum^{n}_{i=1}\mathbbm{1}_{\{\tilde{x}_{j}^{i}>\mathbb{P}_{n}
\tilde{X}_{j}\}}/n$, serving as the estimation of $\mathbb{P}(\tilde{X}_{j}>E[\tilde{X}_{j}])$.

Let $\bar{\Phi}(z_{1},z_{2};\rho)=\mathbb{P}(Z_{1}>z_{1},Z_{2}>z_{2})$, where $(Z_{1},Z_{2})^{T}$ follows a bivariate normal distribution with mean zero, variance one and covariance $\rho$. We define

$$ $\displaystyle\tau_{j_{1},j_{2}}$ $\displaystyle=\mathbb{P}(\tilde{x}_{j_{1}}^{i}>E[\tilde{X}_{j_{1}}],\tilde{x}_ {j_{2}}^{i}>E[\tilde{X}_{j_{2}}])=\bar{\Phi}(h_{j_{1}},h_{j_{2}};\sigma_{j_{1} ,j_{2}}).$ (3) $$

That is, the proportion of discretized variables larger than their mean can be expressed as a function of underlying covariance.
This equation serves as the key of estimating latent covariance based on the discretized observations. Specifically, we can substitute those true parameters with their estimation and construct the bridge equation to get the estimated covariance:

###### Definition 3.2 (Bridge Equation for A Discretized-Variable Pair) .

For discretized variables $\tilde{X}_{j_{1}}$ and $\tilde{X}_{j_{2}}$, the bridge equation is defined as:

$$ $\displaystyle\hat{\tau}_{j_{1},j_{2}}=\hat{P}(\tilde{X}_{j_{1}}>\mathbb{P}_{n} \tilde{X}_{j_{1}},\tilde{X}_{j_{2}}>\mathbb{P}_{n}\tilde{X}_{j_{2}})=\frac{1}{ n}\sum^{n}_{i=1}\mathbbm{1}_{\{\tilde{x}^{i}_{j_{1}}>\mathbb{P}_{n}\tilde{X}_{ j_{1}},\tilde{x}_{j_{2}}^{i}>\mathbb{P}_{n}\tilde{X}_{j_{2}}\}}=T(\sigma_{j_{1 },j_{2}};\{\hat{h}_{j_{1}},\hat{h}_{j_{2}}\}),$ $\displaystyle\text{and the function }T(\sigma_{j_{1},j_{2}};\{\hat{h}_{j_{1}}, \hat{h}_{j_{2}}\}):=\bar{\Phi}(\hat{h}_{j_{1}},\hat{h}_{j_{2}};\sigma)=\int_{x _{1}>\hat{h}_{j_{1}}}\int_{x_{2}>\hat{h}_{j_{2}}}\phi(x_{j_{1}},x_{j_{2}}; \sigma)d{x_{j_{1}}}d{x_{j_{2}}},$ $$

where $\phi$ is the probability density function of a bivariate normal distribution,
$\hat{h}_{j_{1}},\hat{h}_{j_{2}}$ can be simply calculated using Def. 3.1.

Following the same intuition, we can directly apply the same bridge equation to estimate the covariance of mixed pairs. The only difference is there is no need to estimate the boundary $\hat{h}_{j}$ for the continuous variable. Instead, we can incorporate its true mean of zero into the equation.

###### Definition 3.3 (Bridge Equation for A Continuous-Discretized-Variable Pair) .

For one continuous variable $X_{j_{1}}$ and one discretized variable $\tilde{X}_{j_{2}}$, the bridge function is defined as follows:

$$ $\displaystyle\hat{\tau}_{j_{1},j_{2}}=\hat{P}(X_{j_{1}}>0,\tilde{X}_{j_{2}}> \mathbb{P}_{n}\tilde{X}_{j_{2}})=\frac{1}{n}\sum^{n}_{i=1}\mathbbm{1}_{\{x^{i} _{j_{1}}>0,\tilde{x}_{j_{2}}^{i}>\mathbb{P}_{n}\tilde{X}_{j_{2}}\}}=T(\sigma_{ j_{1},j_{2}};\{0,\hat{h}_{j_{2}}\}),$ $$

and the function $T(\cdot)$ has the same form of Def. 3.2.

###### Definition 3.1 .

###### Definition 3.2 (Bridge Equation for A Discretized-Variable Pair) .

###### Definition 3.3 (Bridge Equation for A Continuous-Discretized-Variable Pair) .

##### A Bridge Equation for A Continuous-Variable Pair

When there is no discretized transformation, the sample covariance of $X_{j_{1}}$ and $X_{j_{2}}$ provides a consistent estimation.
In this context, the function $T$ acts merely as an identity mapping.

###### Definition 3.4 (A Bridge Equation for A Continuous-Variable Pair) .

For two continuous variables $X_{j_{1}}$ and $X_{j_{2}}$ , the bridge equation is defined as:

$$ $\displaystyle\hat{\tau}_{j_{1},j_{2}}:=\hat{\sigma}_{j_{1},j_{2}}=\frac{1}{n} \sum_{i=1}^{n}x_{j_{1}}^{i}x_{j_{2}}^{i}-\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i }\frac{1}{n}\sum_{i=1}^{n}x_{j_{2}}^{i}=T(\sigma_{j_{1},j_{2}};\emptyset).$ $$

For two continuous variables $X_{j_{1}}$ and $X_{j_{2}}$, the analytic solution of the estimated covariance can be simply obtained using Def. 3.4.

###### Definition 3.4 (A Bridge Equation for A Continuous-Variable Pair) .

##### Calculation of Estimated Covariance

For the continuous case,
the analytic solution of $\hat{\sigma}_{j_{1},j_{2}}$ can be simply obtained using Def. 3.4.
For the cases involving the discretized variable as proposed in Def. 3.2 and Def. 3.3,
we can rely on the property that variance $\bm{\Sigma}$ only contains $1$ among the diagonal, which implies the covariance $\sigma_{j_{1},j_{2}}$ should vary from $-1$ to $1$. Thus, we can calculate the estimated covariance by solving the objective

$$ $\min_{\sigma_{j_{1},j_{2}}}||\hat{\tau}_{j_{1},j_{2}}-T(\sigma_{j_{1},j_{2}}; \{\hat{h}_{j_{1}},\hat{h}_{{j}_{2}}\})||^{2}\quad s.t.-1<\sigma_{j_{1},j_{2}}<1.$ (4) $$

The $\hat{\tau}_{j_{1},j_{2}}$ is a one-to-one mapping with calculated $\hat{\sigma}_{j_{1},j_{2}}$, $\hat{h}_{j_{1}}$ and $\hat{h}_{j_{2}}$, which is proved in App. A.3

### 3.2 Unconditional Independence Test

The estimation of covariance $\hat{\sigma}_{j_{1},j_{2}}$ can be effectively solved using the designed bridge equation. Now, we focus on deriving the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$. These results is used as an unconditional independence test in the presence of the discretization. Moreover, Thm. 3.5, Lem. 3.6, Lem. 3.7 and Lem. 3.8 will be leveraged in the derivation process of the CI test in Section 3.3.
The detailed derivation steps for both unconditional independence test and CI test are relatively complicated, therefore, we will provide a general intuition. For a complete derivation, please refer to the App. A.4.

Assume we are interested in the true parameter $\theta_{0}$. We denote $\hat{\theta}$ as its estimation which is close to $\theta_{0}$, and $f(\theta)$ is a continuous function. By leveraging Taylor expansion, we have

$$ $f(\hat{\theta})=f(\theta_{0})+f^{\prime}(\theta_{0})(\hat{\theta}-\theta_{0}),$ (5) $$

which directly constructs the relationship between the estimated parameter with the true one. Rearrange the term, we get $\hat{\theta}-\theta_{0}=(f(\hat{\theta})-f(\theta_{0}))/f^{\prime}(\theta_{0})$. If the denominator is a constant and the numerator can be expressed as a sum of i.i.d samples, we can see $\hat{\theta}-\theta_{0}$ will be asymptotically normal according to the central limit theorem (Van der Vaart, ).

Let $\psi_{\hat{\theta}}=[f^{1}_{\hat{\theta}}(\cdot),f^{2}_{\hat{\theta}}(\cdot),f
^{3}_{\hat{\theta}}(\cdot)]^{T}$ contains a group of functions parameterized by $\hat{{\theta}}$ (For discretized pairs, $\hat{\theta}=(\hat{\sigma}_{j_{1},j_{2}},\hat{h}_{j_{1}},\hat{h}_{j_{2}})$). Define $\mathbb{P}_{n}\psi_{\hat{\theta}}$ as sample mean of these functions evaluated at $n$ sample points. Similarly, $\mathbb{P}_{n}\psi_{\hat{\theta}}{\psi_{\hat{\theta}}}^{T}$ is defined as sample mean of the outer product ${\psi_{\hat{\theta}}}{\psi_{\hat{\theta}}}^{T}$. The notation $P{\psi_{\hat{\theta}}}:=E\mathbb{P}_{n}{\psi_{\hat{\theta}}}$ denotes the expectations of the functions in ${\psi_{\hat{\theta}}}$. Furthermore, let ${\psi_{\hat{\theta}}}^{\prime}$ denote the derivative of the functions contained in ${\psi_{\hat{\theta}}}$.
We now provide the main result of derived distribution $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ under the hull hypothesis that test pairs are independent.

###### Theorem 3.5 (Independence Test) .

In our settings, under the null hypothesis that two observed variables indexed with ${j_{1}}$ and ${j_{2}}$ are statistically independent under our framework, i.e., $\sigma_{j_{1},j_{2}}=0$, the independence can be tested using the statistic

$$ $\hat{\sigma}_{j_{1},j_{2}}=T^{-1}(\hat{\tau}_{j_{1},j_{2}};\hat{\theta}).$ $$

This statistic is approximated to follow a normal distribution, as detailed below:

$$ $\hat{\sigma}_{j_{1},j_{2}}\overset{\text{approx}}{\sim}N\left(0,\frac{1}{n}(( \mathbb{P}_{n}\psi^{\prime}_{\hat{\theta}})^{-1}\mathbb{P}_{n}\psi_{\hat{ \theta}}\psi_{\hat{\theta}}^{T}(\mathbb{P}_{n}{\psi_{\hat{\theta}}^{\prime T}} )^{-1})_{1,1}\right),$ (6) $$

where the specific form of $\psi_{\hat{\theta}}$ are presented in Lem. 3.6,Lem. 3.7 and Lem. 3.8.

We now provide the specific forms of $\psi_{\hat{\theta}}$. Since the variables being tested for independence can be both discretized, only one being discretized, or neither being discretized. This results in different forms of $\psi_{\hat{\theta}}$ consequently differs across these scenarios. Let $Z_{j_{1}}$ and $Z_{j_{2}}$ be any two random variables indexed by ${j_{1}}$ and ${j_{2}}$. Let $\hat{\sigma}_{j_{1},j_{2}}^{i}={z}^{i}_{j_{1}}\cdot{z}^{i}_{j_{2}}-\mathbb{P}_
{n}{Z}_{j_{1}}\cdot\mathbb{P}_{n}{Z}_{j_{2}}$ denote the sample covariance based on a $i$-th pairwise observation of the variables $Z_{j_{1}}$ and $Z_{j_{2}}$. Let $\hat{\tau}^{i}_{j_{1}}=\mathbbm{1}_{\{z^{i}_{j_{1}}>\mathbb{P}_{n}Z_{j_{1}}\}}$ and $\hat{\tau}^{i}_{j_{2}}=\mathbbm{1}_{\{Z^{i}_{j_{2}}>\mathbb{P}_{n}Z_{j_{2}}\}}$, each calculated based on $i$-th observations of the variables $Z_{j_{1}}$ and $Z_{j_{2}}$, respectively. Let $\hat{\tau}^{i}_{j_{1},j_{2}}$ be $\hat{\tau}^{i}_{j_{1}}\cdot\hat{\tau}^{i}_{j_{2}}$. We further denote $\bar{\Phi}(\cdot)=1-\Phi(\cdot)$.
The different forms of $\psi_{\hat{\theta}}$ that arise in different cases are defined as follows:

###### Lemma 3.6 .

($\psi_{\hat{\theta}}$ for A Continuous-Variable Pair).
For two continuous variables $X_{j_{1}}$ and $X_{j_{2}}$,

$$ $\displaystyle\psi_{\hat{\theta}}$ $\displaystyle:=\hat{\sigma}_{j_{1},j_{2}}^{i}-\hat{\sigma}_{j_{1},j_{2}}.$ (7) $$

###### Lemma 3.7 ( $\psi_{\hat{\theta}}$ for A Discretized-Variable Pair) .

For discretized variables $\tilde{X}_{j_{1}}$ and $\tilde{X}_{j_{2}}$,

$$ $\displaystyle\psi_{\hat{\theta}}:=\begin{pmatrix}\hat{\tau}^{i}_{j_{1},j_{2}}- T(\hat{\sigma}_{j_{1},j_{2}};\{\hat{h}_{j_{1}},\hat{h}_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{1}}-\bar{\Phi}(\hat{h}_{j_{1}})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(\hat{h}_{j_{2}})\end{pmatrix}.$ (8) $$

###### Lemma 3.8 ( $\psi_{\hat{\theta}}$ for A Continuous-Discretized-Variable Pair) .

For one discretized variable $\tilde{X}_{j_{2}}$ and one continuous variable $X_{j_{1}}$,

$$ $\displaystyle\psi_{\hat{\theta}}$ $\displaystyle:=\begin{pmatrix}\hat{\tau}^{i}_{j_{1},j_{2}}-T(\hat{\sigma}_{j_{ 1},j_{2}};\{0,\hat{h}_{j_{2}})\}\\ \hat{\tau}^{i}_{j_{1}}-\bar{\Phi}(\hat{h}_{j_{2}})\\ \end{pmatrix}.$ (9) $$

Derivation of forms of $\psi_{\hat{\theta}}$ for different cases and their corresponding distribution defined in Eq equation 6 can be found in App. A.5, App. A.6, App. A.7.
Up to this point, our discussion has been confined to the case of covariance $\sigma_{j_{1},j_{2}}$, the indicator of unconditional independence.
In the next section, we will present the results of our CI test.

###### Theorem 3.5 (Independence Test) .

###### Lemma 3.6 .

###### Lemma 3.7 ( ψ θ ^ subscript 𝜓 ^ 𝜃 \psi_{\hat{\theta}} italic_ψ start_POSTSUBSCRIPT over^ start_ARG italic_θ end_ARG end_POSTSUBSCRIPT for A Discretized-Variable Pair) .

###### Lemma 3.8 ( ψ θ ^ subscript 𝜓 ^ 𝜃 \psi_{\hat{\theta}} italic_ψ start_POSTSUBSCRIPT over^ start_ARG italic_θ end_ARG end_POSTSUBSCRIPT for A Continuous-Discretized-Variable Pair) .

### 3.3 Conditional Independence (CI) Test

To construct a CI test of our model, we are interested at two things:
calculation of the estimated precision coefficient $\hat{\omega}_{j,k}$ and the derivation of the corresponding distribution $\hat{\omega}_{j,k}-\omega_{j,k}$.
In the following, we first build $\beta_{j,k}$, which is obtained using nodewise regression and show it serves as a surrogate of testing for $\omega_{j,k}=0$, we then construct the formulation of $\hat{\beta}_{j,k}-\beta_{j,k}$ as the combination of formulation of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ and show it will also be asymptotically normal.

##### Nodewise Regression for CI

To utilize covariance for testing CI, it is necessary to establish a relationship between the estimated covariance and a metric capable of reflecting CI.
To achieve this, we employ the nodewise regression which effectively builds the connection between covariance and precision matrix.
Suppose we can access observations $\{\bm{x}^{1},\bm{x}^{2},\dots,\bm{x}^{n}\}$
from latent continuous variables $\bm{X}=(X_{1},\dots,X_{p})\sim N(0,\bm{\Sigma})$, nodewise regression will do regression on every dimension with all other dimensions as predictors.

$$ $x_{j_{1}}^{i}=\sum_{j_{1}\neq j_{2}}x_{j_{2}}^{i}\beta_{j}+\epsilon_{j_{1}}^{i}.$ (10) $$

It can be shown that there are deterministic relationships between the regression coefficients and the covariance and precision matrices of $\bm{X}$, as illustrated below and proved in App. A.8.1.

$$ $\displaystyle\beta_{j}=\bm{\Sigma}_{-j-j}^{-1}\bm{\Sigma}_{-jj}\in\mathbb{R}^{ p-1},\quad\beta_{j,k}=-\frac{\omega_{j,k}}{\omega_{j,j}},\quad j\not=k,$ (11) $$

where $\bm{\Sigma}_{-j-j}$ is the submatrix of $\bm{\Sigma}$ without $j$th column and $j$th row, and the $\bm{\Sigma}_{-jj}$ is the vector of $j$th column without $j$th row. $\beta_{j,k}\in\mathbb{R}$ is the surrogate of $\omega_{j,k}$ to capture the independence relationship of $X_{j}$ with $X_{k}$ conditioning on other variables. We can use Def.  3.2, Def. 3.3 and Def. 3.4 to get the estimation $\hat{\bm{\Sigma}}_{-j-j}$ and $\hat{\bm{\Sigma}}_{-jj}$ and thus get the estimation $\hat{\beta}_{j}$.

##### Statistical Inference for β j , k subscript 𝛽 𝑗 𝑘 \beta_{j,k} italic_β start_POSTSUBSCRIPT italic_j , italic_k end_POSTSUBSCRIPT

Nodewise regression offers a robust solution for the estimation problem. A pertinent inquiry pertains to the construction of the distribution of $\hat{\beta}_{j}-\beta_{j}$. It is crucial to recognize that the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ is already established. Therefore, if we can conceptualize $\hat{\beta}_{j}-\beta_{j}$ as a linear combination of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$, the problem is directly solved, i.e., the $\hat{\beta}_{j}-\beta_{j}$ is linear combination of dependent Gaussian variables. The underlying relationship between these variables is as follows:

$$ $\hat{\beta}_{j}-\beta_{j}=-\hat{\bm{\Sigma}}_{-j-j}^{-1}\left((\hat{\bm{\Sigma }}_{-j-j}-\bm{\Sigma}_{-j-j})\beta_{j}-(\hat{\bm{\Sigma}}_{-jj}-\bm{\Sigma}_{- jj})\right).$ $$

The derivation is provided in App. A.8.2. For ease of notation, we further express the distribution of the difference between the estimated covariance and the true covariance as

$$ $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n}\xi^{ i}_{j_{1},j_{2}}.$ (12) $$

The specific form of $\xi_{j_{1},j_{2}}^{i}$ is given in App. A.5, A.6, A.7 respectively for different cases.
For notational convenience, we express $\hat{\bm{\Sigma}}_{-j-j}-\bm{\Sigma}_{-j-j}=\frac{1}{n}\sum_{i=1}^{n}\Xi^{i}_{
-j,-j}$ and $\hat{\bm{\Sigma}}_{-jj}-\bm{\Sigma}_{-jj}=\frac{1}{n}\sum_{i=1}^{n}\Xi^{i}_{-j
,j}$, where $\xi_{j_{1},j_{2}}$ is the element of the matrix $\Xi$ at the position indexed by $(j_{1},j_{2})$.
We now propose the statistic and its asymptotic distribution for the CI test in the following theorem.

###### Theorem 3.9 (Conditional Independence test) .

In our settings, under the null hypothesis that $X_{j}$ and $X_{k}$ are conditional statistically independent given a set of variables $\bm{Z}$, i.e., $\beta_{j,k}=0$, the statistic

$$ $\hat{\beta}_{j,k}=(\hat{\bm{\Sigma}}_{-j-j}^{-1}\hat{\bm{\Sigma}}_{-jj})_{[k]},$ (13) $$

where $[k]$ denotes the element corresponding to the variable $X_{k}$ in $\hat{\bm{\Sigma}}_{-j-j}^{-1}\hat{\bm{\Sigma}}_{-jj}$.
The statistic $\hat{\beta}_{j,k}$ has the asymptotic distribution:

$$ $\hat{\beta}_{j,k}\sim N(0,{a^{[k]}}^{T}\frac{1}{n^{2}}\sum_{i=1}^{n}vec(B^{i}_ {-j})vec(B^{i}_{-j})^{T})a^{[k]}),$ $$

$$ $\text{ where }B^{i}=\begin{bmatrix}{\Xi_{-j,j}^{i}}^{T}\\ \Xi_{-j,-j}^{i}\end{bmatrix},\quad a^{[k]}_{l}=\begin{cases}\left(\hat{\bm{ \Sigma}}_{-j-j}^{-1}\right)_{[k],l},&\text{for }l\in\{1,\dots,p-1\}\\ \sum_{q=1}^{n}\left(\hat{\bm{\Sigma}}_{-j-j}^{-1}\right)_{[k],l}\left(\tilde{ \beta}_{j}\right)_{q},&\text{for }l\in\{p,\dots,p^{2}-p\}\end{cases}$ $$

and $\tilde{\beta}_{j}$ is $\beta_{j}$ whose $\beta_{j,k}=0$; $vec$ is row-wise vectorization of a matrix.

In practice, we can plug in the estimation of regression parameter $\hat{\beta}_{j}$ and set $\hat{\beta}_{j,k}=0$ as the substitution of $\tilde{\beta}_{j}$ to calculate the variance and do the CI test.
Specifically, we can obtain the $\hat{\beta}_{j,k}$ using equation 13 where the estimated covariance terms can be calculated by solving the bridge equation Eq. 2. Under the null hypothesis that $\beta_{j,k}=0$ (conditional independence), we can take the calculated $\hat{\beta}_{j,k}$ into the distribution defined in Thm. 3.9 and obtain the p-value. If the p-value is smaller than the predefined significance level $\alpha$ (normally set at 0.05), we will infer the tested pairs are conditionally dependent; otherwise, we do not.
The detailed derivation of the Thm. 3.9 can be found in App. A.8.2.

###### Theorem 3.9 (Conditional Independence test) .

### 3.4 Discussion of assumption

Discretization of continuous variables inevitably leads to information loss in the original data. Compared to the original distributional information, the recovered covariance matrix is naturally less accurate. Given this, constructing a valid statistical inference procedure, rather than solely relying on estimated covariance values for drawing conditional independence conclusions, is desirable.

One major limitation of DCT is its reliance on the assumption that latent continuous variables follow a multivariate normal distribution. Violations of this assumption can lead to erroneous conclusions. For instance, consider a scenario where the relationship between latent variables is nonlinear, such as $X_{j_{1}}=X_{j_{2}}^{2}$. In this case, the covariance $\sigma_{j_{1},j_{2}}$ equals zero despite a deterministic dependency between $X_{j_{1}}$ and $X_{j_{2}}$. Consequently, even if the correlation is perfectly estimated, the model fails to capture the true underlying relationship, leading to incorrect inferences.

Nevertheless, although the theoretical framework of DCT requires latent continuous variables to follow a multivariate Gaussian distribution,
experimental results in various settings, even in situations in which this assumption is violated, demonstrate the usefulness and robustness of DCT, suggesting the development of this technique is essential to causal discovery from discretized continuous data.
Further details of the empirical validations are provided in Appendix C.

Figure: Figure 3: Comparison of results of Type I and calibrated Type II error ($1-\textup{power}$) for all three types of tested data (continuous, mixed, discrete) and different number of samples and cardinality of conditioning set. The suffix attached to a test’s name denotes the cardinality of discretization; for example, "Fisherz_4" signifies the application of the Fisher-z test to data discretized into four levels. Chi-square test is only applicable for the discrete case.
Refer to caption: x1.png

## 4 Experiments

We applied the proposed method DCT to synthetic data to evaluate its practical performance and compare it with Fisher-Z test (Fisher, ) (for all three data types) and Chi-Square test (F.R.S., ) (for discrete data only) as baselines. Specifically, we investigated its Type I and Type II error and its application in causal discovery. The experiments investigating its robustness, performance in denser graphs and effectiveness in a real-world dataset can be found in App. C.

### 4.1 On the Effect of the Cardinality of Conditioning Set and the Sample Size

Our experiment investigates the variations in Type I and Type II error (1 minus power) probabilities under two conditions. In the first scenario, we focus on the effects of modifying the sample size, denoted as $n=(100,500,1000,2000)$, while conditioning on a single variable. In the second, the sample size is held constant at 2000, and we vary the cardinality of the conditioning set, represented as $D=(1,2,\ldots,5)$. It is assumed that every variable within this conditioning set is effective, i.e., they influence the CI of the tested pairs. We repeat each test $1500$ times.

We use $Y,W$ to denote the variables being tested and use $Z$ to denote the variables being conditioned on. The discretized versions of the variables are denoted with a tilde symbol (e.g., $\tilde{Z}$).
For both conditions, we evaluate three distinct types of observations of tested variables: continuous observations for both variables ($Y,W$), discrete observations for both variables ($\tilde{Y},\tilde{W}$) and a mixed type ($\tilde{Y},W$). The variables in the conditioning set will always be discretized observations ($\tilde{Z}$).

To see how well the derived asymptotic null distribution approximates the true one, we verify if the probability of Type I error aligns with the significance level $\alpha$ preset in advance. We generate true continuous multivariate Gaussian data $Y,W$ from $Z_{i}$ (single $i=1$ for the first scenario, and summed over $n$ for the second), structured as $a_{i}Z_{i}+E$ and $\sum_{i=1}^{n}a_{i}Z_{i}+E$, where $a_{i}$ is sampled from $U(0.5,1.5)$ and $E$ follows a standard normal distribution, independent of all other variables. This ensures $Y\perp\!\!\!\perp W|Z$. The data are then discretized into $K=(2,4,8,12)$ levels, with boundaries randomly set based on the variable range. The first column in Fig. 3 (a) (b) shows the resulting probability of Type I errors at the significance level $\alpha=0.05$ compared with other methods.

A good test should have as small a probability of Type II error as possible, i.e., a larger power. To test the power of our DCT, we generate the continuous multivariate Gaussian data $Z_{i}$ from $Y,W$; constructed as $Z_{i}=a_{i}Y+b_{i}W+E$, where $a_{i},b_{i}$ are sampled from $U(0.5,1.5)$ and $E$ follows a standard normal distribution independent with all others, i.e., $Y\not\perp\!\!\!\perp W|Z$. The same discretization approach is applied here. One should note that directly comparing the p-value with a common predefined significance level is unfair since all baselines tend to produce very small p-values. Therefore, all tests are calibrated(^1^11Calibration is the process of empirically finding the decision threshold to match the desired significance level, ensuring accurate control of Type I error.) in this experiment.
The second column in Fig. 3 (a) and (b) correspondingly shows the calibrated Type II error as the number of samples and the cardinality of the conditioning set change, compared to other methods.

From Fig. 3 (a), we note that the Type I error rates with our derived null distribution are well-approximated at 0.05 across all three data types in both scenarios. In contrast, other testing methods show significantly higher Type I error rates, increasing with the number of samples and the size of the conditioning set. This indicates that such methods are more prone to erroneously concluding that tested variables are conditionally dependent.
Additionally, while alternative tests demonstrate considerable power with smaller sample sizes, our approach requires a sample size of 1000 to achieve satisfactory power, particularly in mixed and continuous cases. A possible explanation for this phenomenon is that our method binarizes discretized data, which may not effectively utilize all observations. This aspect warrants further investigation in future research. Moreover, our test shows remarkable stability in response to changes in the number of conditioning sets.

Figure: Figure 5: Experimental result of skeleton discovery on synthetic data for changing sample size (a) and changing number of nodes (b). Fisherz_nodis is the Fisher-z test applied to original continuous data. We evaluate $F_{1}$ ($\uparrow$), Precision ($\uparrow$), Recall ($\uparrow$) and SHD ($\downarrow$).

### 4.2 Application in Causal Discovery

Causal discovery aims at looking for the true causal structure from the data. Under the assumption of causal Markov condition that the causal relationships among variables can be expressed by a Directed Acyclic Graph (DAG) $\mathcal{G}$ and its statistical independence is entailed in this graphic model, faithfulness ensures that the statistical independencies observed in the data can be reliably used to infer the causal structure. Given both assumptions, constraint-based causal discovery, e.g., PC algorithm (Spirtes et al., ) recovers the graph structure relying on testing the conditional independence of observation. Apparently, in the presence of discretization, the failures of testing conditional independence will seriously impair the resulting DAG.

To evaluate the efficacy of the DCT, we construct the true DAG $\mathcal{G}$ utilizing the Bipartite Pairing (BP) model as detailed in (Asratian et al., ), with the number of edges being one fewer than the number of nodes. The subsequent generation of true multivariate Gaussian data involves assigning causal weights drawn from a uniform distribution $U\sim(0.5,2)$ and incorporating noise via samples from a standard normal distribution for each variable.

Following this, we binarize the data, setting the threshold randomly based on each variable’s range. Our experiment is divided into two scenarios: In the first, we set the number of samples $n=5000$, with the number of nodes $p$ varying across 4, 6, 8, and 10. In the second scenario, we fix the number of nodes at $p=8$ and explore sample sizes $n=(500,1000,5000,10000)$.

A comparative analysis is performed using the PC algorithm integrated with various testing methods. Specifically, we compare DCT against the Fisher-z test applied to discretized data, the Chi-Square test, and the Fisher-z test on the original continuous data, the latter serving as a theoretical upper bound. Since the PC algorithm only returns a completed partially directed acyclic graph (CPDAG), we apply the same orientation rules from Dor & Tarsi (), as implemented by Causal-DAG (Chandler Squires, ), to convert a CPDAG into a DAG for easier comparison. We evaluate both the undirected skeleton and the directed graph using structural Hamming distance (SHD), F1 score, precision, and recall as evaluation metrics. For each setting, we run 10 graph instances with different seeds and report the mean and standard deviation for skeleton discovery in Figure 5 and DAG discovery in Figure LABEL:fig:_cd_dag in Appendix B.

According to the result, DCT exhibits performance nearly on par with the theoretical upper bound across metrics such as F1 score, precision, and Structural Hamming Distance (SHD) when the number of variables ($p$) is small and the sample size ($n$) is large. Despite a decline in performance as the number of variables increases with a smaller sample size, DCT significantly outperforms both the Fisher-Z test and the Chi-square test. Notably, in almost all settings, the recall of DCT is lower than that of the baseline tests, which is a reasonable outcome since these tests tend to infer conditional dependencies, thereby retaining all edges given the discretized observations. For instance, a fully connected graph, would achieve a recall of 1.

## 5 Conclusion

In this paper, we present a new testing method tailored for scenarios commonly encountered in real-world applications, where variables, though inherently continuous, are only observable in their discretized forms. Our method distinguishes itself from existing CI tests by effectively mitigating the misjudgment introduced by discretization and accurately recovering the CI relationships of latent continuous variables. We substantiate our approach with theoretical results and empirical validation, underscoring the effectiveness of our testing methods.

## Appendix A Proof and Derivations

### A.1 Proof of Thm. 2.1

If the $X_{1},X_{2}$ and $X_{3}$ are jointly Gaussian and $X_{1}\perp\!\!\!\perp X_{3}|X_{3}$, we have

$$ $Cov(X_{1},X_{3}|X_{2})=0.$ (14) $$

To test if $X_{1},X_{3}$ are conditional independent given $\tilde{X}_{2}$, we are interested if $Cov(X_{1},X_{3}|\tilde{X}_{2})$ equals zero. Using the law of total covariance, we have

$$ $Cov(X_{1},X_{3}|\tilde{X}_{2})=E[Cov(X_{1},X_{3}|X_{2},\tilde{X}_{2})|\tilde{X }_{2}]+Cov(E[X_{1}|X_{2},\tilde{X}_{2}],E[X_{3}|X_{2},\tilde{X}_{2}]|\tilde{X} _{2}).$ (15) $$

Since $\tilde{X}_{2}$ is the deterministic function of $X_{2}$, $\tilde{X}_{2}$ will be conditional independent with $X_{1}$ and $X_{3}$ given $X_{2}$. Therefore,

$$ $Cov(X_{1},X_{3}|X_{2},\tilde{X}_{2})=Cov(X_{1},X_{3}|X_{2})=0.$ (16) $$

The first term of equation 15 is zero. We now focus on the second term. Similarly, we have

$$ $E[X_{1}|X_{2},\tilde{X}_{2}]=E[X_{1}|X_{2}],\quad E[X_{3}|X_{2},\tilde{X}_{2}] =E[X_{3}|X_{2}],$ (17) $$

due to the conditional independence.
One can see

$$ $Cov(X_{1},X_{3}|X_{2},\tilde{X}_{2})=Cov(E[X_{1}|X_{2}],E[X_{3}|X_{2}]|\tilde{ X}_{2}).$ (18) $$

Without loss of generality, we assume mean of $X_{1},X_{2}$ and $X_{3}$ are zero. Then $E[X_{1}|X_{2}]$ and $E[X_{3}|X_{2}]$ are scaled versions of $X_{2}$. The original equation becomes

$$ $Cov(X_{1},X_{3}|X_{2},\tilde{X}_{2})=c\cdot Var(X_{2}|\tilde{X}_{2}),$ (19) $$

where $c$ is a constant. $Var(X_{2}|\tilde{X}_{2})$ will never be zero since $\tilde{X}_{2}$ is the discretized version of $X_{2}$ (information loss). Thus, $X_{1}\not\perp\!\!\!\perp X_{3}|\tilde{X}_{2}$, which completes the proof.

### A.2 Proof of θ ^ ⁢ → 𝑝 ⁢ θ 0 ^ 𝜃 𝑝 → subscript 𝜃 0 \hat{\theta}\overset{p}{\to}\theta_{0} over^ start_ARG italic_θ end_ARG overitalic_p start_ARG → end_ARG italic_θ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT

###### Lemma A.1 .

For the estimation $\hat{\theta}$ which is calculated using bridge equation 3.4 3.2 3.3, as a zero of $\Psi_{n}$ defined in equation 32,equation 39, equation 42 , will converge in probability to $\theta_{0}=(\sigma_{j_{1},j_{2}},h_{j_{1}},h_{j_{2}}),(\sigma_{j_{1},j_{2}},h_
{j_{2}}),(\sigma_{j_{1},j_{2}})$ respectively.

Proof We first focus on the most challenging one where both variables are discrete. According to the law of large numbers, for the estimated boundary $\hat{h}_{j_{1}}$ and $\hat{h}_{j_{2}}$ whose calculations are defined as $\hat{h}_{j}=\Phi^{-1}(1-\hat{\tau}_{j})$, we should have

$$ $n\rightarrow\infty,\quad\hat{\tau}_{j}=\frac{1}{n}\sum_{i=1}^{n}\mathbbm{1}_{ \{\tilde{x}_{j}^{i}>\mathbb{P}_{n}\tilde{X}_{j}\}}\overset{p}{\to}\mathbb{P}( \tilde{X}_{j}>E[\tilde{X}_{j}]).$ (20) $$

Recall the definition $\mathbb{P}(\tilde{X}_{j}>E[\tilde{X}_{j}])=1-\Phi(h_{j})$, according to continuous mapping theorem (Vaart, ), as long as the function $\Phi^{-1}(1-\cdot)$ is continuous, we should have
$\hat{h}_{j}\overset{p}{\to}h_{j}$. And thus $\hat{h}_{j_{1}}\overset{p}{\to}h_{j_{1}}$, $\hat{h}_{j_{2}}\overset{p}{\to}h_{j_{2}}$.

We have $\hat{\tau}_{j_{1},j_{2}}=\bar{\Phi}(\hat{h}_{j_{1}},\hat{h}_{j_{2}},\hat{
\sigma}_{j_{1},j_{2}})$ and the estimation $\hat{\sigma}_{j_{1},j_{2}}$ can be obtained through solving the function. Similarly, we also have

$$ $n\rightarrow\infty,\quad\hat{\tau}_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n} \mathbbm{1}_{\{\tilde{x}_{j_{1}}^{i}>\mathbb{P}_{n}\tilde{X}_{j_{1}}\}} \mathbbm{1}_{\{\tilde{x}_{j_{2}}^{i}>\mathbb{P}_{n}\tilde{X}_{j_{2}}\}} \overset{p}{\to}\mathbb{P}(\tilde{x}_{j_{1}}^{i}>E[\tilde{X}_{j_{1}}],\tilde{x }_{j_{2}}^{i}>E[\tilde{X}_{j_{2}}])=\tau_{j_{1},j_{2}}.$ (21) $$

Similarly, according to the continuous mapping theorem, we have $\hat{\sigma}_{j_{1},j_{2}}\overset{p}{\to}\sigma_{j_{1},j_{2}}$. Thus, the parameter $(\hat{\sigma}_{j_{1},j_{2}},\hat{h}_{j_{1}},\hat{h}_{j_{2}})\overset{p}{\to}(
\sigma_{j_{1},j_{2}},h_{j_{1}},h_{j_{2}})$.

Apparently, the result above could easily extend to the mixed case where we fix $\hat{h}_{1}=h_{1}=0$. Using the same procedure, we should have $(\hat{\sigma}_{j_{1},j_{2}},\hat{h}_{j_{2}})\overset{p}{\to}(\sigma_{j_{1},j_{
2}},h_{j_{2}})$.

For the continuous case whose estimated variance is calculated as $\hat{\sigma}_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i}x_{j_{2}}^{i}
-\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i}\frac{1}{n}\sum_{i=1}^{n}x_{j_{2}}^{i}.$, according to law of large numbers, we should have

$$ $n\rightarrow\infty,\quad\hat{\sigma}_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n}x_ {j_{1}}^{i}x_{j_{2}}^{i}-\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i}\frac{1}{n}\sum _{i=1}^{n}x_{j_{2}}^{i}\overset{p}{\to}E(X_{j_{1}}X_{j_{2}})-E(X_{j_{1}})E(X_{ j_{2}})=\sigma_{j_{1},j_{2}}.$ (22) $$

###### Lemma A.1 .

### A.3 Proof of one-to-one mapping between τ ^ j 1 , j 2 subscript ^ 𝜏 subscript 𝑗 1 subscript 𝑗 2 \hat{\tau}_{j_{1},j_{2}} over^ start_ARG italic_τ end_ARG start_POSTSUBSCRIPT italic_j start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_j start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT with σ ^ j 1 , j 2 subscript ^ 𝜎 subscript 𝑗 1 subscript 𝑗 2 \hat{\sigma}_{j_{1},j_{2}} over^ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT italic_j start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_j start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT

###### Lemma A.2 .

For any fixed $\hat{h}_{j_{1}}$ and $\hat{h}_{j_{2}}$,
$T(\sigma_{j_{1},j_{2}};\{\hat{h}_{j_{1}},\hat{h}_{j_{2}}\})=\int_{x_{1}>\hat{h
}_{j_{1}}}\int_{x_{2}>\hat{h}_{j_{2}}}\phi(x_{j_{1}},x_{j_{2}};\sigma)d{x_{j_{
1}}}d{x_{j_{2}}}$,
is a strictly monotonically increasing function on $\sigma\in(-1,1)$.

Proof To prove the lemma, we just need to show the gradient $\frac{\partial T(\sigma_{j_{1},j_{2}};\{\hat{h}_{j_{1}},\hat{h}_{j_{2}}\}}{
\partial\sigma}>0$ for $\sigma\in(-1,1)$.

$$ $\frac{\partial T(\sigma_{j_{1},j_{2}});\{\hat{h}_{j_{1}},\hat{h}_{j_{2}}\}}{ \partial\sigma}==\frac{1}{2\pi\sqrt{(1-\sigma^{2})}}{\exp\left(-\frac{(\hat{h} _{j_{1}}^{2}-2\sigma\hat{h}_{j_{1}}\hat{h}_{j_{2}}+\hat{h}_{j_{2}}^{2})}{2(1- \sigma^{2})}\right)},$ (23) $$

which is obviously positive for $\sigma\in(-1,1)$. Thus, we have one-to-one mapping between $\hat{\tau}_{j_{1}j_{2}}$ with the calculated $\hat{\sigma}_{j_{1},j_{2}}$ for fixed $\hat{h}_{j_{1}}$ and $\hat{h}_{j_{2}}$.

###### Lemma A.2 .

### A.4 Proof of Thm. 3.5

In this section, we provide the proof of Thm. 3.5, which utilizes a regular statistical tool: Z-estimator (Vaart, ). Specifically,
we are interested in the parameter $\theta$ and we have it estimation $\hat{\theta}$.
Let $\bm{x_{1}},\dots,\bm{x_{n}}$ are sampled from some true distribution $P$, we can construct the function characterized by the parameter $\theta$ related the $\bm{x}$ as $\psi_{\theta}(\bm{x})$. As long as we have $n$ observations, we can construct the function as follows

$$ $\Psi_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\psi_{\theta}(\bm{x}_{i})=\mathbb{P} _{n}\psi_{\theta}.$ (24) $$

We further specify the form

$$ $\Psi(\theta)=\int\psi_{\theta}(\bm{x})d\bm{x}=P\psi_{\theta}.$ (25) $$

Assume the estimator $\hat{\theta}$ is a zero of $\Psi_{n}$, i.e., $\Psi_{n}(\hat{\theta})=0$ and will converge in probability to $\theta_{0}$, which is a zero of $\Psi$, i.e., $\Psi(\theta_{0})=0$. Expand $\Psi_{n}(\hat{\theta})$ in a Taylor series around $\theta_{0}$, we should have

$$ $0=\Psi_{n}(\hat{\theta})=\Psi_{n}(\theta_{0})+(\hat{\theta}-\theta_{0})\Psi_{n }^{\prime}(\theta_{0})+\frac{1}{2}(\hat{\theta}-\theta_{0})^{2}\Psi_{n}^{ \prime\prime}(\theta_{0}).$ (26) $$

Rearrange the equation above, we have

$$ $\displaystyle\hat{\theta}-\theta_{0}$ $\displaystyle=-\frac{\Psi_{n}(\theta_{0})}{\Psi_{n}^{\prime}(\theta_{0})+\frac {1}{2}(\hat{\theta}-\theta_{0})^{2}\Psi_{n}^{\prime\prime}(\theta_{0})}$ (27) $\displaystyle=-\frac{\frac{1}{n}\sum_{i=1}^{n}\psi_{\theta}(\bm{x}_{i})}{\Psi_ {n}^{\prime}(\theta_{0})+\frac{1}{2}(\hat{\theta}-\theta_{0})^{2}\Psi_{n}^{ \prime\prime}(\theta_{0})}.$ $$

According to the central limit theorem, the numerator will be asymptotic normal with variance $P\psi_{\theta_{0}}^{2}/n$ as the mean $\Psi(\theta_{0})=0$ is zero. The first term of denominator $\Psi_{n}^{\prime}(\theta_{0})$ will converge in probability to $\Psi^{\prime}(\theta_{0})$ according to the law of large numbers. The second term $\hat{\theta}-\theta_{0}=o_{P}(1)$. (^2^22We will not provide proof of this in this paper; however, interested readers may refer to (Vaart, )) As long as the denominator converges in probability and the numerator converges in distribution, according to Slusky’s lemma, we have

$$ $\sqrt{n}(\hat{\theta}-\theta_{0})\rightsquigarrow N\left(0,\frac{P\psi_{\theta _{0}}^{2}}{(P\psi_{\theta_{0}}^{\prime})^{2}}\right).$ (28) $$

Extend into the high-dimensional case we should have

$$ $\displaystyle\hat{\theta}-\theta_{0}$ $\displaystyle=-(\Psi_{n}^{\prime}(\theta_{0}))^{-1}\Psi_{n}(\theta_{0}),$ (29) $$

where the second order term is omitted, further assume the matrix $P\psi^{\prime}_{\theta_{0}}$ is invertible, we have

$$ $\sqrt{n}(\hat{\theta}-\theta_{0})\rightsquigarrow N\left(0,(P\psi^{\prime}_{ \theta_{0}})^{-1}P\psi_{\theta_{0}}\psi_{\theta_{0}}^{T}(P\psi^{\prime T}_{ \theta_{0}})^{-1}\right),$ (30) $$

Specifically, in our case $\theta_{0}$ = ($\sigma_{j_{1},j_{2}},\bm{\Lambda}$), where $\bm{\Lambda}$ is another parameter set influencing the estimation of $\sigma_{j_{1},j_{2}}$ (will discuss case in case in later proof). In the practical scenario, we only have access to the estimated parameter $\hat{\theta}$ and the empirical distribution $\mathbb{P}_{n}$, thus we have

$$ $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}\overset{\text{approx}}{\sim}N \left(0,((\mathbb{P}_{n}\psi^{\prime}_{\hat{\theta}})^{-1}\mathbb{P}_{n}\psi_{ \hat{\theta}}\psi_{\hat{\theta}}^{T}(\mathbb{P}_{n}{\psi_{\hat{\theta}}^{ \prime T}})^{-1})_{1,1}\right).$ (31) $$

Under the null hypothesis of independent, $\sigma_{j_{1},j_{2}=0}$.
We provide the proof that $\hat{\theta}\overset{p}{\to}\theta_{0}$ of our case in App. A.2. Thus, $\mathbb{P}_{n}\psi_{\hat{\theta}}$, the function parameterized by $\hat{\theta}$, should also converge in $\mathbb{P}_{n}\psi_{\hat{\theta}_{0}}$ when $n\overset{}{\to}\infty$. Besides, by the law of large numbers, $\mathbb{P}_{n}\psi_{\hat{\theta}_{0}}$ will converge to $P\psi_{\hat{\theta}_{0}}$. Thus, the equation above will converge to equation 30 when $n\rightarrow\infty$.

### A.5 Derivation of Lem. 3.7

Let’s first focus on the most challenging case where both variables are discretized observations and our interested parameter will include $\hat{\theta}=(\hat{\sigma}_{j_{1},j_{2}},\hat{h}_{j_{1}},\hat{h}_{j_{2}})$ (Although we only care about the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$, the estimation of boundary $\hat{h}_{j_{1}}$and $\hat{h}_{j_{2}}$ will influence the estimation of $\hat{\sigma}_{j_{1},j_{2}}$, thus we need to consider all of them).

The next step will be to construct an appropriate criterion function $\psi$ such that $\Psi_{n}(\hat{\theta})=\mathbf{0}$. Given $n$ observations $\{\tilde{\bm{x}}^{1},\tilde{\bm{x}}^{2},\dots,\tilde{\bm{x}}^{n}\}$, which are discretized version of $\{\bm{x}^{1},\bm{x}^{2},\dots,\bm{x}^{n}\}$ we should have

$$ $\Psi_{n}(\hat{\theta})=\begin{pmatrix}\Psi_{n}(\hat{\sigma}_{j_{1},j_{2}})\\ \Psi_{n}(\hat{h}_{j_{1}})\\ \Psi_{n}(\hat{h}_{j_{2}})\end{pmatrix}=\frac{1}{n}\sum_{i=1}^{n}\psi_{\hat{ \theta}}(\tilde{\bm{x}}^{i})=\frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}\hat{\tau }_{j_{1},j_{2}}^{i}-T(\hat{\sigma}_{j_{1},j_{2}};\{\hat{h}_{j_{1}},\hat{h}_{j_ {2}}\})\\ \hat{\tau}^{i}_{j_{1}}-\bar{\Phi}(\hat{h}_{j_{1}})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(\hat{h}_{j_{2}})\end{pmatrix}=\mathbf{0}.\\$ (32) $$

$$ $\Psi_{n}(\theta_{0})=\begin{pmatrix}\Psi_{n}(\sigma_{j_{1},j_{2}})\\ \Psi_{n}(h_{j_{1}})\\ \Psi_{n}(h_{j_{2}})\\ \end{pmatrix}=\frac{1}{n}\sum_{i=1}^{n}\psi_{\theta_{0}}(\tilde{\bm{x}}^{i})= \frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}\hat{\tau}_{j_{1},j_{2}}^{i}-T(\sigma_ {j_{1},j_{2}};\{h_{j_{1}},h_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{1}}-\bar{\Phi}(h_{j_{1}})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(h_{j_{2}})\end{pmatrix}.$ (33) $$

The difference between the estimated parameter with the true parameter can be expressed as

$$ $\displaystyle\hat{\theta}-\theta_{0}=\begin{pmatrix}\hat{\sigma}_{j_{1},j_{2}} -\sigma_{j_{1},j_{2}}\\ \hat{h}_{j_{1}}-h_{j_{1}}\\ \hat{h}_{j_{2}}-h_{j_{2}}\end{pmatrix}=-\frac{1}{n}\sum_{i=1}^{n}$ $\displaystyle\begin{pmatrix}\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{ \partial\sigma_{j_{1},j_{2}}}&\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{ \partial h_{j_{1}}}&\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial h_{ j_{2}}}\\ \frac{\partial\Psi_{n}(h_{j_{1}})}{\partial\sigma_{j_{1},j_{2}}}&\frac{ \partial\Psi_{n}(h_{j_{1}})}{\partial h_{j_{1}}}&\frac{\partial\Psi_{n}(h_{j_{ 1}})}{\partial h_{j_{2}}}\\ \frac{\partial\Psi_{n}(h_{j_{2}})}{\partial\sigma_{j_{1},j_{2}}}&\frac{ \partial\Psi_{n}(h_{j_{2}})}{\partial h_{j_{1}}}&\frac{\partial\Psi_{n}(h_{j_{ 2}})}{\partial h_{j_{2}}}\end{pmatrix}^{-1}$ $\displaystyle\cdot\ \ \begin{pmatrix}\hat{\tau}_{j_{1},j_{2}}^{i}-T(\sigma_{j_ {1},j_{2}};\{h_{j_{1}},h_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{1}}-\bar{\Phi}(h_{j_{1}})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(h_{j_{2}})\end{pmatrix},$ (34) $$

where the specific form of each entry of the gradient matrix is expressed as

$$ $\displaystyle\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial\sigma_{j_{ 1},j_{2}}}$ $\displaystyle=-\frac{1}{2\pi\sqrt{(1-{\sigma}_{j_{1},j_{2}}^{2})}}{\exp\left(- \frac{(h_{j_{1}}^{2}-2{\sigma}_{j_{1},j_{2}}h_{j_{1}}h_{j_{2}}+h_{j_{2}}^{2})} {2(1-{\sigma}_{j_{1},j_{2}}^{2})}\right)};$ (35) $\displaystyle\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial h_{j_{1}}}$ $\displaystyle=\int_{h_{j_{2}}}^{\infty}\frac{1}{2\pi\sqrt{1-\sigma_{j_{1},j_{2 }}^{2}}}\exp\left(-\frac{h_{j_{1}}^{2}-2\sigma_{j_{1},j_{2}}h_{j_{1}}x_{2}+x_{ 2}^{2}}{2(1-\sigma_{j_{1},j_{2}}^{2})}\right)dx_{2};$ $\displaystyle\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial h_{j_{2}}}$ $\displaystyle=\int_{h_{j_{1}}}^{\infty}\frac{1}{2\pi\sqrt{1-\sigma_{j_{1},j_{2 }}^{2}}}\exp\left(-\frac{h_{2}^{2}-2\sigma_{j_{1},j_{2}}h_{j_{2}}x_{1}+x_{1}^{ 2}}{2(1-\sigma_{j_{1},j_{2}}^{2})}\right)dx_{1};$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{1}})}{\partial\sigma_{j_{1},j_{2}}}$ $\displaystyle=0;$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{1}})}{\partial h_{j_{1}}}$ $\displaystyle=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{h_{j_{1}}^{2}}{2}\right);$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{1}})}{\partial h_{j_{2}}}$ $\displaystyle=0;$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{2}})}{\partial\sigma_{j_{1},j_{2}}}$ $\displaystyle=0;$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{2}})}{\partial h_{j_{1}}}$ $\displaystyle=0;$ $\displaystyle\frac{\partial\Psi_{n}(h_{j_{2}})}{\partial h_{j_{2}}}$ $\displaystyle=\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{h_{j_{2}}^{2}}{2}\right).$ $$

For simplicity of notation, we define

$$ $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n}\xi^{ i}_{j_{1},j_{2}},$ (36) $$

where the specific form is of $\{\xi^{i}_{j_{1},j_{2}}\}$ is defined in equation A.5. We should note that $\{\xi^{i}_{j_{1},j_{2}}\}$ are i.i.d random variables with mean zero (this property will be the key to the derivation of inference of CI). As long as our estimation $\hat{\theta}$ converge in probability to $\theta_{0}$ as proved in A.2, we have

$$ $\sqrt{n}(\hat{\theta}-\theta_{0})\rightsquigarrow N\left(0,((P\psi^{\prime}_{ \theta_{0}})^{-1}P\psi_{\theta_{0}}\psi_{\theta_{0}}^{T}(P\psi^{\prime T}_{ \theta_{0}})^{-1})\right),$ (37) $$

where $\psi_{{\theta}_{0}}$ is defined in equation 33. However, in practice, we don’t have access to either $P$ or $\theta_{0}$. In this scenario, we can plug in the empirical distribution of $\mathbb{P}_{n}\psi_{\hat{\theta}}$ to get the estimated variance, i.e., the actual variance used in the calculation of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ is

$$ $\frac{1}{n}\left((\mathbb{P}_{n}\psi^{\prime}_{\hat{\theta}})^{-1}\mathbb{P}_{ n}\psi_{\hat{\theta}}\psi_{\hat{\theta}}^{T}(\mathbb{P}_{n}{\psi_{\hat{\theta} }^{\prime T}})^{-1}\right)_{1,1}.$ (38) $$

### A.6 Derivation of Lem. 3.8

Use the same line of procedure as in the derivation of Lem. 3.7, for mixed pair of observations where $X_{j_{1}}$ is continuous and $\tilde{X}_{j_{2}}$ is discrete, we can construct the criterion function

$$ $\Psi_{n}(\hat{\theta})=\begin{pmatrix}\Psi_{n}(\hat{\sigma}_{j_{1},j_{2}})\\ \Psi_{n}(\hat{h}_{j_{2}})\end{pmatrix}=\frac{1}{n}\sum_{i=1}^{n}\psi_{\hat{ \theta}}(\tilde{\bm{x}}^{i})=\frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}\hat{\tau }_{j_{1},j_{2}}^{i}-T(\hat{\sigma}_{j_{1},j_{2}};\{0,\hat{h}_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(\hat{h}_{j_{2}})\end{pmatrix}=\mathbf{0}.\\$ (39) $$

$$ $\Psi_{n}(\theta_{0})=\begin{pmatrix}\Psi_{n}(\sigma_{j_{1},j_{2}})\\ \Psi_{n}(h_{j_{2}})\\ \end{pmatrix}=\frac{1}{n}\sum_{i=1}^{n}\psi_{\theta_{0}}(\bm{\tilde{x}}^{i})= \frac{1}{n}\sum_{i=1}^{n}\begin{pmatrix}\hat{\tau}_{j_{1},j_{2}}^{i}-T(\sigma_ {j_{1},j_{2}};\{0,h_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(h_{j_{2}})\end{pmatrix}.$ (40) $$

The difference between the estimated parameter with the true parameter can be expressed as

$$ $\hat{\theta}-\theta_{0}=\begin{pmatrix}\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1 },j_{2}}\\ \hat{h}_{j_{2}}-h_{j_{2}}\end{pmatrix}=-\frac{1}{n}\sum_{i=1}^{n}\begin{ pmatrix}\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial\sigma_{j_{1},j_ {2}}}&\frac{\partial\Psi_{n}(\sigma_{j_{1},j_{2}})}{\partial h_{j_{2}}}\\ \frac{\partial\Psi_{n}(h_{j_{2}})}{\partial\sigma_{j_{1},j_{2}}}&\frac{ \partial\Psi_{n}(h_{j_{2}})}{\partial h_{j_{2}}}\end{pmatrix}^{-1}\begin{ pmatrix}\hat{\tau}_{j_{1},j_{2}}^{i}-T(\sigma_{j_{1},j_{2}};\{0,h_{j_{2}}\})\\ \hat{\tau}^{i}_{j_{2}}-\bar{\Phi}(h_{j_{2}}).\end{pmatrix},$ (41) $$

where the specific form of each entry of the gradient matrix can be found in equation 35. Using exactly the same procedure, we should have the same formation of the variance calculated as equation 38 with a different definition of $\psi_{\theta_{0}}$ and $\psi_{\hat{\theta}}$ defined in equation 40 equation 39.

### A.7 Derivation of Lem. 3.6

Use the same line of procedure as in derivation of Lem. 3.7, for a continuous pair of variables, we can construct the criterion function

$$ $\Psi_{n}(\hat{\theta})=\Psi_{n}(\hat{\sigma}_{j_{1},j_{2}})=\frac{1}{n}\sum_{i =1}^{n}x_{j_{1}}^{i}x_{j_{2}}^{i}-\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i}\frac{ 1}{n}\sum_{i=1}^{n}x_{j_{2}}^{i}-\hat{\sigma}_{j_{1},j_{2}}=0.$ (42) $$

$$ $\Psi_{n}(\theta_{0})=\Psi_{n}(\sigma_{j_{1},j_{2}})=\frac{1}{n}\sum_{i=1}^{n}x _{j_{1}}^{i}x_{j_{2}}^{i}-\frac{1}{n}\sum_{i=1}^{n}x_{j_{1}}^{i}\frac{1}{n} \sum_{i=1}^{n}x_{j_{2}}^{i}-\sigma_{j_{1},j_{2}}.$ (43) $$

Denote $\frac{1}{n}\sum_{i=1}^{n}x^{i}_{j_{1}}$ as $\bar{x}_{j_{1}}$ and $\frac{1}{n}\sum_{i=1}^{n}x^{i}_{j_{2}}$ as $\bar{x}_{j_{2}}$. We should have

$$ $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}=\frac{1}{n}\sum_{i=1}^{n}x_{j_ {1}}^{i}x_{j_{2}}^{i}-\bar{x}_{j_{1}}\bar{x}_{j_{2}}-\sigma_{j_{1},j_{2}}.$ (44) $$

According to equation 28, we have

$$ $\sqrt{n}(\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}})\rightsquigarrow N \left(0,\frac{P\psi_{\theta_{0}}^{2}}{(P\psi_{\theta_{0}}^{\prime})^{2}}\right).$ (45) $$

where $(P\psi_{\theta_{0}}^{\prime})^{2}=1$. In practical calculation, we have the variance

$$ $\frac{1}{n}\mathbb{P}_{n}\psi_{\hat{\theta}}^{2}/(\mathbb{P}_{n}\psi_{\hat{ \theta}}^{\prime})^{2}=\frac{1}{n^{2}}\sum_{i=1}^{n}(x_{j_{1}}^{i}x_{j_{2}}^{i }-\bar{x}_{j_{1}}\bar{x}_{j_{2}}-\hat{\sigma}_{j_{1},j_{2}})^{2}.$ (46) $$

### A.8 Proof of Thm. 3.9

#### A.8.1 Proof of Relation between Σ Σ \Sigma roman_Σ , Ω Ω \Omega roman_Ω with β 𝛽 \beta italic_β

Consider our latent continuous variables $\bm{X}=(X_{1},\dots,X_{p})\sim N(0,\bm{\Sigma})$ and do nodewise regression

$$ $X_{j}=X_{-j}\beta_{j}+\epsilon_{j}.$ (47) $$

We can divide its covariance $\bm{\Sigma}$ and its precision matrix $\Omega=\bm{\Sigma}^{-1}$ into $X$ and $Y$ part in our regression:

$$ $\begin{array}[]{cc}\bm{\Sigma}=\begin{pmatrix}\bm{\Sigma}_{jj}&\bm{\Sigma}_{j- j}\\ \bm{\Sigma}_{-jj}&\bm{\Sigma}_{-j-j}\end{pmatrix}&\bm{\Omega}=\begin{pmatrix} \bm{\Omega}_{jj}&\bm{\Omega}_{j-j}\\ \bm{\Omega}_{-jj}&\bm{\Omega}_{-j-j}\end{pmatrix}\\ &\end{array}.$ (48) $$

Just like regular linear regression, we can get

$$ $n\rightarrow\infty,\quad\beta_{j}=\bm{\Sigma}_{-j-j}^{-1}\bm{\Sigma}_{-jj}.$ (49) $$

From the invertibility of a block matrix

$$ $\begin{bmatrix}A&B\\ C&D\end{bmatrix}^{-1}=\left[\begin{array}[]{cc}(A-BD^{-1}C)^{-1}&-(A-BD^{-1}C) ^{-1}BD^{-1}\\ -D^{-1}C(A-BD^{-1}C)^{-1}&D^{-1}+D^{-1}C(A-BD^{-1}C)^{-1}BD^{-1}\end{array} \right].$ (50) $$

If $A$ and $D$ is invertible, we will have

$$ $\begin{bmatrix}A&B\\ C&D\end{bmatrix}^{-1}=\begin{bmatrix}(A-BD^{-1}C)^{-1}&0\\ 0&(D-CA^{-1}B)^{-1}\end{bmatrix}\begin{bmatrix}I&-BD^{-1}\\ -CA^{-1}&I\end{bmatrix}.$ (51) $$

Thus, we can get:

$$ $\displaystyle\bm{\Omega}_{jj}=(\bm{\Sigma}_{jj}-\bm{\Sigma}_{j-j}\bm{\Sigma}_{ -j-j}^{-1}\bm{\Sigma}_{-jj})^{-1};$ (52) $\displaystyle\bm{\Omega}_{j-j}=-\left(\bm{\Sigma}_{jj}-\bm{\Sigma}_{j-j}\bm{ \Sigma}_{-j-j}^{-1}\bm{\Sigma}_{-jj}\right)^{-1}\bm{\Sigma}_{j-j}(\bm{\Sigma}_ {-j-j})^{-1}.$ $$

Move one step forward:

$$ $-\bm{\Omega}_{jj}^{-1}\bm{\Omega}_{j-j}=\bm{\Sigma}_{j-j}(\bm{\Sigma}_{-j-j})^ {-1}.$ (53) $$

Take transpose for both sides, as long as $\bm{\Omega}$ is a symmetric matrix and $\bm{\Omega}_{-jj}=\bm{\Omega}_{j-j}^{T}$, we will have

$$ $-\bm{\Omega}_{jj}^{-1}\bm{\Omega}_{-jj}=\bm{\Sigma}_{-j-j}^{-1}\bm{\Sigma}_{- jj}=\beta_{j}.$ (54) $$

We should note testing $\bm{\Omega}_{-jj}=0$ is equivalent to testing $\beta_{j}=0$ as the $\bm{\Omega}_{jj}$ will always be nonzero. The variable $\bm{\Omega}_{-jj}$ captures the CI of $X_{j}$ with other variables. As long as the variable $\bm{\Omega}_{jj}$ is just one scalar, we can get

$$ $\beta_{j,k}=-\frac{\omega_{j,k}}{\omega_{j,j}}$ (55) $$

capturing the CI relationship between variable $X_{j}$ with $X_{k}$ conditioning on all other variables.

#### A.8.2 Detailed derivation of inference for β j subscript 𝛽 𝑗 \beta_{j} italic_β start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT

Nodewise regression allows us to use the regression parameter $\beta_{j}$ as the surrogate of $\Omega_{-jj}$. The problem now transfers to constructing the inference for $\beta_{j}$, specifically, the derivation of distribution of $\hat{\beta}_{j}-\beta_{j}$. The overarching concept is that we are already aware of the distribution of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ and we know that there exists a deterministic relationship between $\beta_{j}$ with $\bm{\Sigma}$. Consequently, we can express $\hat{\beta}_{j}-\beta_{j}$ as a composite of $\hat{\sigma}_{j_{1},j_{2}}-\sigma_{j_{1},j_{2}}$ to establish such an inference. Specifically, we have

$$ $\displaystyle\hat{\beta}_{j}-\beta_{j}$ $\displaystyle=\hat{\bm{\Sigma}}^{-1}_{-j-j}\hat{\bm{\Sigma}}_{-jj}-\bm{\Sigma} ^{-1}_{-j-j}\bm{\Sigma}_{-jj}$ (56) $\displaystyle=\hat{\bm{\Sigma}}^{-1}_{-j-j}\left(\hat{\bm{\Sigma}}_{-jj}-\hat{ \bm{\Sigma}}_{-j-j}\bm{\Sigma}^{-1}_{-j-j}\bm{\Sigma}_{-jj}\right)$ $\displaystyle=-\hat{\bm{\Sigma}}^{-1}_{-j-j}\left(\hat{\bm{\Sigma}}_{-j-j} \beta_{j}-\bm{\Sigma}_{-j-j}\beta_{j}+\bm{\Sigma}_{-j-j}\beta_{j}-\hat{\bm{ \Sigma}}_{-jj}\right)$ $\displaystyle=-\hat{\bm{\Sigma}}_{-j-j}^{-1}\left((\hat{\bm{\Sigma}}_{-j-j}- \bm{\Sigma}_{-j-j})\beta_{j}-(\hat{\bm{\Sigma}}_{-jj}-\bm{\Sigma}_{-jj})\right),$ $$

where each entry in matrix $(\hat{\bm{\Sigma}}_{-j-j}-\bm{\Sigma}_{-j-j})$ and $(\hat{\bm{\Sigma}}_{-jj}-\bm{\Sigma}_{-jj})$ denotes
the difference between estimated covariance with true covariance.

Suppose that we want to test the CI of the variable $X_{1}$ with other variables, $j=1$. then

$$ $\displaystyle\hat{\bm{\Sigma}}_{-1-1}-\bm{\Sigma}_{-1-1}$ $\displaystyle=\begin{bmatrix}\hat{\sigma}_{2,2}\dots\hat{\sigma}_{2,p}\\ \dots\\ \hat{\sigma}_{p,2}\dots\hat{\sigma}_{p,p}\end{bmatrix}-\begin{bmatrix}{\sigma} _{2,2}\dots\sigma_{2,p}\\ \dots\\ {\sigma}_{p,2}\dots\sigma_{p,p}\end{bmatrix}$ (57) $\displaystyle:=\frac{1}{n}\sum_{i=1}^{n}\begin{bmatrix}{\xi}^{i}_{2,2}\dots\xi ^{i}_{2,p}\\ \dots\\ {\xi}^{i}_{p,2}\dots\xi^{i}_{p,p}\end{bmatrix},$ (58) $$

where $\{\xi^{i}_{j_{1},j_{2}}\}$ are i.i.d random variables with specific form defined in equation A.5 for discrete case, equation 41 for mixed case and equation 44 in continuous case.
Put them together:

$$ $\displaystyle\hat{\beta}_{1}-\beta_{1}=\begin{bmatrix}\hat{\beta}_{1,2}-\beta_ {1,2}\\ \hat{\beta}_{1,3}-\beta_{1,3}\\ \ldots\\ \hat{\beta}_{1,p}-\beta_{1,p}\end{bmatrix}$ $\displaystyle=-\hat{\bm{\Sigma}}_{-1-1}^{-1}\frac{1}{n}\sum_{i=1}^{n}\left( \begin{bmatrix}\xi^{i}_{2,2}&\xi^{i}_{2,3}&\ldots&\xi^{i}_{2,p}\\ \xi^{i}_{3,2}&\xi^{i}_{3,3}&\ldots&\xi^{i}_{3,p}\\ \ldots&\ldots&\ldots&\ldots\\ \xi^{i}_{p,2}&\xi^{i}_{p,3}&\dots&\xi^{i}_{p,p}\end{bmatrix}\begin{bmatrix} \beta_{1,2}\\ \beta_{1,3}\\ \dots\\ \beta_{1,p}\end{bmatrix}-\begin{bmatrix}\xi^{i}_{2,1}\\ \xi^{i}_{3,1}\\ \dots\\ \xi^{i}_{p,1}\\ \end{bmatrix}\right).$ (59) $$

As $\frac{1}{n}\sum_{i=1}^{n}\xi^{i}_{j_{1},j_{2}}$ is asymptotically normal, the who vector of $\hat{\beta}_{1}-\beta_{1}$ is a linear combination of Gaussian distribution. However, We cannot merely engage in a linear combination of its variance as they are dependent with each other. For example,
if $Y_{1},Y_{2}$ are dependent and we are trying to find out $Var(aY_{1}+bY_{2})$, we should have

$$ $Var(aY_{1}+bY_{2})=\begin{bmatrix}a&b\end{bmatrix}\begin{bmatrix}Var(Y_{1})& Cov(Y_{1},Y_{2})\\ Cov(Y_{1},Y_{2})&Var(Y_{2})\end{bmatrix}\begin{bmatrix}a\\ b\end{bmatrix}.$ (60) $$

Now, suppose we are interested in the distribution of $\hat{\beta}_{1,2}-\beta_{1,2}$, we have

$$ $\hat{\beta}_{1,2}-\beta_{1,2}=\frac{1}{n}\sum_{i=1}^{n}(\hat{\bm{\Sigma}}_{-1- 1}^{-1})_{[2],:}\left(\begin{bmatrix}\xi^{i}_{2,2}&\xi^{i}_{2,3}&\ldots&\xi^{i }_{2,p}\\ \xi^{i}_{3,2}&\xi^{i}_{3,3}&\ldots&\xi^{i}_{3,p}\\ \ldots&\ldots&\ldots&\ldots\\ \xi^{i}_{p,2}&\xi^{i}_{p,3}&\dots&\xi^{i}_{p,p}\end{bmatrix}\begin{bmatrix} \beta_{1,2}\\ \beta_{1,3}\\ \dots\\ \beta_{1,p}\end{bmatrix}-\begin{bmatrix}\xi^{i}_{2,1}\\ \xi^{i}_{3,1}\\ \dots\\ \xi^{i}_{p,1}\\ \end{bmatrix}\right),$ (61) $$

where $(\hat{\bm{\Sigma}}_{-1-1}^{-1})_{[2],:}$ is the row of index of $X_{2}$ of $\hat{\bm{\Sigma}}_{-1-1}^{-1}$ ($[2]$ denotes the index of the variable). For ease of notation, let

$$ $\Xi^{i}_{-1,-1}=\begin{bmatrix}\xi^{i}_{2,2}&\xi^{i}_{2,3}&\ldots&\xi^{i}_{2,p }\\ \xi^{i}_{3,2}&\xi^{i}_{3,3}&\ldots&\xi^{i}_{3,p}\\ \ldots&\ldots&\ldots&\ldots\\ \xi^{i}_{p,2}&\xi^{i}_{p,3}&\dots&\xi^{i}_{p,p}\end{bmatrix},\quad\quad\Xi^{i} _{-1,1}=\begin{bmatrix}\xi^{i}_{2,1}\\ \xi^{i}_{3,1}\\ \dots\\ \xi^{i}_{p,1}\\ \end{bmatrix},$ (62) $$

and let

$$ $B^{i}_{-1}=\begin{pmatrix}\xi^{i}_{2,1}&\xi^{i}_{3,1}&\dots&\xi^{i}_{p,1}\\ \xi^{i}_{2,2}&\xi^{i}_{2,3}&\ldots&\xi^{i}_{2,p}\\ \xi^{i}_{3,2}&\xi^{i}_{3,3}&\ldots&\xi^{i}_{3,p}\\ \ldots&\ldots&\ldots&\ldots\\ \xi^{i}_{p,2}&\xi^{i}_{p,3}&\dots&\xi^{i}_{p,p}\end{pmatrix}=\begin{bmatrix}{ \Xi_{-j,j}^{i}}^{T}\\ \Xi_{-j,-j}^{i}\end{bmatrix}$ (63) $$

as the concatenation of those two matrices. Similarly as equation 60,
The variance is calculated as

$$ $Var\left(\sqrt{n}(\hat{\beta}_{1,2}-\beta_{1,2})\right)={a^{[2]}}^{T}\frac{1}{ n}\sum_{i=1}^{n}vec(B^{i}_{-1})vec(B^{i}_{-1})^{T}a^{[2]},$ (64) $$

where

$$ $a^{[2]}_{l}=\begin{cases}\left(\hat{\bm{\Sigma}}_{-1-1}^{-1}\right)_{[2],l},& \text{for }l\in\{1,\dots,p-1\}\\ \sum_{q=1}^{n}\left(\hat{\bm{\Sigma}}_{-1-1}^{-1}\right)_{[2],l}\left(\beta_{1 }\right)_{q},&\text{for }l\in\{p,\dots,p^{2}-p\}\end{cases}$ (65) $$

$vec(B^{i}_{-1})$ is the row-wise vectorization of matrix $B^{i}_{-1}$, i.e.,

$$ $vec(B^{i}_{-1})=\begin{pmatrix}\xi^{i}_{2,1}\\ \xi^{i}_{3,1}\\ \vdots\\ \xi^{i}_{p,p}\end{pmatrix}\in\mathbb{R}^{p^{2}-p}.$ (66) $$

Thus, the distribution of $\hat{\beta}_{j,k}-\beta_{j,k}$ is

$$ $\hat{\beta}_{j,k}-\beta_{j,k}\sim N(0,{a^{[k]}}^{T}\frac{1}{n^{2}}\sum_{i=1}^{ n}vec(B^{i}_{-j})vec(B^{i}_{-j})^{T})a^{[k]}).$ (67) $$

In practice, we can plug in the estimates of $\beta_{j}$ to estimate the interested distribution and do the CI test by hypothesizing $\beta_{j,k}=0$.

### A.9 Discussion of assumption of zero mean and identity variance

In this section, we engage in a more thorough discussion regarding our assumptions about $\bm{X}$. Specifically, we demonstrate that this assumption of mean and variance does not compromise the generality. In other words, the true model may possess different mean and variance values, but we proceed by treating it as having a mean of zero and identity variance.

The key ingredient allowing us to assume such a model is, the discretization function $g_{j}$ is an unknown nonlinear monotonic function. Suppose the $g_{j}^{\prime}$ maps the continuous domain to a binary variable, and we have the "groundtruth" variable, denoted $X_{j}^{\prime}$, with mean $a$ and variance $b$. Assume the cardinality of the discretized domain is only 2, i.e., our observation $\tilde{X}_{j}$ can only be 0 or 1.
We further have the constant $d_{j}^{\prime}$ as the discretization boundary such that we have the observation

$$ $\tilde{X}_{j}=\mathbbm{1}(g_{j}^{\prime}(X_{j}^{\prime})>d_{j}^{\prime})= \mathbbm{1}(X_{j}^{\prime}>g_{j}^{\prime-1}(d_{j}))$ $$

We can always produce our assumed variable $X_{j}$ with mean 0 and variance 1, such that $X_{j}=\frac{1}{\sqrt{b}}X_{j}^{\prime}-\frac{a}{\sqrt{b}}$ and the same observation with a different nonlinear transformation $g_{j}$ and decision boundary $d_{j}$, such that

$$ $\tilde{X}_{j}=\mathbbm{1}(g_{j}(X_{j})>d_{j})=\mathbf{1}(X_{j}>g_{j}^{-1}(d_{j }))=\mathbbm{1}(X_{j}^{\prime}>\sqrt{b}g_{j}^{-1}(d_{j})+a)$ $$

As long as the observation $\tilde{X}_{j}$ is the same, we should have $\sqrt{b}g_{j}^{-1}(d_{j})+a=g_{j}^{\prime-1}(d_{j})$. Our assumed model $X_{j}$ clearly mimics the "groundtruth" $X_{j}^{\prime}$. Besides, according to Lem. A.2, we have one-to-one mapping between $\hat{\tau}_{j_{1}j_{2}}$ with the estimated covariance for fixed $\hat{h}_{j_{1}},\hat{h}_{j_{2}}$. Thus, as long as the observation is the same, the estimation of covariance $\hat{\sigma}_{j_{1},j_{2}}$ remains unaffected by our assumptions regarding the mean and variance of $\bm{X}$, so do the following inference.

We further conduct casual discovery experiments to empirically validate our statement, which is shown in App. C.3.

## Appendix B Data Generation and Figure of main experiments: causal discovery

##### Data Generation
and Code

We construct the true DAG $\mathcal{G}$ using the Bipartite Pairing (BP) model (Asratian et al., ), with the number of edges being one fewer than the number of nodes. The subsequent generation of true multivariate Gaussian data involves assigning causal weights drawn from a uniform distribution $U\sim(0.5,2)$ and incorporating noise via samples from a standard normal distribution for each variable.
Following this, we binarize the data, setting the threshold randomly based on each variable’s range. The code implementation is based on (Zheng et al., )
.

## Appendix C Additional experiments

### C.1 Linear non-Gaussian and nonlinear

Our model requires that the original data must adhere to the hypothesis of following a multivariate normal distribution, which appears to potentially limit the generalizability. Therefore, it is worthwhile to explore its robustness when such assumptions are violated. In this regard, we conducted several experiments, including scenarios involving linear non-Gaussian and nonlinear Gaussian.

For both cases, we follow the setting of our experiment where there are $p=8$ nodes and $p-1$ edges. We explore the effect of changing sample size $n=(100,500,2000,5000)$. Specifically for linear non-Gaussian case,
we adhere to some of the settings outlined by (Shimizu et al., ), conducting experiments where the original continuous data followed: (1) a Student’s t-distribution with 3 degrees of freedom, (2) a uniform distribution, and (3) an exponential distribution. Each variable is generated as $X_{i}=f(PA_{i})+noise$, where $noise$ follows the distribution in (1), (2), (3) correspondingly and $f$ is an arbitrary linear function. The first three rows of Fig. LABEL:fig_linear_nongaussian_nonlinear and Fig. LABEL:fig_linear_nongauss_nonlinear_dir show the result of the linear non-Gaussian case.

For the nonlinear cases, we follow setting in (Li et al., ), where every variable $X_{i}$ is generated as $X_{i}=f(WPA_{i}+noise)$, $noise\sim N(0,1)$ and $f$ is a function randomly chosen from (a) $f(x)=sin(x)$, (b) $f(x)=x^{3}$, (c) $f(x)=tanh(x)$, and (d) $f(x)=ReLU(x)$. $W$ is a linear function. Similarly, we set the number of nodes at $p=8$ and change the number of samples $n=(500,2000,5000)$. For both cases, we run 10 graph instances with different seeds and report the result of skeleton discovery in Fig. LABEL:fig_linear_nongaussian_nonlinear and DAG in Fig. LABEL:fig_linear_nongauss_nonlinear_dir (The same orientation rules (Dor & Tarsi, ) used in the main experiment are employed to convert a CPDAG (Chandler Squires, ) into a DAG). The last row of Fig. LABEL:fig_linear_nongaussian_nonlinear and Fig. LABEL:fig_linear_nongauss_nonlinear_dir shows the result of the nonlinear case.

Based on the experimental outcomes, DCT demonstrates marginally superior or comparable efficacy in terms of the F1-score, precision, and SHD relative to both the Fisher-Z test and the Chi-square test when dealing with small sample sizes. Nevertheless, as the sample size increases, DCT’s performance clearly surpasses that of the aforementioned tests across all three evaluated metrics, especially in the linear case. Consistent with observations from the main experiment, DCT exhibits a lower recall in comparison to the baseline tests. This discrepancy can be attributed to the baseline tests being prone to incorrectly infer conditional dependence and connect a large proportion of nodes. According to the results, our test shows notable robustness under the case assumptions are violated, confirming its practical effectiveness.

### C.2 Denser graph

DCT primarily works on cases where CI is mistakenly judged as conditional dependence due to discretization. Consequently, its efficacy is more pronounced in scenarios characterized by a relatively sparse graph, as numerous instances are truly conditionally independent. Nevertheless, the investigation of causal discovery with a dense latent graph is essential for evaluating the power of a test, i.e., its ability to successfully reject the null hypothesis when the tested pairs are conditionally dependent.
Thus, we conduct the experiment where $p=8,n=10000$ and changing edges $(p+2,p+4,p+6)$. Similarly, the latent continuous data follows a multivariate Gaussian model and the true DAG $\mathcal{G}$ is constructed using BP model. We run 10 graph instances with different seeds and report the result of the skeleton discovery and DAG in Fig. LABEL:fig_dense_graph.

According to the experiment results, DCT exhibits better performance in terms of the F1-score, precision, and SHD relative to both the Fisher-Z test and the Chi-square test.
As the graph becomes progressively denser, the superiority of the DCT correspondingly diminishes as there are few conditional independent cases in the true DAG. Due to the same reason, The recall remains lower than that of other baseline methods.

### C.3 Multivariate Gaussian with nonzero mean and non-unit variance

We employed a setting nearly identical to the main experiment, with the only difference being the alteration in data generation: instead of using a standard normal distribution, we used a Gaussian distribution with mean sampled from $U(-2,2)$ and variance sampled from $U(0,3)$.
We fix the number of variables as $p=8$ and change the number of samples $n=(100,500,2000,5000)$. The Fig. LABEL:fig_nonmean shows the result and demonstrates the effectiveness of our method.

### C.4 Real-world dataset

To further validate DCT, we employ it on a real-world dataset: Big Five Personality [https://openpsychometrics.org/](https://openpsychometrics.org/_rawdata/), which includes 50 personality indicators and over 19000 data samples. Each variable contains 5 possible discrete values to represent the scale of the corresponding questions, where 1=Disagree, 2=Weakly disagree, 3=Neutral, 4=Weakly agree and 5=Agree, e.g., "N3=1" means "I agree that I worry about things". This scenario clearly suits DCT, where the degree of agreement with a certain question must be a continuous variable while we can only observe the result after categorization. We choose three variables respectively: [N3: I worry about things], [N10: I often feel blue
], [N4: I seldom feel blue]. We then do the casual discovery using PC algorithm with DCT and compare it with the Chi-square test and Fisher-Z test. The result can be found in
Fig. LABEL:real_world.

Based on the experimental outcomes, despite the absence of a groundtruth for reference, we observe that the results obtained via DCT appear more plausible than those derived from Fisher-Z and Chi-square tests. Specifically, DCT suggests the relationship $N_{3}\perp\!\!\!\perp N4|N10$, which is reasonable as intuitively, the answer of ’I often feel blue’ already captures the information of ’I seldom feel blue’. As a comparison, both Fisher-Z and Chi-square return a fully connected graph. The results directly correspond to our illustrative example shown in Fig. LABEL:fig:1, substantiating the necessity of our proposed test.

## Appendix D Related Work

Testing for CI is pivotal in the field of causal discovery (Spirtes et al., ), and
a variety of methods exist for performing CI tests (CI tests). An important group of CI test methods involves the assumption of Gaussian variables with linear dependencies. For example, under this assumption, Gaussian graphical models are extensively studied (Yuan & Lin, ; Peterson et al., ; Mohan et al., ; Ren et al., ). To address CI test under Gaussian assumption, partial correlation serves as a viable method for CI testing (Baba et al., ).
To evaluate the independence of variables $X_{1}$ and $X_{2}$ conditional on $\bm{Z}$, The technique proposed by (Su & White, ) determines CI by comparing the estimations of $p(X_{1}|X_{2},\bm{Z})$ and $p(X_{1}|X_{2})$.

Another approach involves discretizing $\bm{Z}$ and performing independent tests within each resulting bin (Margaritis, ). Our work, however, diverges from these existing methods in two significant ways. Firstly, we are equipped to handle data, where partial variables are discretized. Additionally, we postulate that discrete variables are derived from the transformation of continuous variables in a latent Gaussian model. With the same assumption, the most closely related study is by (Fan et al., ), where the authors developed a novel rank-based estimator for the precision matrix of mixed data. However, their work stops short of providing a CI test for this method. Our research fills this gap, offering the ability to estimate the precision matrix for both discrete and mixed data and providing a rigorous CI test for our methodology.

Recent advancements in CI testing have utilized kernel methods for continuous variables influenced by nonlinear relationships. (Fukumizu et al., ) describes non-parametric CI relationships using covariance operators in reproducing kernel Hilbert spaces (RKHS). KCI test (Zhang et al., ) assesses the partial associations of regression functions linking $x$, $y$, and $z$, while RCI test (Strobl et al., ) aims to enhance the KCI test’s efficiency. In KCIP test (Doran et al., ) employs permutations of samples to emulate CI scenarios. CCI test (Sen et al., ) further reformulates testing into a process that leverages the capabilities of supervised learning models. For discrete variable analysis, the $G^{2}$ test (Aliferis et al., ) and conditional mutual information (Zhang et al., ) are commonly employed.
However, their method cannot deal with our setting where only discretized version of latent variables can be observed.

## Appendix E Resource Usage

All the experiments are run using Intel(R) Xeon(R) CPU E5-2680 v4 with 55 processors. It costs 4 hours to run experiments in Section 3.1.

## Appendix F Limiation and Broader Impacts

##### Limitation

So far, the largest limitation of our method is to treat discretized variables as binary, which wastes the available information. Besides that, the parametric assumption limits its generalizability. However, we need to point out this is pretty normal in CI test fields.

##### Broader Impacts

The goal of our proposed method is to test the conditional independence relationship given discretized observation. This task is essential and has broad applications. We are confident that our method will be beneficial and will not result in negative societal impacts.

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Please see the above comments.
soundness: 3
strengths: The problem studied is highly relevant and significant across various fields. 
The theoretical findings and asymptotic analysis provide valuable insights, and the empirical results are promising.
summary: The authors examine the problem of developing conditional independence (CI) tests when a subset of variables is discretized. This discretization introduces bias into existing CI test methods. They apply the concept of bridge equations to accommodate discretized random variables. Additionally, they demonstrate that the asymptotic distribution of their test under the null hypothesis (conditional independence) is normal and finally evaluate their test empirically.
weaknesses: The method is developed under the assumption of Gaussianity, which is quite limiting. Although the authors argue that this assumption is not restrictive, I find their reasoning unconvincing. The analysis relies on correlation measures (covariance and precision matrices), which only reveal dependency and independence structures in the case of multivariate normal distributions.

Additionally, Equation (5) lacks precision, as a Taylor expansion includes a residual term. Specifically, the last term should be $f'(\theta_*)$  for a $\theta_*$ between $\theta_0$ and $\hat{\theta}$ rather than $f'(\theta_0)$.

The empirical results are limited for showing the generality of the developed method for non-gaussian settings.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. In Theorem 2.1, the authors emphasize the importance of establishing correct CI relationships among the latent continuous variables of interest. However, if we can only observe $\tilde X_2$, why is it necessary to infer the original conditional relationships? Could the authors provide practical examples to clarify this necessity?

2. In line 192, multiple definitions for $\tilde h_j$ are possible, such as substituting $E\tilde X_j$ with other quantities. Considering that $\tilde h_j$ impacts the asymptotic variance of the estimated covariance as noted in Theorem 3.5, is there an optimal choice for $\tilde h_j$?
soundness: 4
strengths: This paper is well-written and easy to follow. The proposed CI test effectively corrects and infers the true CI relationships among latent continuous variables, recognizing that discretized variables may show different conditional independence compared to their continuous counterparts.
summary: This paper tackles the challenge of testing conditional independence (CI) with a focus on accommodating discretized data, a crucial aspect for various applications. The authors introduce bridge equations to estimate covariance and establish asymptotic normality, facilitating an unconditional independence test. For the conditional independence test, they employ nodewise regression to recover precision coefficients. The effectiveness of the proposed method is substantiated through rigorous theoretical analysis and empirical validation.
weaknesses: The authors primarily discretized variables into binary categories, which could potentially lead to some efficiency loss.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 2
questions: Major questions:
1. The authors empirically verify the Type I and Type II errors of the test. I’m curious if the asymptotic Type I error and power of the test could also be theoretically verified as well. I'm assuming the asymptotic power or Type I error will depend on $\mathbf{Λ}$ as well.
2. Have the authors considered providing a quantitative rate theory to support the test’s validity for smaller datasets? Additionally, are there any debiasing corrections that could be applied for small dataset sizes?
3. [Optional] The authors present a comprehensive framework based on the assumption that the random vector follows a multivariate normal distribution. I wonder if they have any thoughts on extending this framework to accommodate non-normal variables.

Minor questions:
1. Line 34: undefined notation $X_{j_1}$. More explanation on what ${j_1}$ means or may consider using $X_{1}$ and $X_{2}$ here.  
2. Line 89-90: same issue as in 1.
3. Line 122: inconsistent notations $\mathrel{\text{{$\perp\mkern-10mu\perp$}}}$ versus $\bot$ defined at Line 36. The former notation is more commonly referred as conditional independence in literature.
4. Line 126: what does DCT stand for? Should introduce the abbreviation beforehand.
5. Line 754: minor typo, should be $X_1\mathrel{\text{{$\perp\mkern-10mu\perp$}}}X_3 | X_2$
6. Line 137: undefined subscript $\sigma_{j_1, j_2}$. Consider using different subscripts for $\sigma_{j_1, j_2}$ and $\omega_{k,l}$ to distinguish them.
7. Line 362: remove ) after $vec(B^i_{-j})^T$.
8. Line 369: as much as I love the theory, I do find the mathematical notations concerning. I don't really see people using $vec$ to denote the rows of a matrix.
soundness: 3
strengths: The paper presents a well-established theoretical foundation and methodology. The authors begin with a straightforward theorem to illustrate the impact of discretization, emphasizing the potential need for correction in CI testing. The theoretical framework is clearly articulated and thorough.
summary: This paper addresses a specific challenge in confidence independence (CI) testing when variables of interest are discretized. In practice, accurately measuring continuous variables is often difficult due to data collection limitations, resulting in a common reliance on discretized values for these variables. To tackle this, the authors propose a unified bridge equation to estimate the covariance between two random variables by binarizing them. They then employ nodewise regression to recover precision coefficients that capture the conditional dependence among latent continuous variables. An appropriate test statistic is proposed, and its asymptotic distribution under the null hypothesis is derived.
weaknesses: The paper needs improvement in writing. There are several instances where the notations are unclear, inconsistent, and occasionally used imprecisely (see Questions for more details).

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1- In the numerical results for the causal discovery example, In Figure 3, with increasing the number of samples F1 score, and precision are reduced for Fisher_z and Chis-q. Usually, increasing the number of samples increases accuracy. What is your explanation? 

2-  The proposed CI test was developed and analytically proven for the Gaussian distribution. Although it has been applied to other distributions in the numerical results, can the proposed method be analytically extended to other continuous distributions or to handle nonlinear relationships?
soundness: 3
strengths: 1- Improving CI test accuracy is an essential and practical problem, particularly for structure learning in Bayesian networks and causal discovery.

2- The proposed method is novel, and the solution provided is relatively simple and efficient.
summary: In this paper, a CI test has been designed specifically for handling the presence of discretization. An appropriate test statistic for the CI of latent continuous variables, based solely on discretized observations, has been derived. The key is to build connections between the discretized observations and the parameters needed for testing the CI of the latent continuous variables. To achieve this, first, the bridge equations have been developed to estimate the covariance of the underlying continuous variables with discretized observations. Then, a node-wise regression has been applied to derive appropriate test statistics for CI relationships from the estimated covariance. By assuming that the continuous variables follow a Gaussian distribution, it has been derived the asymptotic distributions of the test statistics under the null hypothesis of CI.
weaknesses: 1- In (1), the discretization model (function) is unclear. How many discretization levels are used in this formulation? In addition, While the main idea is well presented, the authors have not clearly conveyed the details of the proposed method.

2- Providing a step-by-step explanation of CDT in a numerical example would enhance understanding of the proposed method.

3- The proposed method has improved the accuracy of the CI tests. However, another significant challenge with CI tests is the computational burden, particularly when there are many conditioning variables and a large sample size. Including a complexity analysis or a runtime comparison with existing methods in the numerical results would help demonstrate the impact of the proposed method.

4- Does not the number of discretizing levels affect the accuracy (Type I and Type II error) of the DCT method? In Figure 2 the effect of the number of discretizing levels on Fisherz and Chi2 methods has been checked and shown but the impact of it on CDT has not been reported. Why? I think the number of discretizing levels affects the accuracy of the CDT.  

5- Using practical examples for the application of the CDT in causal discovery can help to indicate the effectiveness of the proposed method in practical and real problems. 

6- In proof of Theorem 2.1, in appendix A.1, “Var(X2| \tilde{X2}) will never be zero since  \tilde{X2} is the discretized version of X2 (information loss).” is intuitively correct but it is not a precise mathematical expression for proof. Is there any reference or more precise explanation for that?

7- According to Figures 5, 6, and 8, CDT outperforms other methods in F1 score and Precision but generally performs worse in Recall.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 2
questions: Minor comments: 

line 73) This is a very interesting idea to use a DAG to invest the discretization problem; however, the authors implicitly assume that the DAG is faithful. The fact that $X_1$ is not d-separated of $X_3$ given $\tilde X_2$in the DAG does not imply that the $X_1\not\perp X_3\mid \tilde X_2$ in the distribution. 
line 73) extra comma
line 89) 1). 2). 3). are not common. Consider 1) 2) 3)
line 100) So all discretized variables are Gaussian? This is a limitation.
[important] line 137) The notation seems to be redundantly complex. Why \sigma_{j_1,j_2} rather than just \sigma_{j,k}, or even simpler \sigma_{jk}?
line 159) better not to start a sentence with an equation.
line 264) Consider using ^\top for transpose.
- Was DCT ever defined? What does it stand for?
line 264) grammar “contain”
line 269) “null hypothesis”
line 270) Be more specific, what are “our settings” and “our framework”?
soundness: 3
strengths: The studied problem is important. With respect to the presentation, the problem is well motivated, formulated, and explained, and the experimental result is also written clearly. Some good effort is put into the experiments.
summary: The paper addresses errors in existing CI tests caused by an implicit discretization of the original continuous variable during measurements. The authors propose DCT to deal with this issue and establish the necessary theoretical foundations under assumptions such as having a Gaussian distribution. The method is tested using simulations.
weaknesses: 1) The main drawback of this paper is the overall presentation. I was not able to grasp the main idea, preventing me from properly evaluating this work. The following provides some details:

a) The authors were careful and patient in explaining detailed notations up to line 190. However, after that, the pace and notational complexity suddenly got of hand. For example, in line 192, the notation 1{ } is not even defined, let alone the meaning of the equation (the sentence in line 194 does not help). It was as if this part of the paper was written by a different author. I had a hard time understanding the paper from this point.

b) Definition 3.1. What is “the estimated boundary”? It should be rigorously defined somewhere. If the intention of this definition is to define “the estimated boundary” then the language is not appropriate. It should be something like this: “the \emph{estimated boundary of variables $X_j$ is defined as ….”. Currently, it reads as the estimated boundary is already defined before. 

c) Definition 3.4. Isn’t this simply the formula for the sample covariance? that is Cov(X,Y) = 1/n \sum_{i=1}^n (X_iY_i) - \mu_X\mu_Y? If so, what’s the point of giving it a separate definition and name (bridge equation)?

d) Theorem 3.5. I am confused by what is the rigorous statement in this theorem? By “this statistic is approximated to follow a normal distribution”, do the authors mean that the statistic converges to a normal distribution as N grows? 

e) Theorem 3.9. Where is the \hat{\beta}_{jk} formally defined? I suppose the answer is not eq 13 as then it should have been out of the theorem. 

f) The Experiment Section is back on track and written clearly. It was unclear though why the addition of W at line 450 is good for assessing the power. 

2) One of the main motivations of this work was its application in Bayesian network (BN) structure learning. Nevertheless, insufficient effort was put into it. As the authors stated, the experiments in Section 4.2 are only for [extremely] sparse DAGs. They do consider denser ones in C2, but that is an only minor improvement, as the number of edges increased from p-1 to at most p+6. True DAGs can be much denser. The authors should provide a fair comparison for this part (also not to move these important parts in the appendix). 

a) I suggest to test several established datasets in BN structure learning such as ASIA, CHILD, WATER, and ANDES. 

b) Similar to the other experiment section, here, non-Gaussian distributions should be tested for structure learning. 

c) Even if the results here are weak, DCT appears to be a promising method that is worthy of publication, provided a major rewriting. 

d) Note that without this structure learning application, this paper perhaps better suits a statistics journal.

---
Overall, a good problem is studied in this paper with a seemingly effective solution. However, the presentation of the technical parts needs a substantial improvement.

Meta-review: Since there is a large span for the ratings of the reviewers (but 4 out of 5 judge the content "good"), I have read the paper and went through the reviews, the rebuttal stage, the updates in the paper.

The key part in the negative opinion of i1w6 is the sparsity of DAG used. In response, the authors have performed a lightweight experiment on denser graphs, which shows quite unbalanced precision and recall compared to the other baselines (and a smaller F1). The key argument of the authors, however, is that denser graphs should not be taken as the focus of the paper. More than the additional experiments, I am happy to consider this argument. I also note that the authors have made a *very* substantial effort of writing / rewriting parts of the paper as a consequence on the rebuttal process, which made reviewer WY1e bump their score up by 3 points.

i1w6's suggestion to submit to a statistical journal is reasonable given the paper's content, but it does not prevent the paper to have merits for ICLR as well, and I do consider it has such merits, mainly because of the broad interest of the question tackled -- and not just for causal inference --.

In conclusion, I believe the authors have made a very substantial effort to answer the reviewers comments and can be taken into consideration for acceptance.

Decision: Accept (Poster)