# Non-parametric Conditional Independence Testing for Mixed Continuous-
Categorical Variables: A Novel Method and Numerical Evaluation

Oana-Iuliana Popescu1, Andreas Gerhardus1, Jakob Runge1,2

###### Abstract

Conditional independence testing (CIT) is a common task in machine learning,
e.g. for variable selection, and a main component of constraint-based causal
discovery. While most current CIT approaches assume that all variables are
numerical or all variables are categorical, many real-world applications
involve mixed-type datasets that include numerical and categorical variables.
Non-parametric CIT can be conducted using conditional mutual information (CMI)
estimators combined with a local permutation scheme. Recently, two novel CMI
estimators for mixed-type datasets based on k-nearest-neighbors (k-NN) have
been proposed. As with any k-NN method, these estimators rely on the
definition of a distance metric. One approach computes distances by a one-hot
encoding of the categorical variables, essentially treating categorical
variables as discrete-numerical, while the other expresses CMI by entropy
terms where the categorical variables appear as conditions only. In this work,
we study these estimators and propose a variation of the former approach that
does not treat categorical variables as numeric. Our numerical experiments
show that our variant detects dependencies more robustly across different data
distributions and preprocessing types.

## Introduction

Conditional independence testing (CIT) is a central component of constraint-
based causal discovery frameworks; e.g., in algorithms such as PC and FCI
(Spirtes et al. 2000), and is used to infer causal relations from purely
observational data. The performance of causal discovery algorithms depends
heavily on the performance of the CIT and its robustness toward different
types, distributions, and sample sizes of the data. A good CIT approach
achieves high statistical power to detect true conditional dependence while
simultaneously controlling false positives at the desired level. Most current
CIT approaches assume that either all variables are numerical or that all
variables are categorical. Still, many real-world applications involve mixed-
type datasets, e.g., datasets with variables such as gender and height in
medicine or weather regime types and continuous temperature in climate
science.

In this work, we consider the two recent $k$-NN estimators of CMI for mixed-
type data. The estimator of Mesner and Shalizi (2019) transforms categorical
variables by one-hot encoding and then measures distances on the resulting
product space of mixed continuous-discrete variables. Zan et al. (2022) frame
mixed-type datasets as consisting of quantitative and qualitative variables
and rewrite the CMI as a linear combination of entropies where the qualitative
variables appear as conditions only, thus requiring distance notions on the
quantitative subspaces. CMI is of relevance for non-parametric CIT since the
conditional independence $X\perp\\!\\!\\!\perp Y~{}|~{}Z$ holds if and only if
$I(X;Y|Z)=0$ (Gray 2011) and thus is a non-parametric measure for conditional
independence. To construct CIT, CMI estimators can be combined with a local
permutation scheme (similar to Runge (2018)) and be formulated as a
statistical test of the null hypothesis $H_{0}:\,X\perp\\!\\!\\!\perp Y\mid Z$
as done in Zan et al. (2022). No general analytical results for the finite or
asymptotic distribution of CMI under conditional independence are known.

We study CIT based on the previously mentioned estimators from an empirical
perspective. We first outline the challenges these estimators face, some of
which have been discussed in the respective works, for example, the curse of
dimensionality that leads to increased bias. We briefly discuss why the CMI
estimation performance can affect the outcomes of the corresponding CIT. We
investigate how the two estimators and their respective CITs perform under
different choices of hyperparameters, data distributions, and combinations of
variable types and dimensionalities. To reduce the effect of the challenges
that the two estimators face on CIT, we propose a variant based on the
estimator of Mesner and Shalizi (2019) that does not rely on one-hot encoding
of categorical variables. In summary, our main contributions are (1) a new
$k$-NN estimator for the CMI of mixed-type data that is a variant of the
estimator of Mesner and Shalizi (2019), (2) an empirical evaluation of the
three CMI estimators, and (3) an extensive and systematic numerical evaluation
of CIT performance based on the three CMI estimators in combination with a
local permutation scheme.

## Background and related work

After the preliminaries, we give an extensive summary of related work in order
to highlight the subtle differences between the different estimators and to
build upon further below.

### Preliminaries

Let $X:\Omega\to\mathcal{X}$, $Y:\Omega\to\mathcal{Y}$ and
$Z:\Omega\to\mathcal{Z}$ be (vectors of) random variables with
$dim(\mathcal{X})=m_{X}$, $dim(\mathcal{Y})=m_{Y}$, and
$dim(\mathcal{Z})=m_{Z}$. We demand that
$\mathcal{X}=\mathcal{X}_{1}\times\ldots\mathcal{X}_{m_{X}}$ where
$\mathcal{X}_{a}$ with $1\leq a\leq m_{X}$ is $\mathbb{R}$ or a discrete set;
similarly for $\mathcal{Y}$ and $\mathcal{Z}$. Let $P_{XYZ}$ be the
probability measure on $\mathcal{X}\times\mathcal{Y}\times\mathcal{Z}$ induced
by the joint vector $(X,Y,Z)$. We assume that the conditional probability
measure $P_{XY|Z}$ exists and is absolutely continuous wrt to the product
measure $P_{X|Z}\times P_{Y|Z}$. As discussed in Mesner and Shalizi (2019),
these assumptions are fulfilled if every component of $(X,Y,Z)$ is either
discrete (i.e., absolutely continuous wrt to the counting measure) or non-
singular continuous (i.e., absolutely continuous wrt to the Lebesgue measure)
or a mixture of these two cases (for simplicity, from here on we refer to
“non-singular continuous” as “continuous”). The CMI $I(X;Y|Z)$ of $X$ and $Y$
given $Z$ can be defined as, see Gray (2011),

$I(X;Y|Z)=\int\log\left(\frac{dP_{XY|Z}}{d(P_{X|Z}\times
P_{Y|Z})}\right)dP_{XY|Z}\,,$ (1)

where the argument of the $\log$ is the Radon-Nikodym derivative of $P_{XY|Z}$
wrt $P_{X|Z}\times P_{Y|Z}$. If all components of $(X,Y,Z)$ are discrete, then
the rhs of eq. (1) reduces to the familiar form

$\sum_{x,y,z}p_{XYZ}(x,y,z)\log\frac{p_{XY|Z}(x,y|z)}{p_{X|Z}(x|z)p_{Y|Z}(y|z)}$
(2)

in terms of the probability mass functions (pmfs) $p_{\cdot}$. If all
components are continuous, then in eq. (2) integrals replace sums and
probability density functions (pdfs) $f_{\cdot}$ replace pmfs.

We distinguish three types of discrete random variables $V$: First, the values
of $V$ can be a discrete subset of $\mathbb{R}$ such that distance notions on
$\mathbb{R}$ (e.g., $L^{p}$-distance) other than the discrete metric are
semantically meaningful (“discrete numeric”). Second and third, the values of
$V$ can be on an ordinal or a nominal/categorical scale (“non-numeric”). From
a measure-theory perspective, these cases are equivalent. When speaking of
“mixed continuous-categorical variables” (in short: “mixed variables”), we
refer to any of the following cases: First, all of $X$, $Y$, $Z$ are either
fully discrete or fully continuous, with at least one of them fully discrete
and another fully continuous. Second, at least one of $X$, $Y$, $Z$ contains
both a discrete and a continuous component, but no component of $(X,Y,Z)$ is a
mixture variable. Third, at least one component $V$ of $(X,Y,Z)$ is a mixture
variable; that is, this component $V$ itself is neither discrete nor
continuous. In this paper, we focus on the first two of these cases.

### CMI estimation in the fully continuous case

KL estimator for differential entropy Let $W:\Omega\to\mathcal{W}$ be a
(vector of) continuous random variables with $\mathcal{W}=\mathbb{R}^{m_{W}}$
and let $w_{1},\,\ldots,w_{n}$ be _iid_ observations of $W$. The Kozachenko
and Leonenko (1987) (KL) estimator of the differential entropy $H(W)$ is the
sample average

$\hat{H}^{KL}(W)=-\frac{1}{n}\sum_{i\,=\,1}^{n}\log\hat{f}_{W}(w_{i}).$ (3)

The local density estimates$\hat{f}_{W}(w_{i})$ are calculated under the
assumption that $f_{W}$ is locally constant within an $L^{p}$-ball
$B(w_{i},\rho_{i})$ of radius $\rho_{i}$ around $w_{i}$ where $\rho_{i}$ is
the $L^{p}$-distance of $w_{i}$ to its $k$-th nearest neighbor (not counting
$w_{i}$ itself) for some positive integer $k$. Since $W$ is continuous, this
$k$-th nearest neighbor is unique with probability one. The local constancy
assumption implies that the probability $P_{i}$ of the event $w\in
B(w_{i},\rho_{i})$ is $P_{i}\equiv p_{w}(w_{i})\cdot
V_{m_{W},p}\cdot\rho_{i}^{m_{W}}$, where $V_{m_{W},p}$ is the volume of the
unit-ball in the $L^{p}$-metric. Using that $\mathbb{E}[\log
P_{i}]=\psi(k)-\psi(n)$ with the digamma function $\psi(x)$, see Kraskov,
Stögbauer, and Grassberger (2004), and approximating $E[\log\rho_{i}]$ by a
sample average, eq. (3) then takes the form

$\hat{H}^{KL}(W)=\psi(n)-\psi(k)+\log
V_{m_{W},p}+\frac{m_{W}}{n}\sum_{i\,=\,1}^{n}\log\rho_{i}\,.$ (4)

KSG estimator for mutual informationKraskov, Stögbauer, and Grassberger (2004)
estimate the MI $I(X;Y)=H(X)+H(Y)-H(X,Y)$ by estimating the three individual
entropies with the KL estimator (4). The authors heuristically argue that the
errors incurred by the local constancy assumptions approximately cancel out in
the combined estimator if all entropy estimates use the same local length
scales $\rho_{i}$. Thus, they equip $\mathcal{X}\times\mathcal{Y}$ with the
maximum metric
$d_{\mathcal{X}\times\mathcal{Y}}(\cdot,\cdot)=\max\\{d_{\mathcal{X}}(\cdot,\cdot),\,d_{\mathcal{Y}}(\cdot,\cdot)\\}$
and define $\rho_{i}$ as in the KL estimate (4) of $H(X,Y)$, and use the same
radii $\rho_{i}$ for estimating $H(X)$ and $H(Y)$. The estimator takes the
form

$\hat{I}^{KSG}(X;Y)=\frac{1}{n}\sum_{i\,=\,1}^{n}[\psi(k)+\psi(n)-\psi(k_{X,i}+1)-\\\
-\psi(k_{Y,i}+1)]\,,$ (5)

where $k_{X,i}$ and $k_{Y,i}$ are defined by ($W$ is placeholder for $X$ and
$Y$ and $w$ is placeholder for $x$ and $y$)

$k_{W,i}=|\\{w_{j}~{}|~{}\|w_{j}-w_{i}\|<\rho_{i},\,j\neq i\\}|\,,$ (6)

as the number of points $x_{j}\neq x_{i}$ (resp. $y_{j}\neq y_{i}$) within the
_open_ ball $B(x_{i},\rho_{i})$ (resp. $B(y_{i},\rho_{i})$) in $\mathcal{X}$
(resp. $\mathcal{Y}$). The terms with $\log\rho_{i}$ cancel out due to using
the same radii $\rho_{i}$ in all three entropy estimates. Since
$\mathcal{X}\times\mathcal{Y}$ is equipped with the maximum metric, the volume
terms cancel out too.

FP estimator of conditional mutual information Using the same rationale,
(Frenzel and Pompe 2007) extend the KSG estimator to CMI, with $k_{Z,i}$,
$k_{XZ,i},k_{YZ,i}$ as in eq. (6) and $\rho_{i}$ as in the KL estimate of
$H(X,Y,Z)$:

$\hat{I}^{FP}(X;Y|Z)=\frac{1}{n}\sum_{i\,=\,1}^{n}[\psi(k)+\psi(k_{Z,i}+1)-\\\
-\psi(k_{X,i}+1)-\psi(k_{Y,i}+1)]\,.$ (7)

### CMI estimation in the mixed variables case

GKOV estimator of mutual information Gao et al. (2017) propose an estimator
for mixed MI $I(X;Y)$ under the assumption that both $\mathcal{X}$ and
$\mathcal{Y}$ are Euclidean spaces, thus implicitly requiring that either the
discrete values are numeric with a semantically meaningful notion of distance
or that they have been mapped to a real space (that is, ignoring the
conceptual problem of a semantically non-meaningful $L^{p}$-distance). The
GKOV estimator builds on KSG and the observation that, in the mixed case, the
distance $\rho_{i}$ of $(x_{i},y_{i})$ to its $k$-th nearest neighbor in
$\mathcal{X}\times\mathcal{Y}$ can be $\rho_{i}=0$ with non-zero probability.
Gao et al. (2017) consider the event $\rho_{i}=0$ to indicate that point
$(x_{i},y_{i})$ is “discrete”. Their estimator takes the form

$\hat{I}^{GKOV}(X;Y)=\frac{1}{n}\sum_{i\,=\,1}^{n}[\psi(\tilde{k}^{\prime}_{i})+\log(n)-\\\
-\log(\tilde{k}_{X,i}+1)-\log(\tilde{k}_{Y,i}+1)]\,,$ (8)

where $\tilde{k}^{\prime}_{i}=k$ if $\rho_{i}>0$ and
$\tilde{k}^{\prime}_{i}=\tilde{k}_{XY,i}$ if $\rho_{i}=0$ with

$\tilde{k}_{W,i}=|\\{w_{j}~{}|~{}\|w_{j}-w_{i}\|\leq\rho_{i},\,j\neq i\\}|\,.$
(9)

As opposed to eq. (6), eq. (9) uses the non-strict inequality
$\|w_{j}-w_{i}\|\leq\rho_{i}$. The combination of $\psi(\cdot)$ and
$\log(\cdot)$ terms is ad-hoc and ultimately justified by their consistency
proof.

MS estimator of conditional mutual information Mesner and Shalizi (2019)
propose an estimator that slightly modifies the CMI generalization of the GKOV
estimator. The modifications are motivated by the observation that, besides
$\rho_{i}=0$, there is also a non-zero probability that different pairs of
points have the same distance. Thus, the $k$-th nearest neighbor of
$(x_{i},y_{i},z_{i})$ is non-unique with non-zero probability, and a non-
unique $k$-th nearest neighbour is equivalent to $k<\tilde{k}_{XYZ,i}$ with
$\tilde{k}_{XYZ,i}$ as defined by eq. (9) for $W=XYZ$. Instead of
$\rho_{i}=0$, the event $k<\tilde{k}_{XYZ,i}$ is considered to indicate that
$(x_{i},y_{i},z_{i})$ is a “discrete” point. Specifically, their estimator
takes the form

$\hat{I}^{MS}(X;Y|Z)=\frac{1}{n}\cdot\\\
\sum_{i\,=\,1}^{n}\underbrace{\left[g(\tilde{k}_{XYZ,i})+g(\tilde{k}_{Z,i})-g(\tilde{k}_{XZ,i})-g(\tilde{k}_{YZ,i})\right]}_{\equiv\,\hat{\xi}^{MS}_{i}(X;Y|Z)}\,$
(10)

where $g(\cdot)=\psi(\cdot)$ if $\tilde{k}_{XYZ,i}=k$ and
$g(\cdot)=\log(\cdot)$ if $\tilde{k}_{XYZ,i}>k$.111In their paper, Mesner and
Shalizi (2019) define their estimator by additionally computing a maximum of
the estimate with $0$, motivated by the fact that $I(X;Y|Z)\geq 0$. However,
their implementation does seem not actually apply this, and preliminary
experiments of ours show that the maximum with 0 can be detrimental for CIT.
We do not apply the maximum with $0$ in any experiments. The authors prove
consistency of their estimator, and also show that it suffers from the curse
of dimensionality: For fixed $m_{X}$ and $m_{Y}$, if the dimension $m_{Z}$ of
$\mathcal{Z}$ increases to infinity and $H(Z)/m_{Z}$ is non-zero in this
limit, then $\hat{I}^{MS}(X;Y|Z)$ converges to $0$ in probability as
$m_{Z}\to\infty$. The MS estimator equips the discrete components of
$\mathcal{X}$, $\mathcal{Y}$, $\mathcal{Z}$ with the discrete metric, which is
equivalent to a one-hot encoding of the components and again raises the
conceptual problem that the corresponding distance notions might not be
semantically meaningful. In their experiments, the authors heuristically set
$k=n/10$ where $n$ is the sample size.

ZMADG estimator of conditional mutual information (Zan et al. 2022) assumes
the absence of mixture variables and proposes an estimator for this case of
mixed variables CMI $I(X;Y|Z)$ that avoids defining a distance between
qualitative components. They split $X$, $Y$, and $Z$ in their respective
quantitative components $X^{t},Y^{t},Z^{t}$ and qualitative components
$X^{l},Y^{l},Z^{l}$ and express the CMI as

$I(X;Y|Z)=H(X^{t},Z^{t}|X^{l},Z^{l})+H(Y^{t},Z^{t}|Y^{l},Z^{l})-\\\
-H(X^{t},Y^{t},Z^{t}|X^{l},Y^{l},Z^{l})-H(Z^{t}|Z^{l})+H(X^{l},Z^{l})+\\\
+H(Y^{l},Z^{l})-H(X^{l},Y^{l},Z^{l})-H(Z^{l})\,,$ (11)

where the first four terms on the rhs are (conditional) differential entropies
and the last four are (conditional) entropies. The (conditional) entropies are
estimated with the standard plug-in estimator using empirical frequencies,
while (conditional) differential entropies are calculated using the KL
estimator on each subset of the samples defined by fixed values of the
qualitative components and then averaging according to the empirical
frequencies of the qualitative values. The parameter $k$ of the KL estimates
is set to $k=\max\left\\{\lfloor n_{cluster}/10\rfloor,\,1\right\\}$ with
$n_{cluster}$ the number of samples in the respective subsets determined by
the values of the qualitative components (i.e., $k$ is separately chosen for
each subset of samples). As a sum of consistent estimators, the estimator is
consistent. The ZMADG estimator does not seem to suffer from the curse of
dimensionality as the MS estimator, but, as we will further discuss below, we
believe it incurs higher variance.

### Non-parametric CIT using CMI and a local permutation scheme

To statistically test the null hypothesis $H_{0}:\,X\perp\\!\\!\\!\perp
Y~{}|~{}Z$ of conditional independence from finite samples
$(x_{i},y_{i},z_{i})$, a distribution of the estimate $\hat{I}(X;Y|Z)$ under
the null hypothesis (the so-called null distribution) or an approximation
thereof is needed. If $X\perp\\!\\!\\!\perp Y~{}|~{}Z$, then the component
values $x_{i}$ and $y_{i}$ within the subset of samples determined by the
value $z_{i}$ can be permuted arbitrarily without changing the distribution of
the estimated CMI, i.e., setting $\tilde{x}_{i}=x_{\sigma(i)}$ with a
permutation $\sigma$ such that for all $i$ both $x_{i}$ and $x_{\sigma(i)}$
are in the subset of samples determined by $z_{i}$, the estimators
$\hat{I}(X;Y|Z)$ and $\hat{I}(\tilde{X};Y|Z)$ have the same distribution.
Since this equality holds for any such permutation, a null distribution can be
obtained. For fully discrete $Z$, the subset of samples determined by $z_{i}$
are all samples $(x_{j},y_{j},z_{j})$ with $z_{j}=z_{i}$. For fully continuous
$Z$, Runge (2018) uses a $k$-NN approach to determine the subsets of samples
for which $z_{j}\approx z_{i}$ according to the $L^{\infty}$-distance. Zan et
al. (2022) adapt this method to the mixed data case: The sample
$(x_{j},y_{j},z_{j})$ with $z_{j}=(z^{t}_{j},z^{l}_{j})$, where $z^{t}_{j}$ is
the quantitative and $z^{l}_{j}$ the qualitative component, is part of the
subset of samples determined by $z_{i}=(z^{t}_{i},z^{l}_{i})$ if and only if
$z^{l}_{j}=z^{l}_{i}$ and $z^{t}_{j}\approx z^{t}_{i}$. Supplementaryi2
Material (SM) Sec. A describes how p-values are obtained.

## Proposed novel estimator

We first discuss the problems of the MS and ZMADG estimators that motivate us
to introduce a novel estimator. We then formally define this novel estimator
and presents its theoretical guarantees.

### Motivation: Problems of the MS and ZMADG estimators

We highlight three issues of the MS estimator (Mesner and Shalizi 2019).
First, it suffers from the conceptual problem that—because the $k$-NNs can
come from different _clusters_ (defined as the subsets of samples points with
equal values of the discrete variable)—it implicitly assumes local constancy
across different clusters. However, different clusters might be entirely
unrelated to each other. For example, it could be that dependence exists in
only one of the clusters. Despite this fact, the MS estimator might estimate
the local contribution of a point by combining neighbours from both the
cluster with and without dependence. Not only does this give rise to
conceptual complications, but it can also be expected to negatively affect
statistical power. Second, due to the one-hot encoding of non-numeric discrete
variables, the MS estimator is (as opposed to CMI) not invariant under scaling
all variables with a common factor. Third, as discussed in Mesner and Shalizi
(2019), the MS estimator is biased towards $0$ in high-dimensional settings.
To exemplify, say the continuous and numeric discrete variables are scaled to
$[0,1]$ in preprocessing. Then, due to one-hot encoding and the
$L^{\infty}$-metric, the maximum distance between any two sample points is
$1$. Thus, if the clusters of the $i$-th sample point $w_{i}$ contains at most
$k$ points, then $\rho_{i}=1$ (since there are not enough points in the
cluster), which in turn implies $k_{i,XZ}=k_{i,YZ}=k_{i,XYZ}=k_{i,Z}=n$ and
hence $\xi^{MS}_{i}(X;Y|Z)=0$. Zan et al. (2022) discuss further cases in
which the MS estimators suffers from local zero estimates. Generally, a bias
towards zero can affect CI test performance because it can lead to false
conclusions of independence.

The ZMADG estimator Zan et al. (2022) reduces these problems by considering
each discrete cluster individually and adaptively reducing $k$ (in the
estimation of entropies). However, this approach can lead to another issue
that, unfortunately, has not yet been discussed or investigated in detail:
Since the CMI estimator is a sum of up to $8$ entropy estimators, the CMI
estimator might suffer from higher variance than the MS estimator, leading to
increased CIT error rates.

### Definition and intuition of the proposed novel estimator

To address these problems, we introduce a novel CMI estimator MS0-∞ that
combines ideas from the MS and ZMADG estimators. Specifically, MS0-∞ can be
understood as a variant of MS with the following two modifications.

First, instead of one-hot encoding non-numeric variables, we keep the original
space $\mathcal{X}\times\mathcal{Y}\times\mathcal{Z}$ and equip it with the
$0-\infty$ “metric”222Formally, the $0-\infty$ “metric” is not a metric due to
the value $+\infty$. However, we use this formulation to highlight the
similarity with MS. defined as

$||w_{i}-w_{j}||_{0-\infty}=\begin{cases}||w_{i,c}-w_{j,c}||_{L^{\infty}}&\text{if
}w_{i,d}=w_{j,d}\\\ \infty&\text{otherwise}\end{cases}\,,$ (12)

where we split the point $w_{k}=(x_{k},y_{k},z_{k})$ into its numeric
component $w_{k,c}=(x_{k,c},y_{k,c},z_{k,c})$ and its non-numeric component
$w_{k,d}=(x_{k,d},y_{k,d},z_{k,d})$; similarly for the subspaces
$\mathcal{X}\times\mathcal{Z}$, $\mathcal{Y}\times\mathcal{Z}$, and
$\mathcal{Z}$. That is, if $w_{i}$ and $w_{j}$ are in the same cluster (i.e.,
$w_{i,d}=w_{j,d}$), then their distance is finite and measured by the
$L^{\infty}$-distance, else their distance is $\infty$.

Second, we adopt the heuristic to adaptively set $k=\lfloor
k_{c}\cdot(n_{cl,\min}-1)\rfloor$, where $0<k_{c}<1$ is a hyperparameter and
$n_{cl,\min}=\min_{i\in[[1,n]]}|\\{w_{j}:||w_{i}-w_{j}||_{L^{\infty}}\neq\infty\\}|$
is the number of points in the “smallest” cluster. The necessity of such a
heuristic stems from the fact that, unlike in the infinite sample case, in
practice some clusters might contain less than $k+1$ points. In Sec. D of the
SM, we compare multiple heuristics and motivate our final choice.

To formally specify our estimator, we first define the counts

$\tilde{k}_{W,i}^{0-\infty}=|\\{w_{j}~{}|~{}\|w_{j}-w_{i}\|_{0-\infty}\leq\rho_{i},\,j\neq
i\\}|\,,$ (13)

where $\rho_{i}<\infty$ is the $0-\infty$ distance of $w_{i}$ to its $k$-th
nearest neighbour (here, since $\rho_{i}$ is finite, this distance equals the
$L^{\infty}$-distance) and $W$ stands for $XYZ$, $XZ$, $YZ$ or $Z$. In terms
of these counts, our estimator reads

$\hat{I}^{0-\infty}(X;Y|Z)=\frac{1}{n}\cdot\\\
\sum_{i\,=\,1}^{n}\underbrace{\left[g(\tilde{k}_{XYZ,i}^{0-\infty})+g(\tilde{k}_{Z,i}^{0-\infty})-g(\tilde{k}_{XZ,i}^{0-\infty})-g(\tilde{k}_{YZ,i}^{0-\infty})\right]}_{\equiv\,\hat{\xi}^{0-\infty}_{i}(X;Y|Z)},$
(14)

where $g(\cdot)=\psi(\cdot)$ if $\tilde{k}^{0-\infty}_{XYZ,i}=k$ and
$g(\cdot)=\log(\cdot)$ if $\tilde{k}^{0-\infty}_{XYZ,i}>k$.

While the modifications that define our estimator MS0-∞ might appear minor at
first, they indeed address the above explained problems of the MS and ZMADG
estimators: First, our estimator by construction restricts all nearest
neighbours of a point to the cluster of that point. Hence, our estimator does
not assume local constancy across different clusters. Second, our estimator is
invariant under a common scaling of all variables. Third, there seem to be
fewer cases than for the MS estimator in which our estimator has local zero
estimates: For example, in the case discussed for the third problem of the MS
estimator. A discussion of all cases in which MS and MS0-∞ have local zero
estimates is, however, out of scope. Thus, an empirical evaluation of the bias
towards zero is called for. Fourth, unlike the ZMADG estimator, our estimator
is not the sum of up to $8$ entropy terms but retains the same general form as
the MS estimator. Thus, our estimator is not expected to incur increased
variance, which is another hypothesis subject to empirical evaluation.

### Theoretical guarantees

We provide theoretical guarantees of our estimator under the assumptions that
(1) there are at most finitely many clusters as defined by the non-numeric
components of $XYZ$ and (2) all numeric components of $XYZ$ have a finite
range. We are confident that the theoretical guarantees also hold without the
second assumption and that, to prove them, only mild adaptions of the
corresponding proofs in Mesner and Shalizi (2019) are needed. However, we
consider such an adaption to be out of scope here and leave it to future work.

The presented theoretical guarantees concern the $k$-NN limit
$\lim_{n\to\infty}$ with $k\to\infty$ and $\tfrac{k}{n}\to 0$,333Note that our
above heuristic choice of $k$ does not lead to $k/n\to 0$ as $n\to\infty$.
However, that above choice should be considered as heuristic for how to choose
$k$ for finite $n$, whereas the convergence results for $n\to\infty$ require
$k/n\to 0$. where $n$ is the sample size, and are based on the following
Lemma.

######  Lemma 1.

Let $q$ be a positive integer, and let $X^{\prime}Y^{\prime}Z^{\prime}$ be
obtained by applying a common non-constant affine function
$h:\mathbb{R}\to\mathbb{R}$ to all numeric components of $XYZ$ such that the
ranges of all numeric components of $X^{\prime}Y^{\prime}Z^{\prime}$ are
contained within the open interval $(0,1)$.444Such a function $h$ exists due
to the second assumption. Then, the difference
$\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to the constant $0$ in $L^{q}$-norm, that is,

$\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]=0\,.$
(15)

Proof sketch. Let $i_{1},\,\ldots,i_{q}$ be arbitrary integers within
$[[1,n]]$. It suffices to show that $\mathbb{E}[\chi]$ with
$\chi=\prod_{\alpha\,=\,1}^{q}\,|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})|$
converges to $0$. To this end, we first note that
$|\hat{\xi}^{0-\infty}_{i}(X;Y|Z)|$ and
$|\hat{\xi}^{(MS)}_{i}(X^{\prime};Y^{\prime}|Z^{\prime})|$ are bounded by
$2\log(n)$ and, thus, $|\chi|\leq 4^{q}\log(n)^{q}$. Next, consider the $i$-th
sample point $w_{i}=(x_{i},y_{i},z_{i})$. This point belongs to a certain
cluster $C(w_{i})$ as determined by the non-numeric part $w_{i,d}$ (if $w_{i}$
has only numeric components, then that cluster is the entire space). The
probability $\mathbb{P}(C(w_{i}))$ that an arbitrary point lies in $C(w_{i})$
is non-zero because else $w_{i}$ would not have been in $C(w_{i})$. Then,
since $k/n\to 0$, the probability that there are at least $k+1$ many points in
$C(w_{i})$ exponentially converges to one as $n\to\infty$ according to the
Chernoff Bound. Thus, using the union bound, also the probability
$\mathbb{P}(A)$ of the event $A$ that there are at least $k+1$ many points in
all of the (not necessarily distinct) $p$ clusters
$C(w_{i_{1}}),\,\ldots,\,C(w_{i_{p}})$ exponentially converges to one.
Equivalently, $\mathbb{P}(A^{c})$ is exponentially suppressed as $n\to\infty$.
Finally, one can show that
$\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)=\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$
conditioned on $A$, and thus $\chi=0$ conditioned on $A$. Thus
$\mathbb{E}[\chi]=\mathbb{E}[\chi~{}|~{}A^{c}]\cdot\mathbb{P}(A^{c})\leq
4^{q}\log(n)^{q}\cdot\mathbb{P}(A^{c})$, which implies
$\lim_{n\to\infty}\mathbb{E}[\chi]=0$ due to the exponential suppression of
$\mathbb{P}(A^{c})$. $\square$

Writing the difference $\hat{I}^{0-\infty}(X;Y|Z)-I(X;Y|Z)$ as
$[\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})]+[\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})-I(X;Y|Z)]$
and using $I(X;Y|Z)=I(X^{\prime};Y^{\prime}|Z^{\prime})$, Lemma 1 transfers
the convergence results of the MS estimator to our estimator. Specifically, we
get the following.

######  Theorem 2.

Our CMI estimator $\hat{I}^{0-\infty}(X;Y|Z)$ is $L^{1}$-consistent in the
$k$-NN limit, that is

$\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-I(X;Y|Z)\right|\right]=0\,.$
(16)

######  Theorem 3.

Assume that, in addition to the requirements of the $k$-NN limit,
$\tfrac{\left[k\cdot\ln(n)\right]^{2}}{n}\to 0$ as $n\to\infty$. Then, our CMI
estimator is $L^{2}$-consistent, that is,

$\lim_{n\to\infty}\,Var\left[\hat{I}^{0-\infty}(X;Y|Z)\right]=0\,.$ (17)

In particular, our estimator is asymptotically unbiased and converges in
probability to the true CMI. We provide all proofs in the SM Sec. B.

## Numerical evaluation of the CMI estimators

Here, we empirically study and compare bias and variance of the MS, ZMADG and
our MS0-∞ estimator.

Experimental setup We consider four models, partly taken from Mesner and
Shalizi (2019) and Zan et al. (2022) for reproducibility. Since Mesner and
Shalizi (2019) and Zan et al. (2022) do not mention any transformations of the
continuous variables, we do not apply any preprocessing. We evaluate the mean
and variances of the estimates on $100$ realizations qualitatively using
violin plots that capture the estimates’ mean and range, and quantitatively
using statistical tests. Here, we present and study one model. Further results
and an evaluation of computational runtimes are presented in Sec. E and F of
the SM.

“Independent $Z$” ((Mesner and Shalizi 2019)): $X$ is discrete uniform
$X\sim\mathcal{U}(\\{0,\ldots,c\\})$ with $1\leq c\in\mathbb{N}$,
$Y\sim\mathcal{U}\left([X,X+2]\right)$ is continuous uniform, and
$Z=(Z_{1},\ldots,Z_{d})$ is discrete with $Z_{i}\sim Ber\left(0.5\right)$ for
$1\leq i\leq d$. The ground truth is
$I(X;Y|Z)=\ln{c}-\frac{c-1}{c}\cdot\ln{2}$. In our experiments, we set the
$c=5$, vary the sample size $n\in\\{300,600,1000,2000\\}$ and $d\in\\{1,3\\}$.
We vary the $k$ parameter using $k_{c}\in\\{0.01,0.1,0.2,0.3\\}$: for MS,
$k=k_{c}\cdot n$; for MS0-∞, $k=k_{c}\cdot n_{cl,\min}$; for ZMADG,
$k=n_{subset}\cdot k_{c}$.

Figure 1: Distribution of the CMI estimates for the "Independent $Z$" model.
Each row shows the results for different $k_{c}$. The ground truth is shown as
the dashed line.

Results: In the presented violin plots for the "Independent $Z$" model (Figure
3), we observe that for both dimensionalities $d=1$ or $d=3$, all estimators
perform comparably well only for $k_{c}=0.01$. For $k_{c}>0.01$, the MS
estimates show a bias towards zero, in contrast to the MS0-∞ estimates. The
ZMADG estimator performs well, but, as anticipated, it has the highest
variance among the estimators, especially for small $k_{c}$ and small $n$, and
we observe marginal bias for $k_{c}>0.1$. Our approach also suffers from
slight bias, e.g., for $d=1$ and $k_{c}=0.1$, and $d=3$ and smaller $n$ or
$k_{c}\neq 0.1$. For the other three models, we briefly summarize: MS has good
performance except for the case when $X,Y$ are Gaussian and confounded by $Z$,
while ZMADG performs well except when $X,Y$ are uniformly distributed and
confounded by $Z$. ZMADG generally has higher variance than MS0-∞ and MS,
especially for smaller $n$ and higher dimensionality. MS0-∞ has stable
performance across all models but suffers from slight bias for large $k_{c}$
and small $n$. In summary, we confirm the above expectations and observe that
MS0-∞ draws a compromise between the strengths and weaknesses of MS and ZMADG:
It alleviates bias towards zero at the expense of higher bias in cases with
many discrete dimensions and small sample size. While not our primary focus,
our MS0-∞ estimator can also handle mixture variables, i.e., variables that
contain partially continuous and categorical samples, more robustly across
$k_{c}$ values compared to the MS estimator, as shown by preliminary results
in SM Sec. D.

## Evaluation of Conditional Independence Tests

Figure 2: False positive rate (FPR, ideally under 0.05, log-scale) and true
positive rate (TPR, higher is better, with 1 best) with standard error bars as
the $95\%$ confidence interval (computed as described in SM Sec. E) for the
following models with standardization (upper two rows) and scaling to $(0,1)$
(bottom two rows) for the continuous variables and coupling factor $w$ as
depicted: "Confounder" (left), "Independent Z"  (center) and "Cluster-
dependent confounder". The yellow areas indicate an advantage of the MS
estimator, and the red areas indicate an advantage of the MS0-∞ estimator.

We study whether the CIT controls the false positive rate (FPR) and retains
statistical power, measured by the true positive rate (TPR) at a fixed
significance level of $0.05$ in various mixed-type data setups. We focus on a
controlled environment to obtain evaluation metrics with low error, which are
difficult to obtain otherwise. Zan et al. (2022) evaluate CIT for the ZMADG
and MS estimators on synthetic and real-world data, computing acceptance rates
using simulated fork, chain, and collider structures to evaluate FPR and TPR.
The evaluation is done using models resembling those used in Sec. Numerical
evaluation of the CMI estimators, involving different combinations of
univariate mixed and discrete $X,Y,Z$. However, these models are not ideal for
systematical CIT evaluation in a scenario close to the real world, where
difficulties such as weak dependence occur. Moreover, the authors use rank
transformations on the continuous variables, which, as they mention, can put
the MS estimator at a disadvantage due to scaling. We apply standardization,
re-scaling to $(0,1)$, and rank transformation of the continuous variables for
a comprehensive comparison.

Experimental setup Our data-generating models are inspired by the post-
nonlinear model (Zhang and Hyvärinen 2009). We consider multiple causal
structures where $Z$ is a confounder, part of a chain structure, or
independent of $X$ and $Y$. We describe the individual data generation
functions below. In all models, the coefficients $\beta_{\cdot}$ are randomly
drawn as $\beta_{\cdot}\sim\mathcal{U}([-1,1])$. We introduce dependence
between $X$ and $Y$ using an additional noise term $\eta_{W}$ that influences
both $X$ and $Y$, where $\eta_{w}\sim\mathcal{N}(0,1)$ and the coupling factor
$w$ that defines the dependency strength is $w=0$ for independence and $w>0$
for dependence. The random variable $Z=(Z_{1},\ldots,Z_{m})$ is mixed-type
with $dim_{c}$ continuous dimensions and $dim_{d}$ discrete dimensions, and
each discrete component has $n_{c}$ categories. The noise terms
$\eta_{x},\eta_{Y}$ of $X$ and $Y$ follow $\mathcal{N}(0,1)$.

"Confounder": Here, $Z_{j}\sim Bin(n_{c}-1,0.5)$ for the discrete components
$Z_{1},\ldots,Z_{dim_{d}}$ and $Z_{j}\sim\mathcal{N}(0,1)$ for the continuous
components $Z_{dim_{d}+1},\ldots,Z_{m}$, and $V\in\\{X,Y\\}$ is calculated as
follows:

$V=\\!\sum_{j\in{1,...,dim_{d}}}\beta_{i}\cdot
l^{-1}(Z_{j})+\\!\sum_{j\in{dim_{c}+1,...,m}}\\!\beta_{j}\cdot
Z_{j}+\eta_{v}+w\cdot\eta_{w}$ (18)

Here, $l^{-1}(x)=\frac{e^{x}}{1+e^{x}}$ is the inverse logit function.

"Independent $Z$": Here, $Z$ has only discrete components
$Z_{1},\ldots,Z_{dim_{d}=m}$ where $Z_{j}\sim Bin(n_{c}-1,0.5)$. $X$ and $Y$
are continous univariate, and are computed as follows:

$X=\eta_{x}\\!+w\cdot\eta_{w},\quad Y=\eta_{y}\\!+w\cdot\eta_{w}.$ (19)

"Cluster-dependent confounder": Here, $Z$ is discrete univariate $Z\sim
Bin(n_{c}-1,0.5)$. $X$ and $Y$ are continous univariate. For $Z=0$, $X$ and
$Y$ are generated according to the "Confounder"-model with coupling factor
$w>0$, while for $Z\neq 0$ the same model is used with $w=0$. Thus, the
variables are only dependent for the cluster formed by $Z=0$.

"Chain": $X$ and $Y$ are continuous, and $Z$ is discrete, and all variables
are univariate. The model is defined as follows:

$\displaystyle X=\eta_{x}+w\cdot\eta_{w},$ (20) $\displaystyle
Z=\left[\sigma_{\sim}(\beta_{x}\cdot
X,(n_{c}-1))+\eta_{z}\right]\mod(n_{c}-1),$ $\displaystyle Y=\beta_{y}\cdot
l^{-1}(Z)+\eta_{y}+w\cdot\eta_{w}$

Here, $\sigma_{\sim}(\beta_{x}\cdot X,n_{c}-1)$ denotes sampling from the
multinomial distribution with $n_{c}-1$ categories where the $a$-th category
(starting the count at 0) has probability $\tfrac{e^{a\cdot
x}}{\sum_{a^{\prime}=0}^{n_{c}-2}e^{a^{\prime}\cdot x}}$. The noise $\eta_{z}$
follows $\eta_{z}\sim Bin(2,0.7)$.

Results We present results for a sample size of $1000$, with coupling factor
for all models $w=0.5$ except for "Cluster-dependent confounder" where
$w=0.75$ for $Z=0$. The number of classes is $n_{c}=3,\,4$ for the
"Confounder," "Independent $Z$," and "Chain" models, and $n_{c}=2,3,4$ for
"Cluster-dependent confounder". We evaluate
$k_{c}=0.01,\,0.1,\,0.2,\,0.3,\,0.5$ to set $k$ as in Sec. Numerical
evaluation of the CMI estimators. For "Confounder" and "Independent Z," we
vary $dim_{d}=1,2$. We generate p-values with $300$ permuted surrogates using
$k_{perm}=5$ and repeat each experiment $100$ times. Here, we present results
with standardization and scaling to $(0,1)$, and postpone the results with
rank transformation to SM Sec.E. All code and experimental results for the CMI
and CIT evaluation will be made public upon acceptance (see Sec. G in SM).

For the "Confounder"-model, we show results where the confounder $Z$ has one
continuous and one discrete dimension (further results in SM. Sec. E). For the
standardized case, we found that MS and MS0-∞ perform best (Fig. 2, upper
left). The scaling-related problem of MS becomes apparent when variables are
scaled to $(0,1)$, which leads to a decrease in TPR as $k_{c}$ increases (Fig.
2, bottom left). We observe another effect of this problem for $n_{c}=4$ and
$k_{c}=0.5$: concurrently with the low TPR, the FPR of MS suddenly drops to 0
due to the observed and permutation statistics both being equal to 0, which
results in a p-value equal to or close to 1. Our MS0-∞ CIT performs robustly
even in this case. ZMADG typically gives satisfying results yet has lower TPR
than MS and MS0-∞ and does not control FPR for $k_{c}=0.01$. For the
"Independent $Z$"-model (Fig. 2, center), ZMADG has either low TPR or high
FPR. MS and MS0-∞ perform better for both preprocessing types, yet we
sometimes observe an elevated FPR. The scaling-related issues of MS persist
for this model with scaling to $(0,1)$ as well, while MS0-∞ again performs
consistently well (except for $k_{c}=0.01$, where FPR is not controlled).
Notably, we observe that the FPR of MS0-∞ drops to $0$ for $k_{c}=0.2$ and
$k_{c}=0.5$. However, the corresponding increase in TPR indicates that these
drops are not stemming from bias towards $0$. For all previous models, MS and
MS0-∞ show similar performance for the standardized and rank transformation
case, even as $k_{c}$ increases, since, in these models, there is dependence
in each cluster. The results for the "Cluster-dependent confounder"-model
(Fig. 2, right) show how MS performance suffers if data distributions differ
between clusters. In this case, MS0-∞ identifies dependence more accurately
than MS while controlling the FPR (with some exceptions) irrespective of the
preprocessing method, especially for $n_{c}=3$. ZMADG again suffers from
either high FPR or low TPR. Summarizing, we observe that among the three
estimators, MS and MS0-∞ have superior performance. Nonetheless, MS0-∞ seems
to be the most robust estimator across the different models and preprocessing
types. However, the hyperparameter $k_{c}$ has a crucial influence: In
general, our approach benefits from higher $k_{c}$ (and thus higher $k$),
possibly due to a reduction in variance, while for MS, the opposite holds. We
provide recommendations on setting $k_{c}$ in the SM.

## Discussion and Conclusion

Understanding the performance of CIT on heterogeneous data is pivotal for
causal discovery and relevant across many fields, such as medicine or Earth
Sciences. In this work, we evaluated the $k$-nearest neighbor CMI estimators
of Mesner and Shalizi (2019) (MS) and Zan et al. (2022) (ZMADG) for mixed-type
data and discussed their challenges: ZMADG suffers from high variance while MS
treats categorical variables as numeric via one-hot encoding, leading to the
conceptual problem of mixing categories. We proposed a modification of the MS
estimator that mitigates the latter issue. We compared the bias and variance
of the estimators and the corresponding CIT performance on synthetic models
replicating realistic settings occurring in causal discovery.

As anticipated, the ZMADG estimator has a low bias but high variance,
particularly for small $k$ and $n$. The MS estimator has low variance but
suffers from bias towards zero for larger $k$ and larger numbers of discrete
variables. Our MS0-∞ reduces bias compared to MS and variance compared to
ZMADG. Surprisingly, the ZMADG CIT test does not perform consistently,
obtaining satisfying results only for the "Confounder" model. For other data
models, it does not control false positives, and we suspect these problems
occur in case of weaker dependencies. MS and MS0-∞ perform robustly, with MS
slightly outperforming when dependence holds in all clusters, and the
dimensionality is higher. However, their performance remains comparable
considering the individual optimal $k_{c}$ value. When the continuous
variables are scaled to $(0,1)$, the scaling-related problems of MS lead to
underperformance for higher $k_{c}$, contrary to our estimator. Our method has
superior performance for the case when data distributions differ between the
clusters. Thus, both from a theoretical and an empirical perspective, we
recommend our approach as the most robust estimator for mixed-type data
scenarios. While our estimator suffers from the curse of dimensionality less
than the MS estimator, the curse of dimensionality still applies to our
estimator. Hence, sufficient samples per cluster are essential for reliable
outcomes.

Lastly, while our analysis covered a range of scenarios, an evaluation of
real-world data and causal discovery is beyond the scope of this paper and is
left for future work.

## Acknowledgements

This work was partly funded by the European Union’s Horizon 2020 research and
innovation programme (project XAIDA, Grant No. 101003469) and by the European
Research Council (ERC) (project Causal Earth, Grant No. 948112). We thank Tom
Hochsprung for his valuable comments.

## References

  * Bonferroni (1935) Bonferroni, C. E. 1935.  Il calcolo delle assicurazioni su gruppi di teste.
  * Frenzel and Pompe (2007) Frenzel, S.; and Pompe, B. 2007.  Partial mutual information for coupling analysis of multivariate time series.  _Physical review letters_ , 99 20: 204101.
  * Gao et al. (2017) Gao, W.; Kannan, S.; Oh, S.; and Viswanath, P. 2017.  Estimating Mutual Information for Discrete-Continuous Mixtures.  _CoRR_ , abs/1709.06212.
  * Gray (2011) Gray, R. M. 2011.  _Entropy and information theory_.  Springer Science & Business Media.
  * Kozachenko and Leonenko (1987) Kozachenko, L. F.; and Leonenko, N. N. 1987.  Sample estimate of entropy of a random vector.  _Problems of Information Transmission_ , (23): 95–101.
  * Kraskov, Stögbauer, and Grassberger (2004) Kraskov, A.; Stögbauer, H.; and Grassberger, P. 2004.  Estimating mutual information.  _Physical review. E, Statistical, nonlinear, and soft matter physics_ , 69 6 Pt 2: 066138.
  * Levene (1960) Levene, H. 1960.  Robust test for equality of variances.  _Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling (I. Olkin ed.)_ , 278–292.
  * Mesner and Shalizi (2019) Mesner, O. C.; and Shalizi, C. R. 2019.  Conditional Mutual Information Estimation for Mixed Discrete and Continuous Variables with Nearest Neighbors.
  * Runge (2018) Runge, J. 2018.  Conditional independence testing based on a nearest-neighbor estimator of conditional mutual information.  In Storkey, A.; and Perez-Cruz, F., eds., _International Conference on Artificial Intelligence and Statistics_ , volume 84 of _Proceedings of Machine Learning Research_ , 938–947. PMLR.
  * Spirtes et al. (2000) Spirtes, P.; Glymour, C. N.; Scheines, R.; and Heckerman, D. 2000.  _Causation, prediction, and search_.  MIT press.
  * Wasserman (2010) Wasserman, L. 2010.  _All of statistics : a concise course in statistical inference_.  New York: Springer.  ISBN 9781441923226 1441923225.
  * Wilcoxon (1945) Wilcoxon, F. 1945.  Individual Comparisons by Ranking Methods.  _Biometrics_ , 1: 196–202.
  * Zan et al. (2022) Zan, L.; Meynaoui, A.; Assaad, C. K.; Devijver, E.; and Gaussier, É. 2022.  A Conditional Mutual Information Estimator for Mixed Data and an Associated Conditional Independence Test.  _Entropy_ , 24.
  * Zhang and Hyvärinen (2009) Zhang, K.; and Hyvärinen, A. 2009.  On the Identifiability of the Post-Nonlinear Causal Model.  In McAllester, D., ed., _Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence_ , UAI ’09, 647–655. Arlington, Virginia, USA: AUAI Press.  ISBN 9780974903958.

Non-parametric Conditional Independence Testing for Mixed Continuous-
Categorical Variables: A Novel Method and Numerical Evaluation
SUPPLEMENTARY MATERIAL

## Nearest-Neighbor Permutation Test

Algorithm 1 (Wasserman 2010) describes the procedure for generating the
p-value from a set of $B$ permutations $\Pi=\\{\pi_{j}|j=1,...,B\\}$ of $n$
elements. Each permutation is generated using the nearest-neighbor scheme
described in Sec. 2 of the main paper. Given a test statistic, in our case,
the conditional mutual information (CMI) estimated using one of the estimators
(MS, MS0-∞, or ZMADG), denoted as $T_{CMI}$, we obtain the p-value of the
conditional independence test (CIT) of the respective estimator as follows:

Input: samples $\\{(x_{i},y_{i},z_{i})|i=1,...,n\\}$, permutation set
$\Pi=\\{\pi_{j}|j=1,...,B\\}$, test statistic $T_{CMI}$, nominal level
$\alpha$

Output: p-value

1 For each permutation $\pi_{j}\in\Pi,j=1,...,B$: compute the statistic
$T_{CMI}^{\pi_{j}}$ on the permuted samples

Compare the test statistic on the observed samples, $T_{CMI}$ with the
permuted test statistics and calculate the p-value as
$p=\frac{1}{B}\sum_{\pi_{j}\in\Pi}\mathbf{1}\\{T_{CMI}^{\pi_{j}}\geq
T_{CMI}\\}$ (21)

Algorithm 1 Local Permutation Test

## Proofs

In this part of the Supplementary Material, we formally prove the theoretical
claims made in Sec. 3 of the main paper. Our proofs rely on the following two
assumptions:

######  Assumption 4.

There are at most finitely many clusters as defined by the non-numeric
components of $XYZ$

######  Assumption 5.

All numeric components of $XYZ$ have a finite range.

Note that also (Mesner and Shalizi 2019) implicitly make Assumption 4 by
assuming finite-dimensional random vectors in combination with one-hot
encoding the non-numeric components. As we mentioned in the main text, we are
confident that the theoretical guarantees of our estimator also hold without
Assumption 5 and that, to prove them in this case, only mild adaptions of the
corresponding proofs in Mesner and Shalizi (2019) are needed. However, we
leave such an adaption of the proofs to future work and here do adopt
Assumption 5.

Proof of Lemma 1.

We start by using the triangle inequality to get

$\displaystyle\begin{aligned}
&\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]\\\
=&\lim_{n\to\infty}\,\mathbb{E}\left[\left|\frac{1}{n}\sum_{i\,=\,1}^{n}\hat{\xi}^{0-\infty}_{i}(X;Y|Z)-\frac{1}{n}\sum_{i\,=\,1}^{n}\hat{\xi}^{MS}_{i}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]\\\
\leq&\lim_{n\to\infty}\,\frac{1}{n^{q}}\,\mathbb{E}\left[\left(\sum_{i\,=\,1}^{n}\left|\hat{\xi}^{0-\infty}_{i}(X;Y|Z)-\hat{\xi}^{MS}_{i}(X^{\prime};Y^{\prime}|Z^{\prime})\right|\right)^{q}\,\right]\\\
=&\lim_{n\to\infty}\,\frac{1}{n^{q}}\sum_{i_{1}\,=\,1}^{n}\sum_{i_{2}\,=\,1}^{n}\dots\\\
&\sum_{i_{q}\,=\,1}^{n}\,\mathbb{E}\left[\prod_{\alpha\,=\,1}^{q}\left|\
\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|\right]\,.\end{aligned}$
(22)

Next, let $N_{n}(C_{i})$ be the number of points other than the $i$-th point
that are in the cluster $C_{i}$ of the $i$-th point. This random variable
$N_{n}(C_{i})$ follows the distribution $Bin(n-1,p_{i})$, where
$p_{i}=\mathbb{P}(C_{i})$ is the probability that an arbitrary point belongs
to the cluster $C_{i}$. We have $p_{i}>0$ because else the $i$-th point would
not have been in the cluster $C_{i}$. The Chernoff bound for the lower tail of
the Binomial distribution then gives

$\mathbb{P}(N_{n}(C_{i})\leq k-1)\leq e^{-\tfrac{(n-1)\cdot
p_{i}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot p_{i}}\right]^{2}}\,.$ (23)

Now let $A_{i_{\alpha}}$ be the event that $N_{n}(C_{i_{\alpha}})\leq k-1$ and
let $A_{i_{1},\,\ldots,i_{q}}=\cup_{\alpha\,=\,1}^{q}A_{i_{\alpha}}$. The
probability of $A_{i_{1},\,\ldots,i_{q}}$ is upper bounded according to

$\displaystyle\mathbb{P}(A_{i_{1},\,\ldots,i_{q}})$
$\displaystyle\leq\sum_{\alpha\,=\,1}^{q}\mathbb{P}(N_{n}(C_{i_{\alpha}})\leq
k-1)$ (24) $\displaystyle\leq\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{i_{\alpha}}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{i_{\alpha}}}\right]^{2}}\,.$

Next, we use the law of total expectation to condition the expectation value
on the right-hand-side of the last line of ineq. (22) on the events
$A_{i_{1},\,\ldots,i_{q}}$ and $A_{i_{1},\,\ldots,i_{q}}^{c}$, which gives

$\displaystyle\mathbb{E}\left[\prod_{\alpha\,=\,1}^{q}\left|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|\right]$
(25)
$\displaystyle=\,\,\underbrace{\mathbb{E}\left[\prod_{\alpha\,=\,1}^{q}\left|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|~{}\middle|~{}A_{i_{1},\,\ldots,i_{q}}\right]}_{\geq
0}\cdot\underbrace{\mathbb{P}(A_{i_{1},\,\ldots,i_{q}})}_{\leq\,\sum_{\alpha\,=\,1}^{q}\exp\left\\{-\tfrac{(n-1)\cdot
p_{i_{\alpha}}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{i_{\alpha}}}\right]^{2}\right\\}}$
$\displaystyle\,\,+\mathbb{E}\Bigg{[}\,\,\underbrace{\prod_{\alpha\,=\,1}^{q}\left|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|~{}\Bigg{|}~{}A_{i_{1},\,\ldots,i_{q}}^{c}}_{=\,0}\,\,\Bigg{]}\cdot\underbrace{\mathbb{P}(A_{i_{1},\,\ldots,i_{q}}^{c})}_{\leq\,1}$
$\displaystyle\leq\mathbb{E}\left[\prod_{\alpha\,=\,1}^{q}\left|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|~{}\middle|~{}A_{i_{1},\,\ldots,i_{q}}\right]\cdot\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{i_{\alpha}}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{i_{\alpha}}}\right]^{2}}$

Here, we have used
$\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)=\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$
conditioned on the event $A_{i_{1},\,\ldots,i_{q}}^{c}$. This equality holds
for the following reason: Given $A_{i_{1},\,\ldots,i_{q}}^{c}$, for all
$1\leq\alpha\leq q$ at least $k$ points other than the $i_{\alpha}$-point are
in the cluster $C_{i_{\alpha}}$ of $i_{\alpha}$. Now fix some $\alpha$ with
$1\leq\alpha\leq q$ and let $v_{0},v_{1},\ldots,v_{m}$ be the points in
cluster $C_{i_{\alpha}}$, ordered such that $v_{0}$ is the $i_{\alpha}$-th
point and that $\|v_{a}-v_{0}\|_{L^{\infty}}\leq\|v_{b}-v_{0}\|_{L^{\infty}}$
if $a<b$. Then, $m\geq k$ and the distance $\rho_{i_{\alpha}}$ used by the
estimate $\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)$, here denoted as
$\rho_{i_{\alpha}}^{0-\infty}$, equals $\|v_{k}-v_{0}\|_{L^{\infty}}<\infty$.
Let $h(w)=\lambda\cdot w+c$ with $\lambda\neq 0$ be the non-constant affine
function that transforms the numeric components of $XYZ$ to
$X^{\prime}Y^{\prime}Z^{\prime}$, and for all $0\leq a\leq m$ let
$v_{a}^{\prime}$ denote the transformed version of $v_{a}$. Since the same
transformation $h$ is applied to all numeric components of $XYZ$, we get that
$\|v^{\prime}_{a}-v^{\prime}_{0}\|_{L^{\infty}}\leq\|v^{\prime}_{b}-v^{\prime}_{0}\|_{L^{\infty}}$
if $a<b$. Moreover, since the ranges of numeric components of
$X^{\prime}Y^{\prime}Z^{\prime}$ are contained within $(0,1)$, we get that
$\|v^{\prime}_{m}-v^{\prime}_{0}\|_{L^{\infty}}<1$. Consequently, for the
purpose of $\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$ the
$k$-nearest neighbours of $v^{\prime}_{0}$ are
$v^{\prime}_{1},\,\ldots,v^{\prime}_{k}$ and the distance $\rho_{i_{\alpha}}$
used by the estimate
$\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$, here denoted
as $\rho_{i_{\alpha}}^{MS}$, equals
$\|v^{\prime}_{k}-v_{0}\|_{L^{\infty}}=\lambda\cdot\rho_{i_{\alpha}}^{0-\infty}<1$.
We thus see that also
$\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$ uses only
points within the cluster $C_{i_{\alpha}}$ of the $i_{\alpha}$-th point. Let
$W$ be a wildcard for $XYZ$, $XZ$, $YZ$ and $Z$ and consider the count
$\tilde{k}^{0-\infty}_{W,i_{\alpha}}$ used by the estimate
$\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)$ as defined in eq. (14) in the main
text. By definition, see eqs. (12) and (13) in the main text, this count is
the number of points other than the $i_{\alpha}$-th point $v_{0}$ that in
$W$-space have a $0-\infty$ “distance” of at most
$\rho^{0-\infty}_{i_{\alpha}}$ to $v_{0}$. Since
$\rho^{MS}_{i_{\alpha}}=\lambda\cdot\rho^{0-\infty}_{i_{\alpha}}<1$ and the
$L^{\infty}$-distance in $W^{\prime}$-space (where, for example,
$W^{\prime}=X^{\prime}Y^{\prime}Z^{\prime}$ if $W=XYZ$) is $\lambda$ times the
$L^{\infty}$-distance in $W$-space, the corresponding count
$\tilde{k}_{W^{\prime},i_{\alpha}}$ used by
$\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$, see eqs. (9)
and (10) in the main text, equals $\tilde{k}^{0-\infty}_{W,i_{\alpha}}$. From
the equality
$\tilde{k}^{0-\infty}_{W,i_{\alpha}}=\tilde{k}_{W^{\prime},i_{\alpha}}$ we
conclude that, as claimed,
$\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)=\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})$
conditioned on the event $A_{i_{1},\,\ldots,i_{q}}^{c}$.

To upper bound the remaining conditional expectation, we use the triangle
inequality and the fact that $\ln(a)\geq\psi(a)\geq 0$ for $a\geq 0$ to get

$\displaystyle\left|\hat{\xi}^{0-\infty}_{i}(X;Y|Z)\right|\leq 2\ln(n)\hskip
5.0pt\text{and}$ (26) $\displaystyle\hskip
5.0pt\left|\hat{\xi}^{MS}_{i}(X^{\prime};Y^{\prime}|Z^{\prime})\right|\leq
2\ln(n)\,,$

which implies

$\left|\hat{\xi}^{0-\infty}_{i}(X;Y|Z)-\hat{\xi}^{MS}_{i}(X^{\prime};Y^{\prime}|Z^{\prime})\right|\leq
4\ln(n)\,.$ (27)

We thus find

$\displaystyle\begin{aligned}
\mathbb{E}\left[\prod_{\alpha\,=\,1}^{q}\left|\hat{\xi}^{0-\infty}_{i_{\alpha}}(X;Y|Z)-\hat{\xi}^{MS}_{i_{\alpha}}(X^{\prime};Y^{\prime}|Z^{\prime})\right|~{}\middle|~{}A_{i_{1},\,\ldots,i_{q}}\right]\\\
\leq 4^{q}\ln(n)^{q}\,.\end{aligned}$ (28)

By combining the above results, we get

$\displaystyle\begin{aligned}
&\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]\\\
\leq&\lim_{n\to\infty}\frac{4^{q}\cdot\ln(n)^{q}}{n^{q}}\cdot\\\
&\,\left(\sum_{i_{1}\,=\,1}^{n}\sum_{i_{2}\,=\,1}^{n}\dots\sum_{i_{q}\,=\,1}^{n}\right)\,\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{i_{\alpha}}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{i_{\alpha}}}\right]^{2}}\,.\end{aligned}$ (29)

Let $C^{(1)},\,\ldots,C^{(m)}$ be the list of all clusters with positive
probability, where $m<\infty$ because according to Assumption 4 there are at
most finitely many clusters. Then
$p_{min}=\min_{\gamma\in[[1,m]]}\mathbb{P}(C^{(\gamma)})$, which is
independent of $k$ and $n$, exists and $p_{min}>0$. Noting that
$p_{i_{\alpha}}=\mathbb{P}(C_{i_{\alpha}})>0$ for all $i_{\alpha}$ because
else the $i_{\alpha}$-th point would not have been in the cluster
$C_{i_{\alpha}}$, we get that $1\geq p_{i_{\alpha}}\geq p_{min}>0$ for all
$i_{\alpha}$ and hence

$\frac{k-1}{(n-1)\cdot p_{i_{\alpha}}}\leq\frac{k-1}{(n-1)\cdot p_{min}}$ (30)

for all $i_{\alpha}$. Since $\tfrac{k}{n}\to 0$ and $k\to\infty$ in the $k$-nn
limit and since $p_{i_{\alpha}}$ is independent of $k$ and $n$, we find that
$\tfrac{k-1}{(n-1)\cdot p_{i_{\alpha}}}\to 0^{+}$ for all $i_{\alpha}$ and
$\tfrac{k-1}{(n-1)\cdot p_{min}}\to 0^{+}$. Thus, there is a positive integer
$n_{0}$ such that for all $n\geq n_{0}$ and for all $i_{\alpha}$ the bound

$\left[1-\frac{k-1}{(n-1)\cdot
p_{i_{\alpha}}}\right]^{2}\geq\left[1-\frac{k-1}{(n-1)\cdot
p_{min}}\right]^{2}$ (31)

holds. We then get

$\displaystyle\begin{aligned}
&\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]\\\
\leq&\lim_{n\to\infty}\frac{4^{q}\cdot\ln(n)^{q}}{n^{q}}\cdot\\\
&\,\left(\sum_{i_{1}\,=\,1}^{n}\sum_{i_{2}\,=\,1}^{n}\dots\sum_{i_{q}\,=\,1}^{n}\right)\,\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{min}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{min}}\right]^{2}}\,.\end{aligned}$ (32)

We can now pull the sum of exponentials (that is, the sum over $\alpha$) out
of the product-sum (that is, the sums over $i_{1},\,\ldots,i_{q}$) and get

$\displaystyle\begin{aligned}
&\lim_{n\to\infty}\,\mathbb{E}\left[\left|\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})\right|^{q}\right]\\\
\leq&\\!\begin{aligned}
&\lim_{n\to\infty}\frac{4^{q}\cdot\ln(n)^{q}}{n^{q}}\cdot\\\
&\,\left(\sum_{i_{1}\,=\,1}^{n}\sum_{i_{2}\,=\,1}^{n}\dots\sum_{i_{q}\,=\,1}^{n}\right)\,\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{min}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{min}}\right]^{2}}\end{aligned}\\\ =&\\!\begin{aligned}
&\lim_{n\to\infty}\,\,\frac{4^{q}\cdot\ln(n)^{q}}{n^{q}}\cdot\\\
&\underbrace{\sum_{\alpha\,=\,1}^{q}e^{-\tfrac{(n-1)\cdot
p_{min}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{min}}\right]^{2}}}_{=\,q\cdot\exp\left\\{-\tfrac{(n-1)\cdot
p_{min}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot
p_{min}}\right]^{2}\right\\}}\,\underbrace{\left(\sum_{i_{1}\,=\,1}^{n}\sum_{i_{2}\,=\,1}^{n}\dots\sum_{i_{q}\,=\,1}^{n}\right)}_{n^{q}}\end{aligned}\\\
=&\lim_{n\to\infty}\,4^{q}\cdot\ln(n)^{q}\cdot q\cdot e^{-\tfrac{(n-1)\cdot
p_{min}}{2}\cdot\left[1-\tfrac{k-1}{(n-1)\cdot p_{min}}\right]^{2}}\\\
=&\,0\end{aligned}$ (33)

The last equality follows because the argument of the exponential goes to
$-\infty$ in the $k$-nn limit. ∎

Proof of Theorem 2. The sequence
$X_{n}=\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $X=0$ in $L^{1}$ according to Lemma 1 with $q=1$, and the
sequence
$Y_{n}=\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})-I(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $Y=0$ in $L^{1}$ according to the proof of Theorem 3.1 in (Mesner
and Shalizi 2019) 555That theorem itself only claims asymptotic unbiasedness,
which is strictly weaker than $L^{1}$-convergence, but the proof actually
shows $L^{1}$-convergence. Therefore, the sequence
$Z_{n}=X_{n}+Y_{n}=\hat{I}^{0-\infty}(X;Y|Z)-I(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $X+Y=0$ in $L^{1}$. We conclude the proof by noting the equality
$I(X^{\prime};Y^{\prime}|Z^{\prime})=I(X;Y|Z)$, which follows because CMI is
invariant under componentwise non-constant affine transformations. ∎

Proof of Theorem 3. The sequence
$X_{n}=\hat{I}^{0-\infty}(X;Y|Z)-\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $X=0$ in $L^{2}$ according to Lemma 1 with $q=2$, and the
sequence
$Y_{n}=\hat{I}^{MS}(X^{\prime};Y^{\prime}|Z^{\prime})-I(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $Y=0$ in $L^{2}$ according to Theorem 3.2 in Mesner and Shalizi
(2019) in combination with Theorem 3.1 in Mesner and Shalizi (2019).
Therefore, the sequence
$Z_{n}=X_{n}+Y_{n}=\hat{I}^{0-\infty}(X;Y|Z)-I(X^{\prime};Y^{\prime}|Z^{\prime})$
converges to $X+Y=0$ in $L^{2}$. We conclude the proof by noting the equality
$I(X^{\prime};Y^{\prime}|Z^{\prime})=I(X;Y|Z)$, which follows because CMI is
invariant under componentwise non-constant affine transformations. ∎

## Further Results of the Numerical Evaluation of the CMI Estimators

Here, we present further results for the numerical evaluation of the CMI
estimators.

Experimental setup We keep the experimental setup as described in Sec. 4 of
the main paper and generate further synthetic datasets according to the
following models:

"Chain structure” ((Mesner and Shalizi 2019)): Here,
$X\sim\exp\left(10\right)$, $Z=(Z_{1},\ldots,Z_{d})$ is multivariate with
$Z_{1}\sim Poisson\left(X\right)$ and $Z_{i}\sim\mathcal{N}(0,1)$ for $2\leq
i\leq d$, and $Y\sim Bin\left(Z_{1},\ 0.5\right)$. The ground truth is
$I(X;Y|Z)=0$.

"Confounder with Gaussian $X$ and $Y$" (Zan et al. (2022)): This model
describes a confounder structure with normally distributed $X$ and $Y$, where
$X\sim\mathcal{N}\left(Z,1\right)$, $Y\sim\mathcal{N}\left(Z,\ 1\right)$, and
$Z\sim\mathcal{U}\left(\\{0,\ldots,m\\}\right)$. The ground truth is
$I(X,Y|Z)=0$. As in Zan et al. (2022), $m=9$.

"Confounder with uniform $X$ and $Y$": This model describes a confounder
structure with uniformly distributed $X$ and $Y$, where
$X\sim\mathcal{U}\left(0,Z\right)$, $Y\sim\mathcal{U}\left(Z,Z+1\right)$, and
$Z\sim\mathcal{U}\left(\\{0,1\\}\right)$. The ground truth is $I(X,Y|Z)=0$.

Figure 3: Distribution of the CMI estimates for the "Chain structure" model
with $d=1$ on the left and $d=3$ on the right. Each row shows the results for
different $k_{c}$. Ground truth shown as dashed line.

Results For the "Chain structure" model (Fig. 3) with $d=1$, the MS0-∞ and
ZMADG estimators perform comparably well. Notably, the ZMADG estimator has
slightly higher variance than MS0-∞ and MS, especially for smaller $n$ and
higher $d=3$. For the "Confounder with Gaussian $X$ and $Y$" model (Fig. 4,
left), we observe that MS0-∞ performs best, while the MS estimator
overestimates for $k_{c}\geq 0.1$. The ZMADG estimator again suffers from
higher variance compared to our estimator. For the "Confounder with uniform
$X$ and $Y$" (Fig. 4, right), besides having high variance, ZMADG wrongly
finds a strong conditional dependency between $X$ and $Y$. It seems specific
to the uniform distribution, as estimates are correct for the same model using
normally distributed data. For this model, both MS and our estimator perform
well.

Figure 4: Distribution of the CMI estimation for "Confounder with Gaussian $X$
and $Y$" model on the left and "Confounder with uniform $X$ and $Y$" model on
the right. Ground truth values are shown as the dashed line. Each row shows
results for a different $k_{c}$ value.

### Results of the Numerical Evaluation of the Estimators for the Mixture-type
Variable Case

As mentioned in the main paper, our study did not focus on the mixture-type
variable case. However, our MS0-∞ estimator can be used with mixture-type
variables. We present here a preliminary evaluation of the CMI estimation for
this use case. We keep the experimental setup as for the previous experiments,
but do not evaluate the ZMADG estimator, since this method was not designed
for use with mixture-type variables. We thus evaluate only the MS and MS0-∞
estimators on data generated from a model inspired by Mesner and Shalizi
(2019), defined as follows:

"Mixture" (adapted from Mesner and Shalizi (2019)): Here, $Z$ is discrete with
$Z\sim Bin(1,p=0.3)$. With probability $1-p$, $X$ and $Y$ are drawn from a
multivariate Gaussian with correlation coefficient of $0.6$, and with
probability $p$, $X\sim\mathcal{U}(\\{0,...,4\\})$ and
$Y\sim\mathcal{U}\left([X,X+2]\right)$. The ground truth is
$I(X;Y|Z)=-(1-p)\cdot\ln{(1-0.36)}\cdot
0.5+p\cdot(\ln{5}-\frac{4}{5}\cdot\ln{2})=0.472$.

Figure 5: Distribution of CMI estimates for the "Mixture" model. Each row
shows the results for different $k_{c}$. Ground truth shown as dashed line.

Results As Fig. 5 shows, our approach performs best across all $k_{c}$ values.
Both MS and MS0-∞ suffer from bias towards 0 for $k_{c}\geq 0.2$, however, our
estimator is considerably less affected than the MS estimator.

## Choice of Heuristic for $k$

In contrast to the theoretical setting, where we assume infinite samples are
available, we only have access to a finite number of samples in practice.
Thus, for large enough $k$ with fixed $n$, it is probable that some of the
$k$-th nearest neighbors in the dataset are at $\infty$ distance, i.e.,
originate from a different cluster. Because our approach does not allow
neighbors from different clusters, a heuristic was necessary to deal with this
particular case. We defined and tested three different heuristics. In this
Section, we describe the other two heuristics besides the heuristic described
in Sec. 3 of the main paper. We then motivate our choice for the heuristic in
the main paper, which we based on an empirical comparison of the MS0-∞
estimator’s performance using the different heuristics.

### Two Alternative Heuristics

We start with the presentation of the other two heuristics, namely the
"global" and "cluster-size" heuristics. Note that from now on, we will refer
to the heuristic described in Sec. 3 of the main paper as the "local"
heuristic.

#### "Global" heuristic

The "global" heuristic defines a "global" $k$ as a fraction of the number of
samples $n$. If the distance to the $k$-th nearest neighbor is $\infty$, then,
instead of enforcing $k$ nearest neighbours for all sample points, we allow
for the following adaptiveness: If the $k$-th NN of $w_{i}$ is at distance
$\infty$ from $w_{i}$ (that is, if $n_{cl}^{i}\leq k$ with $n_{cl}^{i}$ the
number of points in the cluster of $w_{i}$), then for this $i$ we replace $k$
by $k_{i}=\lfloor k_{c}\cdot n^{i}_{cl}\rfloor$. Explicitly: For all $i$ let
$k^{0-\infty}_{i}=k$ if $k+1\leq n^{i}_{cl}$ and $k^{0-\infty}_{i}=\lfloor
k_{c}\cdot n^{i}_{cl}\rfloor$ else. Thus, in effect, all considered nearest
neighbours of $w_{i}$ come from the same cluster as $w_{i}$.

#### "Cluster-size" heuristic

This heuristic still uses a "global" $k$ as a fraction of the number of
samples $n$. However, if the distance to the $k$-th nearest neighbor is
$\infty$, the "cluster-size" heuristic deals with this case by simply setting
$k_{i}=n^{i}_{cl}$, where $n^{i}_{cl}$ is the number of samples in the cluster
of point $i$, defined as previously described.

### Numerical Evaluation of the Three Heuristics

Experimental setup We run numerical experiments to compare the bias and
variance of our MS0-∞ estimator using the three different heuristics,
alongside the bias and variance of the MS and ZMADG estimators. We keep the
same experimental setting as in the previous experiments (described in the
Sec. 4 of the main paper and Sec. Further Results of the Numerical Evaluation
of the CMI Estimators of the SM) and evaluate results on the "Independent $Z$"
and "Chain structure" models.

Results The violin plots in Figures 6 and 7 show the results of the CMI
estimation using the three heuristics presented above. We observe that the
"cluster-size" heuristic suffers from bias towards zero. This is expected
because, when $k=n^{i}_{cl}$, the distance to the $k$-th nearest neighbor is
equal to the distance from point $i$ to the farthest point in its respective
cluster. Thus, for the subspaces $XY,XZ$ and $Z$, the number of counted
neighbors is equal to $n^{i}_{cl}$ with high probability, which results in an
estimation equal to or close to 0. The "global" heuristic has the highest
variance across the different heuristics, yet still has lower variance
compared to the ZMADG estimator. The "global" approach also slightly suffers
from bias for higher dimesionality, for example, for the "Chain structure"
model with $d=3$. The best bias-variance trade-off is obtained using the
"local" heuristic, where both bias and variance are smaller compared to the
other two heuristics. This motivates our choice to use this heuristic for the
experiments in the main paper.

### Discussion on the Choice of $k_{c}$

Generally, $k$-NN methods benefit from a value of $k$ that is neither too high
or too low: A small $k$ leads to lower bias at the cost of increased variance,
while larger $k$ values lead to low variance but increased bias (Kraskov,
Stögbauer, and Grassberger 2004). From the CMI estimation results and the
outcome of the conditional independence tests, we observe that a larger
$k_{c}$, e.g. $k_{c}>0.2$ is beneficial for our method, especially for the
case of weak dependency and smaller sample size. However, as the sample size
increases, e.g. $n=2000$, we observe that a smaller $k_{c}$, e.g. $k_{c}=0.2$
performs better. We thus generally recommend $k_{c}\geq 0.1$, with larger
$k_{c}$ when the number of samples is small or there are many clusters, as an
appropriate choice that reduces the bias and variance of the estimation.

Figure 6: CMI estimation results for the "Independent $Z$" model using the
three different heuristics for setting the hyperparameter $k$: the "local"
heuristic (used in the main paper, MS${}_{0}-\infty$ local), the "global"
heuristic (MS${}_{0}-\infty$ global), and the "cluster-size" heuristic
(MS${}_{0}-\infty$ cluster). The ground truth CMI values are indicated by the
dashed line. The first four rows show results for different $k_{c}$ values for
the model with $d=1$, and the last four rows show results for different
$k_{c}$ values for the model with $d=3$.

Figure 7: CMI estimation results for the "Chain structure" model using the
three different heuristics for setting the hyperparameter $k$: the "local"
heuristic (MS${}_{0}-\infty$ local), the "global" heuristic (MS${}_{0}-\infty$
global), and the "cluster-size" heuristic (MS${}_{0}-\infty$ cluster). The
ground truth CMI values are indicated by the dashed line. The first four rows
show results for different $k_{c}$ values for the model with $d=1$, and the
last four rows show results for different $k_{c}$ values for the model with
$d=3$.

## Further Results of the Numerical Evaluation of the CIT

We present further results of the CIT evaluation for the configurations
presented in Sec. 5 of the main paper. As in Sec. 5, we show plots of the FPR
and TPR for different $k_{c}$ values. In the following sections, we also
present additional plots of the true null and the permuted CMI distributions
for the different $k_{c}$ values. These distribution plots allow us to
investigate whether the FPR/TPR reflects the desired behavior of the tests:
The true null and the permuted distributions should have CMI values
distributed around 0. In contrast, in the dependent case, the CMI values
should be larger than 0, and their distribution should have minimal overlap if
the null hypothesis does not hold. While we investigated the CMI distribution
plots of all configurations, we refrain from adding all plots in the SM for
length reasons and exemplify using one set of plots.

### Computation of the Confidence Intervals for the TPR/FPR Plots

Before presenting the results, we first describe how we compute the error bars
of the CIT plots, which represent the $95\%$ confidence interval of the false
positive rate (FPR) and true positive rate (TPR).

The confidence intervals are obtained by modelling the false and true
positives as distributed according to the binomial distribution. We describe
the computation of the confidence interval for the FPR, and obtain the
confidence interval for the TPR analogously. For a given model and a set of
values of the CIT parameters, the probability to obtain a false positive in
the $n_{rep}=100$ repetitions of our experiments is $p_{FP}$ (ideally,
$p_{FP}=\alpha$). Under the assumption that the repetitions are independent,
the random variable that describes the number of false positives, $FP$, is
distributed according to the binomial distribution:

$FP\sim Bin(n_{rep},p_{FP}).$ (34)

We can estimate $p_{FP}$ as the empirical fraction of false positives that we
have obtained in our repetitions: $\hat{p}_{FP}=\frac{\\#FP}{n_{rep}}=FPR$.

We now wish to obtain a confidence interval for $\hat{p}_{FP}$, i.e., find the
lower and upper bounds $p_{L},p_{U}$ of the confidence interval such that
$P(p_{L}(FP)<\hat{p}_{FP}<p_{U}(FP))=1-\alpha$. Since $FP\sim
Bin(n_{repetitions},p_{FP})$, and supposing we have observed $FP=fp$, we
obtain $p_{L}$ and $p_{U}$ by numerically solving the following two equations:

$\displaystyle\sum_{k=fp}^{n_{repetitions}}$
$\displaystyle(p_{L})^{k}\cdot(1-p_{L})^{(n_{repetitions}-k)}$ (35)
$\displaystyle=$ $\displaystyle 1-CDF(fp,n_{repetitions},p_{L})$
$\displaystyle=$ $\displaystyle\frac{1-\alpha}{2}$

$\displaystyle\sum_{k=0}^{fp}$
$\displaystyle(p_{U})^{k}\cdot(1-p_{U})^{(n_{repetitions}-k)}$ (36)
$\displaystyle=$ $\displaystyle CDF(fp,n_{repetitions},p_{U})$
$\displaystyle=$ $\displaystyle\frac{1-\alpha}{2}$

Here $CDF$ is the cumulative distribution function of the binomial
distribution.

### "Confounder" Model

Here, we present further results for the "Confounder" model.

Figure 8: False positive rate (FPR, ideally under 0.05) and true positive rate
(TPR, higher is better, with 1 best) for the "Confounder" model where $Z$ has
$dim_{c}=1,dim_{d}=1$, with rank preprocessing for the continuous variables
and coupling factor in the dependent case $w=0.5$.

For the "Confounder"-model with $dim_{c}=1$ and $dim_{d}=1$ and sample size
$n=1000$, the CIT with rank transformation (Fig. 8) behaves similarly to the
CIT with standardization. Inspecting the distributions of the CMI values (see
Fig. 9 and 10 for $n_{c}=3$ and Fig. 11 and 12 for $n_{c}=4$), we observe that
MS and MS0-∞ estimators behave similarly for $k_{c}>0.1$, while ZMADG suffers
from slight negative bias.

Figure 9: The true null and permuted distributions of the CMI estimated with
the MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$
where $Z$ has $dim_{c}=1,dim_{d}=1$ and $n_{c}=3$ with standardization
preprocessing. The left figure of each plot pair shows the true null
distribution under $H_{0}$ as the orange line with the red dot indicating the
$95\%$ quantile, and the permuted null distribution is shown as the black
line, with the black dot indicating the $95\%$ quantile. On the right, we show
the true distribution under $H_{1}$ as the grey line, and the permuted
distributions as the black, with the black dot indicating the $95\%$ quantile.

Figure 10: The true null and permuted distributions of the CMI estimated with
the MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$
where $Z$ has $dim_{c}=1,dim_{d}=1$ and $n_{c}=3$ with scaling to $(0,1)$ and
rank preprocessing. The distributions are depicted as described in Fig. 9.

Figure 11: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model where $Z$ has
$dim_{c}=1,dim_{d}=1$ and $n_{c}=4$ with standardization preprocessing. The
distributions are depicted as described in Fig. 9.

Figure 12: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$ where
$Z$ has $dim_{c}=1,dim_{d}=1$ and $n_{c}=4$ with scaling to $(0,1)$ and rank
preprocessing. The distributions are depicted as described in Fig. 9.

Figure 13: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$ where
$Z$ has $dim_{c}=1,dim_{d}=2$ and $n_{c}=3$ with standardization
preprocessing. The distributions are depicted as described in Fig. 9.

Figure 14: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Confounder"-model with
$n=1000$ where $Z$ has $dim_{c}=1,dim_{d}=2$, with standardization, scaling to
$(0,1)$ and rank preprocessing for the continuous variables and coupling
factor in the dependent case $w=0.5$.

As the dimensionality of the discrete variable increases to $dim_{d}=2$ (Fig.
14), we observe that our approach performs better with higher $k_{c}$, e.g.,
$k_{c}\geq 0.3$. The CMI distribution plots (see Fig. 13, 15, 16 and 17)
indicate that our approach suffers from positive bias in the case of small
$k_{c}$ and high dimensionality. For $n_{c}=3$ with both standardization and
rank preprocessing, our approach performs slightly worse than MS , but
performs slightly better than ZMADG, which generally suffers from negative
bias. For $n_{c}=4$, our approach performs worse than MS and slightly worse
than ZMADG. Due to the definition of the data model, dependence holds in every
cluster. Thus, the MS estimator performs well in the standardized and rank
preprocessing cases, despite the fact that points from other clusters are
considered neighbors. However, when continuous variables are scaled to
$(0,1)$, the scaling-related problems of MS lead to a rapid decline in
performance as $k_{c}$ increases.

Figure 15: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$ where
$Z$ has $dim_{c}=1,dim_{d}=2$ and $n_{c}=3$ with scaling to $(0,1)$ and rank
preprocessing. The distributions are depicted as described in Fig. 9.

Figure 16: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$ where
$Z$ has $dim_{c}=1,dim_{d}=2$ and $n_{c}=3$ with standardization
preprocessing. The distributions are depicted as described in Fig. 9.

Figure 17: The true null and permuted distributions of CMI estimated with the
MS, MS0-∞ and ZMADG estimators for the "Confounder"-model with $n=1000$ where
$Z$ has $dim_{c}=1,dim_{d}=2$ and $n_{c}=3$ with scaling to $(0,1)$ and rank
preprocessing. The distributions are depicted as described in Fig. 9.

Nevertheless, results for the same model with $dim_{d}=2$ and a larger sample
size of $n=2000$ (see Fig. 18) indicate that the performance of our approach
considerably increases given enough samples. Noticeably, in the case of a
larger sample size, a smaller $k_{c}=0.2$ gives the optimal results. Thus,
when applying our estimator, we recommend users to consider both the number of
samples and the dimensionality (i.e., the number of clusters). To exemplify,
consider the case where $n=1000$ and $dim_{d}=2$: There are approximately
$1000/(2\cdot 3)\approx 166$ samples per cluster for $n_{c}=3$ and
approximately $1000/(2\cdot 4)=125$ samples for $n_{c}=4$. In this case, first
of all, $k_{c}$ should be high enough such that $k$ has is large enough.
Otherwise, as previously discussed in Sec. Choice of Heuristic for $k$, the
estimation can suffer from high variance. Second, the number of samples per
cluster should be high enough to obtain reliable results: Although our
approach improves bias towards 0 compared to MS, it is still subject to the
curse of dimensionality.

Figure 18: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Confounder"-model with
$n=2000$, where $Z$ has $dim_{c}=1,dim_{d}=2$, with standardization, scaling
to $(0,1)$ and rank preprocessing for the continuous variables and coupling
factor in the dependent case $w=0.5$.

For the "Confounder"-case where $Z$ has only discrete dimension, i.e.,
$dim_{c}=0$ and sample size $n=1000$ (Fig. 19 and 20), we observe that ZMADG
underperforms, either having high FPR or low TPR. MS and MS0-∞ perform
similarly for the case when $dim_{d}=1$ and $dim_{d}=2$ and $n_{c}=3$
(considering their respective optimal $k_{c}$). Our method suffers from the
curse of dimensionality for $dim_{d}=2$ and $n_{c}=4$. As always, the scaling-
related problems of MS lead to underpeformance when variables are scaled to
$(0,1)$ and $k_{c}$ has larger values.

Figure 19: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the Confounder-model with
$n=1000$ where $Z$ has $dim_{c}=0,dim_{d}=1$, with standardization, scaling to
$(0,1)$ and rank preprocessing for the continuous variables and where $Z$ has
$dim_{c}=0,dim_{d}=2$ with standardization preprocessing for the continuous
variables. Both configurations have coupling factor in the dependent case
$w=0.5$.

Figure 20: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Confounder"-model with
$n=1000$ where $Z$ has $dim_{c}=0,dim_{d}=2$, with standardization, scaling to
$(0,1)$ and rank preprocessing for the continuous variables and coupling
factor in the dependent case $w=0.5$. Figure 21: False positive rate (FPR,
ideally under 0.05) and true positive rate (TPR, higher is better, with 1
best) for the "Independent $Z$"-model with $n=1000$ where $Z$ has $dim_{d}=1$.
Here we present results with rank preprocessing for the continuous variables
and coupling factor in the dependent case $w=0.5$.

### "Independent $Z$" Model

Here we present further results for the "Independent $Z$"-model with $n=1000$.

For the model with one discrete dimension, the MS and MS0-∞ CITs using rank
preprocessing for the continuous variables (Fig. 21) behave similarly to the
CITs with standardizatio. Considering their optimal $k_{c}$ values, the CIT
using our estimator performs similarly to the CIT using the MS estimator.

Figure 22: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Independent $Z$" model
where $Z$ has $dim_{d}=2$, with standardization, scaling to $(0,1)$ and rank
preprocessing for the continuous variables and coupling factor in the
dependent case $w=0.5$.

For the model with higher dimensionality (Fig. 22), we observe that the
behavior of the CITs is aligned with the behavior of the CITs for the
"Confounder"-model: our MS0-∞ CIT has slightly lower performance than MS, and
the performance gap increases for $n_{c}=4$ due to the curse of
dimensionality. Similar to the "Confounder"-model, the "Independent $Z$"-model
also has the characteristic of dependence across clusters. Consequently, the
performance of MS with standardization or rank preprocessing is not affected
if neighbors are taken from distinct clusters. However, due to the scaling-
related problems of MS, our estimator performs better for the scaling to
$(0,1)$ case with $dim_{d}=2$.

### "Cluster-dependent Confounder" Model

Here, we present additional results for the "Cluster-dependent Confounder"
model with $n=1000$. Consistent with the observations from previous models,
the CITs with rank preprocessing (Fig. 23) behave similarly to the CITs using
standardization. Nevertheless, there is a slightly more significant
performance gap between MS and our MS0-∞, with our method showcasing superior
performance, especially for $n_{c}=2$ and $n_{c}=3$.

Figure 23: False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Cluster-dependent"-model
where $Z$ has $dim_{d}=1$. Here we present results with rank preprocessing for
the continuous variables and coupling factor in the dependent case $w=0.75$.

### "Chain" Model

For the "Chain"-model (Fig. 24 and 25), ZMADG consistently suffers from high
FPR. In contrast, MS and MS0-∞ demonstrate good performance across varying
$k_{c}$ values and dimensionalities regarding TPR. The MS and MS0-∞ CITs
perform similarly for the standardization and rank preprocessing. In cases
where $n_{c}=4$, our approach has slightly lower TPR. Nonetheless, both CITs
control FPR effectively, except when $k_{c}=0.01$. The scaling-related
problems of MS lead to performance issues when variables are scaled to
$(0,1)$: an increase in $k_{c}$ leads to elevated FPR and decreased TPR. Our
approach demonstrates more robust performance across varying $k_{c}$ values.

Figure 24:  False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Chain"-model, with
standardization preprocessing for the continuous variables and coupling factor
in the dependent case $w=0.5$.

Figure 25:  False positive rate (FPR, ideally under 0.05) and true positive
rate (TPR, higher is better, with 1 best) for the "Chain"-model, with scaling
to $(0,1)$ and rank preprocessing for the continuous variables and coupling
factor in the dependent case $w=0.5$.

## Mean, Variance and Computational Time for the CMI Estimation Experiments

### Comparison of the Mean and Variances of the CMI

As outlined in Sec. 3 of the main paper, our estimator addresses the
challenges confronted by the MS and ZMADG estimators. Consequently, we expect
the following outcomes

  * •

for some of the models, there will likely be significant differences between
the bias of the MS and MS0-∞ estimators due to our estimator’s capacity to
reduce the bias towards 0,

  * •

and significant differences between the variance of the MS0-∞ estimator and
the variance of the ZMADG estimator for most of the models. These differences
arise from the definition of our estimator, which, in contrast to the ZMADG
estimator, is not an aggregation of multiple estimators.

Since bias and variance are not scaled metrics, we perform statistical tests
to investigate whether there are statistically significant differences in the
bias and variance among the various CMI estimators when measured on the models
presented in Sec.4 of the main paper and Sec. Further Results of the Numerical
Evaluation of the CMI Estimators of the SM. However, it is not immediately
apparent how exactly to conduct these tests. One possible approach would be to
compare the three estimators pairwise per $k_{c}$ value and sample size.
However, the MS and ZMADG estimators perform very differently across $k_{c}$
values; thus, it would not make for a fair comparison.

Hence, we conduct statistical tests to compare the estimators’ bias and
variances as follows: For each model, we identify the $k_{c}$ value with
optimal performance regarding bias and variance for the $100$ repetitions of
the experiment across all sample sizes. Subsequently, to test our hypothesis,
we compare the bias of the MS estimate versus the bias of the MS0-∞ estimate
and the variance of the MS0-∞ versus the variance of the ZMADG estimator.

To compare the bias, for each data model and each estimator
$estim\in\\{MS,MS_{0-\infty}\\}$, we compute the mean absolute error (in
short, MAE) $MAE_{estim,j}$ for each repetition of the experiment
$j\in\\{1,\ldots,100\\}$ as the absolute difference between the estimated CMI
using estimator $estim$ and the actual ground truth CMI value $I(X;Y|Z)$
corresponding to the given data model:

$MAE_{estim,j}=|\hat{I}^{estim}_{j}(X;Y|Z)-I_{j}(X;Y|Z)|.$ (37)

We then apply the Wilcoxon-Signed-Rank Test (Wilcoxon 1945) to compare bias
between estimators. As previously stated, we expect that for some of the data
models, the MAE of the MS estimator is higher than the MAE of our MS0-∞
estimator. We thus perform one-tailed tests. For each model and each sample
size, we formulate the hypothesis for the comparison of the two estimators as
follows:

$H_{0}$: The median of the difference between the $MAE_{MS}$ and
$MAE_{MS_{0-\infty}}$ is negative.

Vs.

$H_{1}$: The median of the difference between the $MAE_{MS}$ and
$MAE_{MS_{0-\infty}}$ is non-negative.

To compare the variance of the $100$ estimates $\hat{I}^{estim}_{j}(X;Y|Z)$
obtained with the MS0-∞ and ZMADG estimators, we test for equality of variance
using Levene’s Test (Levene 1960). All performed tests are two-tailed tests.
Thus, for each model and each sample size, we test the following hypothesis:

$H_{0}$: The variance of the MS0-∞ estimator and the variance of the ZMADG
estimator are equal.

Vs.

$H_{1}$: The variance of the MS0-∞ estimator and the variance of the ZMADG
estimator are not equal.

To account for repeated testing, we apply Bonferroni correction (Bonferroni
1935) by splitting the significance level by the number of hypotheses, in our
case $48$, and thus reject the null hypothesis for p-values under
$\frac{0.05}{48}=0.001$.

In the Tables 1 to 6 below, we report the obtained p-values for the
statistical tests for each model, measurement, and sample size. We observe
that the MAE of the MS estimator is significantly greater than the MAE of our
MS0-∞ estimator for the "Independent $Z$" with $d=3$ and the "Confounder with
uniform $X$ and $Y$" models . This aligns with our expectations, as we have
hypothesized in the Sec. 3 of the main paper: Our estimator should suffer from
less bias towards $0$. Furthermore, for almost all models and all sample
sizes, the hypothesis of equality of variances of our estimator and the ZMADG
estimator can be rejected, as expected and discussed in Sec. 3 of the main
paper.

p-values for the "Independent $Z$" Model with $d=1$
---
$n$ | Bias MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 0.9993 | 0.0000
600 | 1.0000 | 0.0000
1000 | 1.0000 | 0.0000
2000 | 1.0000 | 0.0000

Table 1: Results of the statistical tests for the "Independent $Z$" model with
$d=1$. We select $k_{c}$ for the individual estimators as follows:
$k_{c,MS}=0.01$, $k_{c,MS_{0-\infty}}=0.2$ and $k_{c,ZMADG}=0.1$.

p-values for the "Independent $Z$" Model with $d=3$
---
$n$ | Bias MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 0.0146 | 0.0000
600 | 0.0000 | 0.0000
1000 | 0.0000 | 0.0000
2000 | 0.0000 | 0.0000

Table 2: Results of the statistical tests for the "Independent $Z$" model with
$d=3$. We select $k_{c}$ for the individual estimators as follows:
$k_{c,MS}=0.01$, $k_{c,MS_{0-\infty}}=0.2$ and $k_{c,ZMADG}=0.1$.

p-values for the "Chain structure" Model with $d=1$
---
$n$ | Bias MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 1.0000 | 0.0000
600 | 1.0000 | 0.0000
1000 | 1.0000 | 0.0000
2000 | 1.0000 | 0.0000

Table 3: Results of the statistical tests for the "Chain structure" model with
$d=1$. We select $k_{c}$ for the individual estimators as follows:
$k_{c,MS}=0.01$, $k_{c,MS_{0-\infty}}=0.2$ and $k_{c,ZMADG}=0.2$.

p-values for the "Chain structure" Model with $d=3$
---
$n$ | Bias MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 0.9996 | 0.0000
600 | 1.0000 | 0.0000
1000 | 1.0000 | 0.0000
2000 | 1.0000 | 0.0000

Table 4: Results of the statistical tests for the "Chain structure" model with
$d=3$. We select $k_{c}$ for the individual estimators as follows:
$k_{c,MS}=0.01$, $k_{c,MS_{0-\infty}}=0.2$ and $k_{c,ZMADG}=0.2$.

p-values for the "Confounder with uniform $X$, $Y$" Model
---
$n$ | Mean MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 0.0000 | 0.0000
600 | 0.0000 | 0.0000
1000 | 0.0000 | 0.0000
2000 | 0.0000 | 0.0000

Table 5: Results of the statistical tests for the "Confounder with uniform $X$
and $Y$" model. We select $k_{c}$ for the individual estimators as follows:
$k_{c,MS}=0.1$, $k_{c,MS_{0-\infty}}=0.1$ and $k_{c,ZMADG}=0.3$.

p-values for the "Confounder with Gaussian $X,Y$" Model
---
$n$ | Bias MS vs. MS0-∞ | Var MS0-∞ vs. ZMADG
300 | 0.5000 | 0.0000
600 | 0.5000 | 0.0000
1000 | 0.5000 | 0.0000
2000 | 0.5000 | 0.0000

Table 6: Results of the statistical tests for the "Confounder with Gaussian
$X$ and $Y$" model. We select $k_{c}$ for the individual estimators as
follows: $k_{c,MS}=0.01$, $k_{c,MS_{0-\infty}}=0.3$ and $k_{c,ZMADG}=0.1$.

### Report on Computational Runtimes

Here, we report the computational runtimes for each model in the CMI
estimation experiments,. The experiments were performed on an Intel(R)
Core(TM) i7-6600U CPU. We recorded the runtimes for each run of the $100$
repetitions of the CMI experiments using the different $k_{c}$ values and
sample sizes $n$. In the Tables 7 to 12, we present the average runtime for
each sample size $n$ in seconds, averaged over the $100$ runs and all $k_{c}$
values. We also include the average runtime of the different heuristics for
our approach: "local", "global", and "cluster". We note that in the tables
below, MS0-∞ stands for the MS0-∞ in combination with the "MS0-∞ local"
heuristic for our estimator in combination with the "local" heuristic, and
"MS0-∞ cluster" stands for our estimator in combination with the "cluster-
size" heuristic.

$n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---
300 | 0.013 | 0.016 | 0.013 | 0.029 | 0.036
600 | 0.033 | 0.023 | 0.034 | 0.03 | 0.046
1000 | 0.069 | 0.046 | 0.076 | 0.05 | 0.046
2000 | 0.249 | 0.142 | 0.271 | 0.143 | 0.05

Table 7: Computational runtimes (in seconds) for the "Independent $Z$" model
(Sec. 4 of the main paper) with $d=1$.

$n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---
300 | 0.021 | 0.018 | 0.02 | 0.025 | 0.059
600 | 0.052 | 0.043 | 0.051 | 0.037 | 0.087
1000 | 0.107 | 0.081 | 0.111 | 0.061 | 0.118
2000 | 0.327 | 0.215 | 0.333 | 0.165 | 0.164

Table 8: Computational runtimes (in seconds) for the "Independent $Z$" model
(Sec. 4 of the main paper) with $d=3$.

$n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---
300 | 0.012 | 0.015 | 0.012 | 0.029 | 0.014
600 | 0.036 | 0.023 | 0.035 | 0.033 | 0.021
1000 | 0.082 | 0.05 | 0.08 | 0.053 | 0.027
2000 | 0.332 | 0.171 | 0.306 | 0.136 | 0.039

Table 9: Computational runtimes (in seconds) for the "Chain structure" model
(Sec. Further Results of the Numerical Evaluation of the CMI Estimators) with
$d=1$.

| $n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---|---
0 | 300 | 0.019 | 0.011 | 0.019 | 0.024 | 0.027
1 | 600 | 0.057 | 0.024 | 0.055 | 0.035 | 0.041
2 | 1000 | 0.116 | 0.041 | 0.113 | 0.052 | 0.054
3 | 2000 | 0.433 | 0.097 | 0.399 | 0.131 | 0.079

Table 10: Computationals runtime (in seconds) for the "Chain structure" model
(Sec. Further Results of the Numerical Evaluation of the CMI Estimators) with
$d=3$.

$n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---
300 | 0.01 | 0.012 | 0.01 | 0.025 | 0.023
600 | 0.03 | 0.023 | 0.029 | 0.027 | 0.023
1000 | 0.068 | 0.051 | 0.068 | 0.042 | 0.026
2000 | 0.253 | 0.183 | 0.252 | 0.115 | 0.037

Table 11: Computational runtimes (in seconds) for the "Confounder with uniform
$X$ and $Y$" model (Sec. Further Results of the Numerical Evaluation of the
CMI Estimators).

$n$ | MS0-∞ | MS0-∞ global | MS0-∞ cluster | MS | ZMADG
---|---|---|---|---|---
300 | 0.009 | 0.013 | 0.009 | 0.027 | 0.022
600 | 0.027 | 0.021 | 0.026 | 0.029 | 0.023
1000 | 0.06 | 0.046 | 0.059 | 0.043 | 0.028
2000 | 0.221 | 0.156 | 0.216 | 0.098 | 0.04

Table 12: Computational runtimes (in seconds) for the "Confounder with
Gaussian $X$ and $Y$" model (Sec. Further Results of the Numerical Evaluation
of the CMI Estimators).

## Remarks on Reproducibility

We intentionally reduced the number of plots and tables in this supplementary
material for length reasons. However, all code to obtain the CIT plots and the
measurements for the statistical significance tests, as well as the
computational time reports can be found in the .zip file accompanying this
supplementary, including the details necessary for replicating our
experiments, such as random seeds.