# Measuring association with recursive rank binning

Chris Salahub  Department of Statistics and Actuarial Science, University of
Waterloo Wayne Oldford  Department of Statistics and Actuarial Science,
University of Waterloo

###### Abstract

Pairwise measures of dependence are a common tool to map data in the early
stages of analysis with several modern examples such as those in Reshef et al.
(2011) and Heller et al. (2016) based on maximized partitions of the pairwise
sample space. Following a short survey of modern measures of dependence, we
introduce a new measure which recursively splits the ranks of a pair of
variables to partition the sample space and computes the $\chi^{2}$ statistic
on the resulting bins. Splitting logic is detailed for splits maximizing a
score function and randomly selected splits. Simulations indicate that random
splitting produces a statistic conservatively approximated by the $\chi^{2}$
distribution without a loss of power to detect numerous different data
patterns compared to maximized binning. Though it seems to add no power to
detect dependence, maximized recursive binning is shown to produce a natural
visualization of the data and the measure. Applying maximized recursive rank
binning to S&P 500 constituent data suggests the automatic detection of tail
dependence.

##  1 Introduction

The modern world is awash with data. Ubiquitous data collection renders data
sets comprised of ever more variables and observations which must be sifted
through for insight. This is true, for example, in finance and genetics, where
real time measurement and the advent of complete sequencing have resulted in
truly massive raw data. In finance, the identification of interesting pairwise
variable relationships can motivate strategy or hedge risk, while the
comparison of many genetic markers to a physical trait confers greater
understanding of inheritance and the nature of some disease. How
“interestingness” is measured in these and many other applications takes many
different forms, such as the scagnostics of Wilkinson et al. (2005).

Perhaps the most widely used measure of interestingness is statistical
dependence. A broad measure of the interestingness of a relationship between
variables, it is defined as the absence of statistical independence and so
includes a wide array of patterns. Consequently, many different measures
devised for specific contexts and patterns of interest have been developed
(e.g. Goodman and Kruskal, 1979; Liebetrau, 1983; Choi et al., 2010; Tjøstheim
et al., 2022) since the early efforts of Galton and Pearson to characterize it
with correlation and the $\chi^{2}$ test (Stigler, 1989; Hald, 1998).

The eminence of statistical dependence is driven by the observation that it is
almost universally interesting in an exploratory context. When sifting through
many variable pairs to select those worthy of further investigation, a
statistical dependency in a pair provides information about how the two
variables should be controlled or modelled in any investigation.

The rise of computers has changed the nature of these proposals drastically,
resulting in more flexible and computationally complex measures than ever
before. Many modern measures of dependence are products of the computer age,
such as algorithmic ones based on nearest neighbours (Dümcke et al., 2014),
prediction (Breiman, 2001), or the maximization of a particular data
transformation (Breiman and Friedman, 1985; Reshef et al., 2011; Jiang et al.,
2015; Liu et al., 2018).

A common feature of many measures of statistical dependence is the evaluation
of functionals that compare the observed joint density to the product of
marginal densities, as these address the problem of dependence directly using
its definition.111The organization here is therefore somewhat different than
is typical. Many works such as Tjøstheim et al. (2022), an up-to-date survey
of meaures of dependence, separate measures based on copulas, kernel
functions, and partitioning. These are all cast here as particular instances
of functionals comparing the joint and product marginal densities to emphasize
similarities in the computation of these conceptually different measures in
practice. In particular, much recent work has addressed _bin-based_ methods
which partition the pairwise sample space, as these automatically perform non-
parametric estimation of the joint and product marginal densities (Reshef et
al., 2011; Heller et al., 2016; Cao et al., 2021).

This work introduces and investigates a new measure of dependence based on a
recursive partitioning, as in Rahman (2018). Section 2 outlines some necessary
notation and defines terms used in Section 3 to introduce a brief summary of
measures of interestingness and dependence. Section 4 follows this general
outline by focusing on bin-based measures of dependence and outlining previous
proposals in a common notation.

The details of a new recursive partitioning method to measure association are
described in Section 5, including a proofs that the coordinates of maximizing
splits occur at observations and a discussion of stop criteria and general
splitting logic. All of this is combined in the pseudo-code algorithms of
Section 6 before the details of a package implementing the method in R are
discussed in Section 7. Simulations carried out in Section 8 characterize the
null distribution of the algorithm for different sample sizes and indicate
that recursive binning powerfully detects many different patterns of
dependence, even when the split choice is made randomly. Finally, the
algorithm is applied to S&P 500 constituent data in Section 9 with results
that suggest maximized splits detect patterns of specific interest such as
tail dependence automatically.

##  2 Preliminaries

Let $X$ and $Y$ be a pair of random variables with distribution functions
$F_{X}(x)$ and $F_{Y}(y)$ over the domains $\mathcal{X}$ and $\mathcal{Y}$
both respectively. Denote the joint distribution of the pair $X,Y$ as
$F_{X,Y}(x,y)$, and the conditional distributions $F_{X|Y}(x|y)$ and
$F_{Y|X}(y|x)$. To remain fully general, note that both of $X$ and $Y$ may be
continuous or discrete. Let $f_{X}(x)$, $f_{Y}(y)$, $f_{X,Y}(x,y)$,
$f_{X|Y}(x|y)$, and $f_{Y|X}(y|x)$ be the corresponding probability densities
or probability mass functions in the continuous or discrete cases. We say that
$X$ and $Y$ are independent and write $X\bot\bot Y$ if and only if

$f_{X,Y}(x,y)=f_{X}(x)f_{Y}(y)$ (1)

or equivalently $f_{X|Y}(x|y)=f_{X}(x)$ and $f_{Y|X}(y|x)=f_{Y}(y)$.

Define a $K$-dimensional copula as a distribution function
$C(u_{1},u_{2},\dots,u_{K})$ over $[0,1]^{k}$ with uniform marginal
distributions for all $u_{k}$ as in Embrechts et al. (2001). By Sklar’s
Theorem, for any $K$ continuous random variables there exists a unique copula
$C$ such that

$F_{X_{1},X_{2},\dots,X_{K}}(x_{1},x_{2},\dots,x_{K})=C(F_{X_{1}}(x_{1}),F_{X_{2}}(x_{2}),\dots,F_{X_{K}}(x_{K})),$
(2)

that is the joint distribution can be summarized by a copula on the marginals
transformed to be uniformly distributed. In other words: the relationship
between the variables is uniquely determined by $C$ and the marginals alone.
Define the independence copula

$C_{I}(u_{1},u_{2},\dots,u_{K})=\prod_{k=1}^{K}u_{k}.$ (3)

In particular, consider the copula of $X$ and $Y$,
$C(F_{X}(x),F_{Y}(y))=C(u,v)$, and the bivariate independence copula,
$C_{I}(u,v)=uv$.

In practice these theoretical quantities are unknown. Instead, all that is
available is a sample of paired observations
$(x_{1},y_{1}),(x_{2},y_{2}),\dots,(x_{n},y_{n})$ of $(X,Y)$. Let
$\mathbf{x}={(x_{1},x_{2},\dots,x_{n})}^{\mkern-1.5mu\mathsf{T}}$ and
$\mathbf{y}={(y_{1},y_{2},\dots,y_{n})}^{\mkern-1.5mu\mathsf{T}}$ be the
observed values of $X$ and $Y$ respectively. For $X$ and $Y$ which can be
ordered, let a subscript $(x)$ indicate a non-decreasing sorting of elements
with respect to $\mathbf{x}$ so that $\mathbf{x}_{(x)}$ is the vector
$\mathbf{x}$ in increasing order and $\mathbf{y}_{(x)}$ is the vector
$\mathbf{y}$ sorted in increasing order of $\mathbf{x}$. Elementwise, follow
the convention

$\mathbf{x}_{(x)}={(x_{(1)},x_{(2)},\dots,x_{(n)})}^{\mkern-1.5mu\mathsf{T}}$

to denote the elements of $\mathbf{x}_{(x)}$ and analogously for
$\mathbf{y}_{(y)}$. Define the rank function on a marginal sample $\mathbf{x}$
as

$r(x;\mathbf{x})=\sum_{i=1}^{n}I_{(-\infty,x]}\left(x_{i}\right),$ (4)

where

$I_{A}(x)=\left\\{\begin{array}[]{rcl}1&{}{}{}{}\hfil&\mbox{if}~{}x\in A\\\
0&&\mbox{otherwise}\end{array}\right.$

is the usual indicator function so that $r(x_{i},\mathbf{x})$ gives the index
of $x_{i}$ in $\mathbf{x}_{(x)}$. Note that this definition assumes $x_{i}\neq
x_{j}$ for all $i\neq j$, the case of no _ties_.

To consider different conventions to address ties, suppose
$x_{1}=x_{2}=\dots=x_{m}=x_{(1)}$ for $m<n$. In this case,
$r(x_{1};\mathbf{x})=r(x_{2};\mathbf{x})=\dots=r(x_{m};\mathbf{x})=m$, that is
all are given the _maximum_ index $m$. It is not clear this must be the case,
as choosing the _minimum_ index $1$ seems equally valid for these first $m$
observations. Another option is _random tie-breaking_ , which randomly assigns
the indices $1,2,\dots,m$ to $x_{1},x_{2},\dots,x_{m}$. Random tie breaking is
of particular interest, as it induces a uniform distribution on the ranks for
tied regions. For applications where complete ranks are necessary, this is
more desirable than the gaps introduced by the maximum or minimum indexing.

We can then define the empirical distribution function for these ordered cases
as

$\widehat{F}_{\mathbf{x}}(x)=\frac{1}{n}r(x;\mathbf{x})=\frac{1}{n}\sum_{i=1}^{n}I_{(-\infty,x]}\left(x_{i}\right),$
(5)

with $\widehat{F}_{\mathbf{y}}(y)$ defined similarly. The empirical copula of
$X$ and $Y$ is defined as

$\widehat{C}(u,v)=\frac{1}{n}\sum_{i=1}^{n}I_{(-\infty,u]}\left(\frac{r(x_{i};\mathbf{x})}{n}\right)I_{(-\infty,v]}\left(\frac{r(y_{i};\mathbf{y})}{n}\right).$
(6)

Finally, define the vector-valued rank function as

$\mathbf{r}(\mathbf{x})={(r(x_{1};\mathbf{x}),r(x_{2};\mathbf{x}),\dots,r(x_{n};\mathbf{x}))}^{\mkern-1.5mu\mathsf{T}}.$
(7)

For discrete $X$ without ordered values, take an arbitrary numbering of the
possible values in $\mathcal{X}$ and let $I=|\mathcal{X}|$. Otherwise, number
such that the assigned values reflect the ordering. If $Y$ is discrete, number
$\mathcal{Y}$ similarly with $J=|\mathcal{Y}|$. Then, for $\mathbf{x}$ and
$\mathbf{y}$, define the observed marginal and joint counts

$o_{i+}=\sum_{j=1}^{J}\sum_{k=1}^{n}I_{\\{(i,j)\\}}\left((x_{k},y_{k})\right)=\sum_{j=1}^{J}o_{ij}$
(8)

and $o_{+j}$ analogously. These give corresponding probability estimates
$\widehat{p}_{ij}$, $\widehat{p}_{i+},$ and $\widehat{p}_{+j}$ when divided by
$n$.

Distinguish three different _types_ of $X$ and $Y$ based on $\mathcal{X}$ and
$\mathcal{Y}$. The three types are:

continuous

some interval or collection of intervals of real numbers;

ordinal categorical

discrete possibilities with an ordering (e.g. income brackets);

nominal categorical

discrete possiblities without an ordering (e.g. country of birth).

The possible combinations of these lead to six unique pairwise combinations.
Following the convention of Lee and Huh (2003), call comparisons of data of
the same type _simple_ comparisons. _Complex_ comparisons refer to the unique
pairings between types: continuous versus ordinal, continuous versus nominal,
and ordinal versus nominal.

Introduce the function

$G:\mathcal{X}\times\mathcal{Y}\mapsto\mathcal{R}$

quantifying some pattern of interest between $X$ and $Y$ with a value in
$\mathcal{R}\subset\mathbb{R}$. For the purpose of this work, any such $G$ is
called a measure of interestingness or simply a _measure_. The term
“interestingness” is used instead of “association” following the tradition of
researchers such as John W. Tukey (see Friedman and Stuetzle (2002)) and to
avoid the implications of “association”, typically used to describe
_statistical dependence_ as in Equation (1).

Consider, for example, Tukey’s scagnostic measures of _structure_ as
reintroduced in Wilkinson et al. (2005). Two different constructed data sets
exemplifying the scagnostics for “stringiness” and “clumpiness” respectively
are shown in Figure 1.

|
---|---
(a) Striated/stringy sample | (b) Clumpy sample
Figure 1: Data exhibiting certain scagnostic structures despite the
independent generation of $x$ and $y$.

Both of these simulated data sets have margins generated independently of each
other, and still have relatively large values for their respective measures
compared to the null distributions in Wilkinson and Wills (2008). While
scagnostics have been expressly developed to capture patterns other than
statistical independence, calling them measures of association confuses the
term.222See, for example Liu et al. (2018) utilizing data analogous to Figure
1(b) as an example of ‘association’ in a paper motivated by statistical
independence. Another example is $\lambda_{b}$ from Goodman and Kruskal
(1979), which is motivated by prediction and not strictly tied to dependence.

Numerous frameworks have been proposed to evaluate candidates for $G$. A
statistical framework is given in Rényi (1959), and this framework has been
subsequently updated by Schweizer and Wolff (1981) and Móri and Székely
(2019). Reimherr and Nicolae (2013) present a slightly different framework
which emphasizes interpretability and delineates three different motivations
for measuring association. Here, the focus is on the _ordering_ $G$ induces on
$\mathcal{X}\times\mathcal{Y}$, as this can accomplish all of the goals from
Reimherr and Nicolae (2013) and is critical to the search for interesting
variables.

The importance of ordering restricts most measures of interestingness to a
range on a finite interval
$\mathcal{R}=[g_{\text{min}},g_{\text{max}}]\subset\mathbb{R}$. The upper
bound is obtained for $X$ and $Y$ exemplary of the pattern of interest, while
$g_{\text{min}}$ is less consistent. For signed measures, such as Pearson’s
correlation coefficient, $g_{\text{min}}$ may still indicate perfect
correspondence to a particular pattern. In unsigned measures $g_{\text{min}}$
suggests no indication of the pattern333Any signed measure $G$ can, of course,
be made unsigned by taking $\lvert{G}\rvert$.. The $X$ and $Y$ which lead to
these extremes are typically not unique, but rather represent a family of
patterns which the measure does not distinguish.

Commonly, $G(X,Y)$ is scaled by
$\max\left\\{\lvert{g_{\text{max}}}\rvert,\lvert{g_{\text{min}}}\rvert\right\\}$
such that $\mathcal{R}$ is restricted to $[-1,1]$ or $[0,1]$. By scaling $G$
by the magnitude of its most extreme value, the ordering it imposes is made
explicit. For any pair $X$ and $Y$, this scaled measure communicates directly
how the pair compares to the perfect example. This scaling can be misleading,
however, when the analyst is unfamiliar with its distribution along its range.

##  3 Measuring dependence

As noted by Cramér (1924); Fairfield Smith (1957); Goodman and Kruskal (1979)
and likely others, measures are devised to capture specific patterns. This is
often desirable for interpretability, as Reimherr and Nicolae (2013) note, and
can lead to descriptive adjectives which evoke a measure, as in the
scagnostics of Wilkinson et al. (2005). Some common patterns and corresponding
measures are:

linearity

which is typically measured by Pearson’s product moment correlation;

monotonicity

captured by Spearman’s $\rho$, the rank version of correlation;

concordance

measured by Kendall’s $\tau$ and Goodman and Kruskal’s $\gamma$ (from Goodman
and Kruskal (1979));

predictibility

measured by Goodman and Kruskal’s $\lambda$ (also from Goodman and Kruskal
(1979));

agreement

addressed by the $\kappa$ measure of inter-rater reliability from Cohen
(1960);

and others listed in Liu et al. (2018); Liebetrau (1983); Agresti (1981); and
Goodman and Kruskal (1979).

The variety of patterns which one might want to measure and the different
types of $X$ and $Y$ have led to a great proliferation of bivariate measures
of interestingness444The 2$\times$2 contingency table, for example, has an
almost overwhelming roster of measures analyzed in Choi et al. (2010)., each
with its own interpretation. Unfortunately, the division of measures by both
type and pattern makes the search for interesting patterns in complex data far
more challenging, with research such as Khamis (2008) and Lee and Huh (2003)
devoted to the task of guiding practitioners to choose the most comparable
measures between types. In practice, correlation is often applied to all
variable pairs regardless of type, despite its well-known shortcomings as
outlined in Tjøstheim et al. (2022); Reshef et al. (2011); and Anscombe
(1973).

Regardless of motivation or type, however, $X$ and $Y$ satisfying Equation (1)
are universally considered uninteresting. As a consequence, complex measures
of interestingness typically attempt to measure departures from statistical
independence. These measures of statistical dependence are functionals of the
form

$D(f_{X,Y}(x,y),f_{X}(x)f_{Y}(y))$ (9)

which compare the joint and marginal product distributions of $X$ and $Y$.
Often, they will attempt to satisfy the desiderata outlined by Rényi (1959) or
Schweizer and Wolff (1981), the latter of which relaxed the axioms of the
former after noting they are unnecessarily restrictive and hard to apply in
practice. Some examples for continuous $X$ and $Y$ include

$\varDelta(X,Y)=\int\int_{\\{(x,y):f_{X,Y}(x,y)\geq
f_{X}(x)f_{Y}(y)\\}}\big{[}f_{X,Y}(x,y)-f_{X}(x)f_{Y}(y)\big{]}dxdy$ (10)

from Silvey (1964) and mutual information

$\mathcal{I}(X,Y)=D_{KL}(f_{X,Y}||f_{X}f_{Y})=\int_{\mathcal{Y}}\int_{\mathcal{X}}f_{X,Y}(x,y)\log\left(\frac{f_{X,Y}(x,y)}{f_{X}(x)f_{Y}(y)}\right)dxdy$
(11)

from Shannon (1948), where $D_{KL}(F||G)$ denotes the Kullback-Leibler
divergence of $G$ from $F$. For discrete $X$ and $Y$ a classic example is the
$\chi^{2}$ statistic for independence from Pearson (1900)555Note that this
formulation of $\mathcal{D}(X,Y)$ is somewhat atypical, as this is a quantity
usually defined on the sample $\mathbf{x},\mathbf{y}$.

$\mathcal{D}(X,Y)\propto\sum_{x\in\mathcal{X}}\sum_{y\in\mathcal{Y}}\frac{[f_{X,Y}(x,y)-f_{X}(x)f_{Y}(y)]^{2}}{f_{X}(x)f_{Y}(y)}.$
(12)

Another possible functional based on the cumulative distribution functions is
given by Hoeffding (1948)

$D_{H}(X,Y)=\int\left(F(x,y)-F(x,\infty)F(\infty,y)\right)^{2}dF(x,y).$ (13)

Central to all of these measures is the factorization definition of
independence from Equation (1). In every case, some comparison is made between
$f_{X,Y}(x,y)$ and $f_{X}(x)f_{Y}(y)$ which is zero when the two are equal
almost everywhere.

For any function of the form in Equation (9), an equivalent measure exists in
the copula space. By Sklar’s Theorem, discussed in Sklar (1996) and summarized
in Embrechts et al. (2001), the dependence between $X$ and $Y$ can be captured
by the marginal distribution functions $F_{X}$ and $F_{Y}$ and their copula
$C$. Just as in Equation (1), $X$ and $Y$ are independent only if their copula
is the independence copula $C_{I}(u,v)=uv$. Schweizer and Wolff (1981)
therefore show that a transform can be applied to any measure formulated in
$F_{X}(x)$, $F_{Y}(y)$, and $F_{X,Y}(x,y)$ to convert it to a functional of
the form

$D^{*}(C(F_{X}(u),F_{Y}(v)),uv).$

Such transforms are attractive because they remove the marginal distributions
of the variables being compared. This allows for non-parametric estimation of
quantities without presuming any particular distribution. Indeed, there are
many proposed measures which utilize empirical copulas, such as those
presented in Ding et al. (2017); Siburg and Stoimenov (2010); and Genest and
Rémillard (2004).

Measures like the distance covariance of Székely and Rizzo (2009) instead
apply a functional of the form in Equation (9) to the characteristic functions
of $X$ and $Y$. The characteristic function of $X$ is defined as

$\phi_{X}(t)=E[e^{itX}],$ (14)

where $E$ is the expectation operator with respect to $X$. The joint
characteristic function, $\phi_{X,Y}(t,s)$, and characteristic function of
$Y$, $\phi_{Y}(s)$, are defined similarly. This is based on an important
result of Equation (1), that $\phi_{X,Y}(s,t)=\phi_{X}(t)\phi_{Y}(s)$ if and
only if $X\bot\bot Y$.

Equation (9) is really just an evaluation of the goodness of fit of
$f_{X}(x)f_{Y}(y)$ to $f_{X,Y}(x,y)$, and this creates an obvious analogy to
empirical distribution function goodness of fit tests as described in Zheng et
al. (2021) and Stephens (1974). Many common measures of functional distance
appear in the literature of independence tests, such as the Kolmogorov-Smirnov
test in Heller et al. (2016) and results for any $L_{P}$ distance in Schweizer
and Wolff (1981). These can be made even more general by replacing Euclidean
distances by kernel distances, as in Liu et al. (2018) or the Hilbert-Schmit
independence criterion (HSIC) of Gretton et al. (2007). While Liu et al.
(2018) maximize over a pre-specified set of kernel distances, Lopez-Paz et al.
(2013) instead introduce a measure which randomly transforms $X$ and $Y$ and
takes the maximum of the applied random transforms.

More unique applications of goodness of fit principles are found in Dümcke et
al. (2014) and Heller et al. (2013). Dümcke et al. (2014) utilize the exact
distribution of nearest neighbour distances under independence to develop two
novel tests. Rather than comparing the joint distribution to a product of
marginals, their tests are based on the deviation between the exact
distribution of nearest neighbours and that observed in a sample. Heller et
al. (2013) make use of local distances about each point in turn to construct a
series of contingency tables and then aggregate the $p$-values gained.

The generation of numerous tables evokes a relevant class of measures based on
the application of measures as in Equation (9) to partitions of the outcome
space $\mathcal{X}\times\mathcal{Y}$, after which the values are aggregated.
Such partitioning essentially allows for local estimation of $f_{X,Y}$ and
$f_{X}f_{Y}$ without a parametric family. Rather, the estimate relies only on
the choice of partition. Another benefit of such partitioning is it permits
the application of the same test to any data type, as partitioning continuous
data produces ordinal categorical data. As a consequence of this potential and
their popularity in recent literature, this work primarily focuses on these
methods.

##  4 Bin-based measures

Though it has become more popular as it has become computationally feasible,
partitioning, or _binning_ , has always had a role in computing measures of
association. See, for example, the early investigations into the $\chi^{2}$
test outlined in Plackett (1983). Rather than using $X$ and $Y$ directly, a
binning function can be applied to their marginal values in order to
discretize them.666In the case where one of $X$ and $Y$ is already
categorical, this reduces to the $K$-sample problem. Some proposals, such as
Heller et al. (2016), switch freely between this $K$-sample problem and the
problem of measuring association. Therefore, discussing a measure which
_induces_ bins only makes sense when at least one of $X$ and $Y$ is
continuous, though the results can be applied to categorical cases.

A univariate binning on $J$ bins is a function
$b:\mathcal{B}\mapsto\\{1,\dots,J\\}$ which partitions its continuous domain
$\mathcal{B}\subseteq\mathbb{R}$ into $J$ distinct parts, or _bins_. Any such
$b$ has a vector-valued version
$\mathbf{b}:\mathcal{B}^{p}\mapsto\\{1,\dots,J\\}^{p}$ such that
$\mathbf{b}(\mathbf{x})={(b(x_{1}),\dots,b(x_{n}))}^{\mkern-1.5mu\mathsf{T}}$
for $\mathbf{x}\in\mathcal{B}^{p}$. Consider applying
$b_{X}:\mathcal{X}\mapsto\\{1,\dots,I\\}$ to $X$ and
$b_{Y}:\mathcal{Y}\mapsto\\{1,\dots,J\\}$ to $Y$. That is apply a binning on
$I$ bins to $X$ and a binning on $J$ bins to $Y$. This is equivalent to an
$I\times J$ grid on $\mathcal{X}\times\mathcal{Y}$, and so the bins can be
indexed by $(i,j)$ to correspond with the values $(b_{X}(X),b_{Y}(Y))$. Define

$\epsilon_{ij}=n\int_{\\{x:b_{X}(x)=i\\}}dF_{X}(x)\int_{\\{y:b_{Y}(y)=j\\}}dF_{Y}(y),$
(15)

the expected count of observations of $n$ which fall into bin $(i,j)$ under
independence. Note that in the case of uniform $X$ and $Y$, this simplifies to
$n$ times the area of the $(i,j)$ bin.

In practice the vectors $\mathbf{x}$ and $\mathbf{y}$ are all that is
observed, so define the analogous sample quantity

$e_{ij}=\frac{1}{n}\sum_{k=1}^{n}I_{\\{i\\}}\left(b_{X}(x_{k})\right)\sum_{l=1}^{n}I_{\\{j\\}}\left(b_{Y}(y_{l})\right).$
(16)

The observed $(i,j)$ bin count is given by

$o_{ij}=\sum_{k=1}^{n}I_{\\{(i,j)\\}}\left(\big{(}b_{X}(x_{k}),b_{Y}(y_{k})\big{)}\right),$
(17)

so $e_{ij}=n\frac{o_{i+}}{n}\frac{o_{+j}}{n}$ using the notation of Equation
(8). Under this binning, $X$ and $Y$ are converted to the contingency table in
Table 1.

| $b_{Y}(y)=1$ | $b_{Y}(y)=2$ | $\dots$ | $b_{Y}(y)=J$ |
---|---|---|---|---|---
$b_{X}(x)=1$ | $o_{11}$ | $o_{12}$ | $\dots$ | $o_{1J}$ | $o_{1+}$
$b_{X}(x)=2$ | $o_{21}$ | $o_{22}$ | $\dots$ | $o_{2J}$ | $o_{2+}$
$\vdots$ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$ | $\vdots$
$b_{X}(x)=I$ | $o_{I1}$ | $o_{I2}$ | $\dots$ | $o_{IJ}$ | $o_{I+}$
| $o_{+1}$ | $o_{+2}$ | $\dots$ | $o_{+J}$ | $n$

,

Table 1: The contingency table imposed on $\mathbf{x}$ and $\mathbf{y}$ by the
binnings $b_{X}(x)$ applied to $\mathbf{x}$ and $b_{Y}(y)$ applied to
$\mathbf{y}$.

Once so binned, the $o_{ij}$ and $e_{ij}$ can be used to in place of the
unknown densities in Equations (10), (11), and (12) nonparametrically. In this
way, a grid is simply a particular kind of two-dimensional histogram.777There
are many other possible tessellations which produce a two-dimensional
histogram density estimate, see Carr et al. (1987) and Scott (1988) For
Silvey’s $\varDelta$ from Equation (10), the analogue on the binned data in
Table 1 is

$\varDelta(\mathbf{b}_{\mathbf{x}},\mathbf{b}_{\mathbf{y}})=\sum\sum_{\\{(i,j):o_{ij}\geq
e_{ij}\\}}\frac{o_{ij}-e_{ij}}{n}=\sum\sum_{\\{(i,j):o_{ij}\geq\frac{1}{n}o_{i+}o_{+j}\\}}\frac{no_{ij}-o_{i+}o_{+j}}{n^{2}};$
(18)

Shannon’s mutual information from Equation (11) becomes the multinomial log-
likelihood ratio

$\mathcal{I}(\mathbf{b}_{\mathbf{x}},\mathbf{b}_{\mathbf{y}})=\sum_{i=1}^{I}\sum_{j=1}^{J}\frac{o_{ij}}{n}\log\left(\frac{o_{ij}}{e_{ij}}\right)=\sum_{i=1}^{I}\sum_{j=1}^{J}\frac{o_{ij}}{n}\log\left(\frac{no_{ij}}{o_{i+}o_{+j}}\right);$
(19)

and Pearson’s $\chi^{2}$ measure from Equation (12) becomes

$\mathcal{D}(\mathbf{b}_{\mathbf{x}},\mathbf{b}_{\mathbf{y}})=\sum_{i=1}^{I}\sum_{j=1}^{J}\frac{(o_{ij}-e_{ij})^{2}}{e_{ij}}=\frac{1}{n}\sum_{i=1}^{I}\sum_{j=1}^{J}\frac{(no_{ij}-o_{i+}o_{+j})^{2}}{o_{i+}o_{+j}},$
(20)

where $\mathbf{b}_{\mathbf{x}}=\mathbf{b}_{X}(\mathbf{x})$ and
$\mathbf{b}_{\mathbf{y}}=\mathbf{b}_{Y}(\mathbf{y})$ are used for brevity.
While $\mathbf{b}_{X}$ and $\mathbf{b}_{Y}$ can be arbitrarily defined,
Equations (18), (19), and (20) depend only on the counts in each cell of Table
1. Therefore, only the values of $\mathbf{x}$ and $\mathbf{y}$ need to be
considered as possible partition boundaries, or _bin edges_. Considering edges
between every point, and not allowing for identical $\mathbf{x}$ and
$\mathbf{y}$ values, this means are $n-1$ possible bin edges to consider in
each of $\mathbf{x}$ and $\mathbf{y}$. Noting that the bin edges can either be
present or absent, there are therefore $2^{n-1}2^{n-1}=4^{n-1}$ possible
contingency tables for a given $\mathbf{x}$ and $\mathbf{y}$.

There is often no a priori exploratory reason to choose a particular binning,
and so methods from Reshef et al. (2011); Jiang et al. (2015); Heller et al.
(2016); and Reshef et al. (2018) have the exploration of this large number of
grids at their core. Reshef et al. (2011) and Reshef et al. (2018), in the
introduction of the maximal information criterion (MIC), propose computing
Equation (19) for all grids such that $IJ\leq n^{0.6}$, scaling these values,
and storing them in a matrix $M$. The maximal value of this matrix is then
taken as the MIC.888It should be noted that even this space is too large to
fully explore, and so only a small number of actual grids are computed by
their algorithm in practice. Their concept of _equitability_ has also
generated considerable controversy. See the discussions in Gorfine et al.
(2012); Kinney and Atwal (2014a); Reshef et al. (2014); Kinney and Atwal
(2014b); and Simon and Tibshirani (2014). In order to contextualize MIC
values, $p$-values based on simulated null data sets are computed.

Jiang et al. (2015) propose a penalized version of Equation (19) to find a
solution. The particular penalty is outlined in Equation (24). Conceptually
this penalized optimization assumes a Poisson distribution on the number of
bins for one of the margins conditioned on the other and then maximizes the
likelihood using a dynamic algorithm. As with the MIC, the null distribution
of this method is determined empirically.

Hoeffding (1948) restricts the exploration only to $I=J=2$, and proposes a sum
of Equation (13) for all unique possible two-by-two grids, which Thas and
Ottoy (2004) show is asymptotically equivalent to Equation (20) evaluated over
all two-by-two grids with a suitable scaling. Heller et al. (2016) extend this
by simply restricting the grid to $m$ divisions on both margins. They
investigate both summation and maximization for either of Equation (20) and
(19) across all possible $m\times m$ grids. In the case of summation, the
values over all possible $m\times m$ grids are computed and the sum is
returned, while the maximization case reports only the largest of these
values. The distributions of these statistics is computed empirically just as
for the MIC and Jiang et al. (2015). The language used, where aggregation is
compared to the Cramer-von Mises criterion and maximization to the Kolmogorov-
Smirnov test, is evocative of methods for testing empirical distribution
functions in Stephens (1974).

To implement these methods, all of Heller et al. (2016); Jiang et al. (2015);
and Reshef et al. (2011) utilize similar recursive algorithms. While Heller et
al. (2016) compute their summation statistic directly, as this reduces to a
counting problem, they borrow the procedure of Jiang et al. (2015) with a
different penalty for the maximization case. Reshef et al. (2011) use a
similar concept to guide their exploration of the space of grids with $IJ\leq
n^{0.6}$.

For a more detailed discussion of this recursive algorithm and its
applications, introduce the notation

$D(\mathbf{b}_{\mathbf{x}},\mathbf{b}_{\mathbf{y}})=\sum_{i=1}^{I}\sum_{j=1}^{J}d(o_{ij},e_{ij}),$
(21)

thereby expressing $D$ over the entire data set as the sum of $d$ applied to
each cell. Introduce the subscript notation

$\mathbf{b}_{\mathbf{x}[1:k]}={(b_{X}(x_{1}),b_{X}(x_{2}),\dots,b_{X}(x_{k}))}^{\mkern-1.5mu\mathsf{T}}$

to denote a binning on $I$ bins applied to the first $k$ elements of
$\mathbf{x}$. Consider the value of Equation (21) applied to the first $k$
observations of $\mathbf{x}$ and $\mathbf{y}$ with a given binning on
$\mathbf{y}$, written

$D_{k}(\cdot|\mathbf{b}_{\mathbf{y}}):\\{1,\dots,I\\}^{k}\mapsto\mathbb{R}$
(22)

for $I\leq k$ and a given binning $\mathbf{b}_{\mathbf{y}}$ on $J$ bins
defined for all of $\mathbf{y}$ (not just the first $k$ observations).

Equation (22) presents an important modification of Equation (21). By viewing
only the first $k$ elements with a _pre-specified_ $\mathbf{b}_{\mathbf{y}}$,
the problem of identifying a binning which optimizes $D_{k}$ is greatly
simplified compared to $D$. Rather than selecting among all possible grids,
only a small number at each step need to be considered. Suppose
$\mathbf{b}^{*}_{\mathbf{x}[1:k]}$ is the binning on the first $k$ elements of
$\mathbf{x}$ which maximizes $D_{k}$, define

$D^{*}_{k}=D_{k}(\mathbf{b}^{*}_{\mathbf{x}[1:k]}|\mathbf{b}_{\mathbf{y}}),$

as the maximal value of Equation (22). It is therefore the maximal value of
$D$ applied to the first $k$ observations of $\mathbf{x}$ and $\mathbf{y}$
given a known binning on $\mathbf{y}$.

Using this general notation, the recursive estimate $\widehat{D}_{k}^{*}$ used
by Reshef et al. (2011); Jiang et al. (2015); and Heller et al. (2016) is

$\widehat{D}^{*}_{k}=\max_{1\leq
i<k}\left[\widehat{D}^{*}_{i-1}+\sum_{j=1}^{J}d\left(\sum_{l=i}^{k}I_{\\{j\\}}\left(b_{Y}(y_{l})\right),\frac{k-i}{n}\sum_{l=1}^{n}I_{\\{j\\}}\left(b_{Y}(y_{l})\right)\right)\right].$
(23)

The edges which give these optimal binnings are given by replacing the $\max$
with an $\operatorname*{arg\,max}$. The arguments inside $d$ are the observed
counts of each $Y$ bin from $i$ to $k$ and the expected counts based on the
marginal distribution and the relative length of the interval from $i$ to $k$,
respectively.

For each $i\in\\{1,\dots,k-1\\}$, this estimate considers two parts. The first
part, $\widehat{D}_{i-1}^{*}$, is the previously computed maximal measure on
the first $i-1$ points, with the convention $\widehat{D}^{*}_{0}=0$. The
second part is a sum of $d$ over all $J$ bins of $\mathbf{b}_{\mathbf{y}}$
with observed counts given by the incidence of each of the $J$ bins in
observations $i$ to $k$ and the expected counts given by the marginal
distribution of $\mathbf{b}_{\mathbf{y}}$ multipled by the length of the
interval from $i$ to $k$.

Therefore, this algorithm chooses to add a bin edge at the $i$ which maximizes
Equation (21) _conditional on previous bin edges_. Additionally, this estimate
requires a _pre-specified_ $\mathbf{b}_{\mathbf{y}}$, and so it does not give
a global optimum over all grids. Indeed, it still requires some ad hoc choice
of binning on $\mathbf{y}$. Jiang et al. (2015) suggests using the slicing
methods of Jiang and Liu (2013) to bin $\mathbf{y}$ while Reshef et al. (2011)
suggests a simple equipartition on $\mathbf{y}$. The former chooses bins to
optimize the difference between the conditional and unconditional variances
while the latter is computationally easy to implement.

The main difference between the implementations in Jiang et al. (2015), Heller
et al. (2016), and Reshef et al. (2011) is the choice of $d$. Jiang et al.
(2015) takes a penalized measure

$d(o_{ij},e_{ij})=\frac{o_{ij}}{n}\log\left(\frac{o_{ij}}{e_{ij}}\right)-\frac{\lambda_{0}}{J}\log
n$ (24)

where $\lambda_{0}$ is a penalty parameter and the given
$\mathbf{b}_{\mathbf{y}}$ is a binning on $J$ bins. This is equivalent to a
particular prior being placed on the bin locations so that the number of bins
follows a Poisson distribution. Heller et al. (2016) instead take

$d(o_{ij},e_{ij})=\frac{o_{ij}}{n}\log\left(\frac{o_{ij}}{e_{ij}}\right)+\frac{\lambda_{0}}{J}\log{n-1\choose
k-1},$ (25)

where $k$ matches the index $\widehat{D}^{*}_{k}$. This modification is
equivalent to a uniform distribution over all possible marginal binnings.
Under both versions of the penalized $d$, splits are then considered at each
$x_{k}$ in turn and the bins which maximize the penalized score are used. In
contrast, Reshef et al. (2011) take the unpenalized

$d(o_{ij},e_{ij})=\frac{o_{ij}}{n}\log\left(\frac{o_{ij}}{e_{ij}}\right)$

and choose to restrict the search space to avoid creating bins which are too
small.

Attempts to expand these procedures to optimize over both margins have also
been made. Chen et al. (2016) suggest a modified version of the MIC which is
guided by the significance of a $\chi^{2}$ test to choose the optimal bins.
Cao et al. (2021) suggests a backwards merging algorithm on
$\mathbf{b}_{\mathbf{y}}$ to relax the perfect equipartition of Reshef et al.
(2011). In both cases, however, an equipartition is the starting point, and so
will inevitably impact the final binning.

Given that a conditional and recursive optimization procedure is already the
norm, it seems natural to use recursive binary splits to generate a binning as
in classification and regression trees and the tree-based binning of Rahman
(2018). Recursive binning has several advantages over marginal binning. For
one, recursive binning in these contexts can be set to optimizes the bin edge
choice over both dimensions simultaneously rather than each alone. It also
produces more flexible bin arrangements than marginal methods, which can only
produce bins aligned along both axes. Finally, recursive splits are adaptive
to patterns in the data which are hidden in projections on either axis,
suggesting that the method may detect patterns marginal bins miss.

The main problem recursive binning creates is the estimation of expected
counts in a bin, which can no longer proceed by taking the product of marginal
distributions. However, by first taking the ranks $\mathbf{r}(\mathbf{x})$ and
$\mathbf{r}(\mathbf{y})$, the expected count can be determined by the area of
the bin directly rather than a marginal product.

##  5 Recursive rank binning to measure association

Consider a pair of variables $X$ and $Y$ realized in a sample of $n$ paired
observations $\mathbf{x}$ and $\mathbf{y}$, as in Section 2. Define the
vectors of marginal ranks $\mathbf{s}=\mathbf{r}(\mathbf{x})$ and
$\mathbf{t}=\mathbf{r}(\mathbf{y})$ with the convention of random tie-breaking
for observations with the same rank to avoid ties. Converting to the ranks is
equivalent to taking empirical CDF transforms of $X$ and $Y$ so that the joint
distribution of $\mathbf{s}$ and $\mathbf{t}$ is the empirical copula of $X$
and $Y$ with a uniform joint distribution under independence. Specifically,
this means the expected number of points in a region of $\\{1,\dots,n\\}^{2}$
is equal to the region’s area divided by $n$ under the null hypothesis of no
dependence. This conversion allows much more flexible partitions to be
considered than marginal partitioning. Here, recursive binary splits are
proposed on $\mathbf{s}$ and $\mathbf{t}$ to take advantage of this
flexibility.

To sketch the algorithm, first consider the objects it acts upon. From the
perspective of this algorithm, a bin is a rectangular subspace of
$\\{1,\dots,n\\}^{2}$ which may contain some pairs from the paired vectors
$\mathbf{s}$ and $\mathbf{t}$. It is defined by its lower and upper bounds in
each dimension, and has implicit features such as its area, its depth, the
number of pairs it is expected to contain under independence, and the number
of pairs actually observed within its bounds. The depth of a bin can be
understood as the number of recursive calls required to produce it from the
initial state: a single bin with bounds of $0$ and $n$ in both dimensions that
contains every observation.

At each step, the algorithm is presented with a partition of
$\\{1,\dots,n\\}^{2}$ into a collection of bins resulting from the preceding
splits made by the algorithm. For each bin, the algorithm must choose whether
it should be split, and if so how it should be split. Under completely
unrestricted splitting, splits could be made through either of the two margins
along any horizontal or vertical line within the bin boundaries. Choosing
among the infinite possible splits is accomplished by finding the split
optimizing a score function reflecting the goal of the binning. Once a split
is chosen for every bin to be split, the algorithm proceeds recursively by
considering the resulting partition in the same way. Two choices made by the
analyst therefore dictate the final bins produced by the algorithm: the
_score_ used to choose splits and the _stop criteria_ which determine whether
a split is made at all.

Using the settings of the following section, this algorithm has a runtime
proportional to $n\log n$ if bin count or area limits are used as the stop
criteria. At each step, it searches through each point in each bin, meaning
all $n$ points are considered. In the worst case, this will only halve the bin
count and area at each step, and so $\log n$ splits are required. In any other
case, not all $n$ points will be considered at every depth. If small bins are
created early in the procedure, the points they contain will be ignored
thereafter.

###  5.1 Splitting bins

Many heuristics exist to choose splits, see Garcia et al. (2012) and Rahman
(2018) for surveys, but previous bin-based measures of association have
focused on the $\chi^{2}$ statistic against independence and the mutual
information (MI) (Reshef et al., 2011; Jiang et al., 2015; Chen et al., 2016;
Heller et al., 2016; Cao et al., 2021). As they are designed to measure the
discrepancy of observed distributions from expected distributions with minimal
assumptions, both are natural choices to measure the dependence present in a
sample.

However, these previous works choose only marginal splits on both dimensions,
and so must be adapted to the recursive binning framework of Rahman (2018) by
defining local versions of both, the _chi score_ and _mi score_ , to select
the optimal split within a bin. It will be proven that, for either score, the
optimal split occurs at the coordinate of a point within a bin. To distinguish
between the scores and the final statistics, _chi_ and _mi_ will be used
exclusively to refer to the scores computed to determine splitting and
$\chi^{2}$ and MI will be used to refer to the final statistics computed over
all bins.

To define these local scores recall Equation (21),

$D(\mathbf{b}_{\mathbf{x}},\mathbf{b}_{\mathbf{y}})=\sum_{i=1}^{I}\sum_{j=1}^{J}d(o_{ij},e_{ij}),$

which expresses a functional measuring statistical dependence over $X$ and $Y$
as the sum of a function evaluated over the expected ($e_{ij}$) and observed
($o_{ij}$) number of points in the partitions created by marginal binnings
$\mathbf{b}_{\mathbf{x}}$ and $\mathbf{b}_{\mathbf{x}}$. By adapting Equation
(21) for the $\chi^{2}$ and MI statistics to the recursive binning framework,
the chi and mi scores are implied by the form of $d(o_{ij},e_{ij})$. As the
bins produced by recursive binary splits are not defined by independent
marginal binnings $\mathbf{b}_{\mathbf{x}}$ and $\mathbf{b}_{\mathbf{y}}$,
they do not have obvious $i,j$ indices. Instead, assume the total number of
bins is $n_{bin}$ and (arbitrarily) index the bins by
$i\in\\{1,\dots,n_{bin}\\}$. The arguments must also be changed to
$\mathbf{s}$ and $\mathbf{t}$ to reflect the absence of marginal bins.
Together, this gives the modified expression

$D(\mathbf{s},\mathbf{t})=\sum_{i=1}^{n_{bin}}d(o_{i},e_{i}),$ (26)

where $o_{i}$ is the number of observations within the $i^{\text{th}}$ bin and
$e_{i}$ is the number expected assuming independence. Letting the area of the
$i^{\text{th}}$ bin be $a_{i}$, this is given by $e_{i}=a_{i}/n$. The chi
score of each bin is then

$d(o_{i},e_{i})=\text{chi}(o_{i},e_{i})=\frac{(o_{i}-e_{i})^{2}}{e_{i}}$ (27)

and the mi score is

$d(o_{i},e_{i})=\text{mi}(o_{i},e_{i})=\frac{o_{i}}{n}\log\frac{o_{i}}{e_{i}}.$
(28)

Either of these scores can be maximized in the proposed recursive binning
algorithm to determine the split coordinate.

####  5.1.1 Maximizing scores

In more detail, denote the $o_{i}$ pairs of ranks in bin $i$ as
$\\{(s_{i1},t_{i1}),(s_{i2},t_{i2}),\dots,(s_{io_{i}},t_{io_{i}})\\}$ and its
$s$ bounds $(l_{s},u_{s}]$ and $t$ bounds $(l_{t},u_{t}]$. Bin $i$ can be
split either by a vertical line at $c_{s}\in(l_{s},u_{s})$ or a horizontal
line at $c_{t}\in(l_{t},u_{t})$ resulting in two new bins with two new chi or
mi scores. Denote the observed and expected values for the bin above $c_{s}$
as $o_{i+}(c_{s})$ and $e_{i+}(c_{s})$ respectively (analogously, those above
$c_{t}$ as $o_{i+}(c_{t})$ and $o_{i+}(c_{t})$), and use the subscript $i-$ in
the same way to indicate the new bin below the split. A split at $c$ changes
the total score measured by $d(\cdot,\cdot)$ for the region
$(l_{s},u_{s}]\times(l_{t},u_{t}]$ by

$\delta_{i}(c,d)=d\big{(}o_{i+}(c),e_{i+}(c)\big{)}+d\big{(}o_{i-}(c),e_{i-}(c)\big{)}-d(o_{i},e_{i}),$
(29)

and so the maximizing split coordinate along a given dimension is

$c^{*}=\operatorname*{arg\,max}_{c}\delta_{i}(c,d)=\operatorname*{arg\,max}_{c}\Big{[}d\big{(}o_{i+}(c),e_{i+}(c)\big{)}+d\big{(}o_{i-}(c),e_{i-}(c)\big{)}\Big{]}.$

Though $e_{i+}(c)$ and $e_{i-}(c)$ vary continuously in the split coordinate
$c$, both of $o_{i+}(c)$ and $o_{i-}(c)$ change only when $c$ corresponds with
the coordinate of a point contained in bin $i$, in other words when
$c_{s}\in\\{s_{i1},s_{i2},\dots,s_{io_{i}}\\}$ or
$c_{t}\in\\{t_{i1},t_{i2},\dots,t_{io_{i}}\\}$. This has important
consequences to selecting splits for both the mi and chi scores.

######  Proposition 1 (The split maximizing the chi score occurs at the
coordinate of a point within a bin.).

The split coordinate $c$ which maximizes $\delta_{i}(c,\text{chi})$ is the
coordinate of one of the points within the bin.

###### Proof.

Without loss of generality, consider a split at $c_{s}\in(s_{ij},s_{i,j+1})$
between the $j$ and $j+1$ horizontal coordinates in bin $i$. As the split is
made between point coordinates $j$ and $j+1$, the observed number of points
above is $o_{i+}=o_{i}-j$ and the observed number below is $o_{i-}=j$ by
definition. Substituting Equation (27) into Equation (29) gives

$\displaystyle\delta_{i}(c_{s},\text{chi})$ $\displaystyle=$
$\displaystyle\frac{(o_{i+}-e_{i+}(c_{s}))^{2}}{e_{i+}(c_{s})}+\frac{(o_{i-}-e_{i-}(c_{s}))^{2}}{e_{i-}(c_{s})}-\frac{(o_{i}-e_{i})^{2}}{e_{i}}$
$\displaystyle=$
$\displaystyle\frac{o_{i+}^{2}}{e_{i+}(c_{s})}-2o_{i+}+e_{i+}(c_{s})+\frac{o_{i-}^{2}}{e_{i-}(c_{s})}-2o_{i-}+e_{i-}(c_{s})-\frac{o_{i}^{2}}{e_{i}}+2o_{i}-e_{i},$
$\displaystyle=$
$\displaystyle\frac{(o_{i}-j)^{2}}{e_{i+}(c_{s})}+\frac{j^{2}}{e_{i-}(c_{s})}-\frac{o_{i}^{2}}{e_{i}}$

as $e_{i+}(c)+e_{i-}(c)=e_{i}$ and $o_{i+}+o_{i-}=o_{i}$.

The derivative of $\delta_{i}(c_{s},\text{chi})$ with respect to $c_{s}$ can
be determined with the help of a few observations. First,
$e_{i+}(c_{s})=(u_{s}-c_{s})(u_{t}-l_{t})/n$ and
$e_{i-}(c_{s})=(c_{s}-l_{s})(u_{t}-l_{t})/$ so that
$\frac{d}{dc}e_{i+}(c_{s})=-\frac{d}{dc}e_{i-}(c_{s})$. Furthermore, both of
the observed counts and the final term in the sum are constant in $c_{s}$ for
the restricted range considered, and so all have a derivative of zero.
Therefore

$\frac{d}{dc}\delta_{i}(c_{s},\text{chi})=-\frac{o_{i+}^{2}}{e_{i+}^{2}(c_{s})}\frac{d}{dc}e_{i+}(c_{s})-\frac{o_{i-}^{2}}{e_{i-}^{2}(c_{s})}\frac{d}{dc}e_{i-}(c_{s})=\left[\frac{o_{i+}^{2}}{e_{i+}^{2}(c_{s})}-\frac{o_{i-}^{2}}{e_{i-}^{2}(c_{s})}\right]\frac{d}{dc}e_{i-}(c_{s})$

which simplifies to

$\frac{d}{dc}\delta_{i}(c_{s},\text{chi})=\frac{n}{u_{t}-l_{t}}\left[\frac{o_{i+}^{2}}{(u_{s}-c_{s})^{2}}-\frac{o_{i-}^{2}}{(c_{s}-l_{s})^{2}}\right]$

with a second derivative

$\frac{d^{2}}{dc^{2}}\delta_{i}(c_{s},\text{chi})=\frac{2n}{u_{t}-l_{t}}\left[\frac{o_{i+}^{2}}{(u_{s}-c_{s})^{3}}+\frac{o_{i-}^{2}}{(c_{s}-l_{s})^{3}}\right].$

As the split coordinate $c_{s}\in(s_{ij},s_{i,j+1})$ is necessarily restricted
to within the bin bounds $(l_{s},u_{s}]$, both $c_{s}-l_{s}>0$ and
$u_{s}-c_{s}>0$ and so $\frac{d^{2}}{dc^{2}}\delta_{i}(c_{s},\text{chi})>0$.
This implies that $\delta_{i}(c_{s},\text{chi})$ is concave up between the
horizontal coordinates of the points within a bin so that any optimum within
these bounds must be minimum. As these are continuous functions, this means
maximum must occur at one of the boundaries of the interval
$(s_{ij},s_{i,j+1})$. As the index $j$ was chosen arbitrarily, this same
argument holds for every interval and so the global maximum must occur at one
of these boundaries. These boundaries are defined by the locations of the
points contained within the bin, so the maximal split must occur at the
coordinate of a point within the bin. The same argument holds identically for
the vertical coordinates, a fact easily seen by switching the subscripts. ∎

######  Proposition 2 (The split maximizing the mi score occurs at the
coordinate of a point within a bin.).

The split coordinate $c$ which maximizes $\delta_{i}(c,\text{mi})$ is the
coordinate of one of the points within the bin.

###### Proof.

Using the same argument as for Proposition 1, this is true as long as the
second derivative of $\delta_{i}(c_{s},\text{mi})$ with respect to $c_{s}$ is
non-negative between adjacent observation coordinates. Substituting Equation
(28) into Equation (29) gives

$\delta_{i}(c_{s},\text{mi})=\frac{o_{i+}}{n}\log\frac{o_{i+}}{e_{i+}(c_{s})}+\frac{o_{i-}}{n}\log\frac{o_{i-}}{e_{i-}(c_{s})}-\frac{o_{i}}{n}\log\frac{o_{i}}{e_{i}}$

with a first derivative

$\frac{d}{dc}\delta_{i}(c_{s},\text{mi})=\frac{o_{i+}}{n(u_{s}-c_{s})}-\frac{o_{i-}}{n(c_{s}-l_{s})}$

and a second derivative

$\frac{d^{2}}{dc^{2}}\delta_{i}(c_{s},\text{mi})=\frac{o_{i+}}{n(u_{s}-c_{s})^{2}}+\frac{o_{i-}}{n(c_{s}-l_{s})^{2}}\geq
0.$

Therefore, the value of $c_{s}$ which maximizes $\delta_{i}(c_{s},\text{mi})$
must correspond to the coordinate of a point within the bin. ∎

This means that splits only need to be considered at the points in
$\\{(s_{i1},t_{i1}),(s_{i2},t_{i2}),\dots,$ $(s_{io_{i}},t_{io_{i}})\\}$ to
maximize the chi and mi scores, rather than considering the continuum of
potential splits in $(l_{s},u_{s}]\times(l_{t},u_{t}]$.

####  5.1.2 Empty bins

If splitting occurs at the coordinates of observed points within bins, the
observation at the split has ambiguous bin membership. Taking the convention
that these observations are counted in the bin below the split, which is
consistent with the initial bin bounds $[0,n]$ in both dimensions, it is
impossible for this algorithm to create empty bins below the observed points.
This is despite the potential utility of empty bins when detecting
association, as large regions without any observations in the rank space are a
strong indication of departures from independence.999The ability to create
empty bins directly is an advantage of the recursive binning algorithm over
marginal methods, which cannot doso. Empty regions marginally only reflect the
marginal distribution of a variable, and converting to the ranks presents
margins without any gaps for both variables.

To remedy this and allow empty bins to be created, a potential split
coordinate below the smallest observations horizontally and vertically is
added, denote these pseudo-observations as $s_{i(1)}-1$ and $t_{i(1)}-1$.
Finally, take the convention that a point is included in the lower bin when a
split occurs at one of its coordinates. Along $s$ this gives the maximizing
split coordinate

$c_{s}^{*}=\operatorname*{arg\,max}_{c_{s}\in\\{s_{min}-1,s_{i1},\dots,s_{io_{i}}\\}}\delta_{i}(c_{s},d)$

and similarly

$c_{t}^{*}=\operatorname*{arg\,max}_{c_{t}\in\\{t_{min}-1,t_{i1},\dots,t_{io_{i}}\\}}\delta_{i}(c_{t},d).$

Of these two maximizing splits, that giving the greater $\delta_{i}(c,d)$ is
chosen to split bin $i$.

####  5.1.3 Controlling minimum bin size

Though not strictly necessary, one may want to control the minimum bin size
produced by splits. This requires some balance to be struck between selecting
the maximal split and keeping bins at a particular size. This is relevant, for
example, if the $\chi^{2}$ statistic is applied to the final bins. Supposing
$n_{bin}$ bins are created by the algorithm, the rank space $[0,n]^{2}$ with a
presumed uniform distribution has been partitioned into $n_{bin}$ mutually
exclusive categories constrained only by the restriction that

$\sum_{i=1}^{n_{bin}}o_{i}=n.$

This setting is similar to the circumstance originally considered by Pearson
in his proposal of the $\chi^{2}$ test, and so it is natural to assume a
$\chi^{2}_{n_{bin}-1}$ distribution for the $\chi^{2}$ statistic applied to
this bins.

Care must be taken with bin size in order to apply this result, however. The
$\chi^{2}_{n_{bin}-1}$ distribution is only asymptotically valid for the
$\chi^{2}$ statistic applied to this data and the fit is better the larger the
expected counts in each bin. Therefore, the $\chi^{2}$ distribution is
typically only applied to the $\chi^{2}$ statistic when the expected number of
points in each partition is greater than or equal to five (Cochran, 1952).
This motivates a floor on the bin size at an expected value of five.

The preceding maximization logic provides no such guarantees, and so
restrictions on the candidate splits must be introduced to control the minimal
bin size. Rather than change the score function, this can be accomplished by
changing the $\delta_{i}(c,d)$ function used to evaluate splits. Consider the
modified function

$\delta_{i}^{\prime}(c,\text{chi},z)=I_{[z,\infty)^{2}}\left({(e_{i+}(c),e_{i-}(c))}^{\mkern-1.5mu\mathsf{T}}\right)\delta_{i}(c,\text{chi})$
(30)

that forces the change of score to be zero if either $e_{i+}(c)<z$ or
$e_{i-}(c)<z$, where $I_{A}\left(x\right)$ is the indicator function of $x\in
A$. This change works because $\delta_{i}(c,\text{chi})\geq 0$, as

$\displaystyle\frac{(o_{i}-e_{i})^{2}}{e_{i}}$ $\displaystyle=$
$\displaystyle\frac{\big{(}[o_{i+}(c)-e_{i+}(c)]+[o_{i-}(c)-e_{i-}(c)]\big{)}^{2}}{e_{i}}$
$\displaystyle\leq$
$\displaystyle\frac{[o_{i+}(c)-e_{i+}(c)]^{2}}{e_{i}}+\frac{[o_{i-}(c)-e_{i-}(c)]^{2}}{e_{i}}$
$\displaystyle\leq$
$\displaystyle\frac{[o_{i+}(c)-e_{i+}(c)]^{2}}{e_{i+}(c)}+\frac{[o_{i-}(c)-e_{i-}(c)]^{2}}{e_{i-}(c)}$

by the triangle inequality and because $e_{i+}(c)\leq e_{i}$ and
$e_{i-}(c)\leq e_{i}$. Therefore, splits producing bins that are too small
will give scores less than or equal to the scores produced by all other
splits. Taking $z=5$ restricts bin splits to follow standard practice.

Controlling the bin size for $\delta_{i}(c,\text{mi})$ requires a different
convention, as

$\displaystyle\frac{o_{i}}{n}\log\frac{o_{i}}{e_{i}}$ $\displaystyle=$
$\displaystyle\frac{o_{i}}{n}\log\left[\frac{o_{i+}(c)}{o_{i}}\frac{o_{i}}{e_{i}}+\frac{o_{i-}(c)}{o_{i}}\frac{o_{i}}{e_{i}}\right]$
$\displaystyle\leq$
$\displaystyle\frac{o_{i}}{n}\left[\frac{o_{i+}(c)}{o_{i}}\log\frac{o_{i}}{e_{i}}+\frac{o_{i-}(c)}{o_{i}}\log\frac{o_{i}}{e_{i}}\right]$
$\displaystyle\leq$
$\displaystyle\frac{o_{i+}(c)}{n}\log\frac{o_{i-}(c)}{e_{i}}+\frac{o_{i+}(c)}{n}\log\frac{o_{i+}(c)}{e_{i}}+\frac{o_{i-}(c)}{n}\log\frac{o_{i-}(c)}{e_{i}}+\frac{o_{i-}(c)}{n}\log\frac{o_{i+}(c)}{e_{i}}$

which may be greater than
$\frac{o_{i+}(c)}{n}\log\frac{o_{i+}(c)}{e_{i+}(c)}+\frac{o_{i-}(c)}{n}\log\frac{o_{i-}(c)}{e_{i-}(c)}$.
However, noting that $\text{mi}(o_{i},e_{i})$ and $\text{chi}(o_{i},e_{i})$
are both independent of the split line $c$, maximization of
$\delta_{i}(c,\text{chi})$ and $\delta_{i}(c,\text{chi})$ depends only on the
post-split scores
$\text{mi}(o_{i+}(c),e_{i+}(c))+\text{mi}(o_{i+}(c),e_{i+}(c))\geq 0$ and
$\text{chi}(o_{i+}(c),e_{i+}(c))+\text{chi}(o_{i+}(c),e_{i+}(c))\geq 0$.
Splitting to control the minimum bin size can therefore proceed with
indicators by taking the larger of

$c_{s}^{*}=\operatorname*{arg\,max}_{c_{s}\in\\{s_{min}-1,s_{i1},\dots,s_{io_{i}}\\}}I_{[z,\infty)^{2}}\left({(e_{i+}(c_{s}),e_{i-}(c_{s}))}^{\mkern-1.5mu\mathsf{T}}\right)\big{[}d(o_{i+}(c_{s}),e_{i+}(c_{s}))+d(o_{i+}(c_{s}),e_{i+}(c_{s}))\big{]}$

and

$c_{t}^{*}=\operatorname*{arg\,max}_{c_{t}\in\\{t_{min}-1,t_{i1},\dots,t_{io_{i}}\\}}I_{[z,\infty)^{2}}\left({(e_{i+}(c_{t}),e_{i-}(c_{t}))}^{\mkern-1.5mu\mathsf{T}}\right)\big{[}d(o_{i+}(c_{t}),e_{i+}(c_{t}))+d(o_{i+}(c_{t}),e_{i+}(c_{t}))\big{]}$

when $d\in\\{\text{chi},\text{mi}\\}$. In the case where all splits are tied
on both margins, the bin is halved on a random margin at

$\text{ceiling}\left(\frac{l+u}{2}\right)$

regardless of the distribution of points within the bin.101010Choosing a
random split is another obvious tie-breaking convention, but can force splits
which violate the minimum size restriction.

###  5.2 Stop criteria

At each recursive step, a bin is split only if it fails to meet a set of stop
criteria. These can include the bin area, number of points in the bin, and the
depth of the bin111111Where “depth” is the number of recursive binary splits
needed to create the bin from the original data.. If, for example, the
criteria are a depth of 5 or $n_{i}<10$, a bin with a depth of 5 or a bin with
10 or fewer points will not be split. Any bin which does not satisfy the stop
criteria is split in two and the splitting algorithm is again called on both
of the resulting bins.

If bin size is restricted such that $e_{i}\geq z$ $\forall
i\in\\{1,\dots,n_{bin}\\}$, the corresponding stop criterion to prevent the
creation of bins below this minimum size is $e_{i}<2z$. When $e_{i}<2z$, any
split will produce at least one bin with an expected count smaller than $z$.
Even if there is no restriction $e_{i}\geq z$ when splitting, it is advisable
to incorporate one in the stop criteria to limit the creation of very small
bins.

In all of the following examples, two stop criteria were held constant.
Splitting was always stopped when $n_{i}=0$ or $e_{i}\leq 10$. To explore the
performance of the algorithm over successive splits, the depth criterion was
varied from 2 to 10. For smaller sample sizes, the maximal 1024 bins implied
by the depth limit of 10 was never achieved due to the stop criteria for area
and empty bins.

##  6 An iterative version

Though it was conceived recursively, an iterative implementation of the
algorithm outlined in Section 5 is detailed here. First, consider the outer
wrapper function to oversee arbitrary splitting and stopping, presented in
Algorithm 1. This wrapper acts on a list of bin objects, each with elements:

  1. bounds:

named list of upper and lower bounds on $s$ and $t$

  2. s:

marginal ranks on $x$ of observations in the bin

  3. t:

marginal ranks on $y$ of observations in the bin

  4. e:

expected number of points in the bin

  5. depth:

number of recursive splits required to create the bin.

While any bins in a list fail to satisfy the stop criteria, this wrapper calls
the splitting function on those bins and combines the resulting new bins with
those already stopped. The new bins are then checked against the stop
criteria. To initialize, all observations are placed in a bin with bounds of
$(0,n]$ in both dimensions.

Algorithm 1 Iterative binning wrapper

Input

$\mathbf{x}$ \- vector of observed values in $x$

$\mathbf{y}$ \- vector of observed values in $y$, paired with $x$

stopper \- function which tests a list of bins against the stop criteria

5: splitter \- function which performs a binary split on a bin

init \- (possibly) different splitting function applied to the initial bin
with all observations

function binner($\mathbf{x}$, $\mathbf{y}$, stopper, splitter, init)

$n$ $\longleftarrow$ length($\mathbf{x}$) $\triangleright$ Compute
preliminaries

10: $\mathbf{s}$ $\longleftarrow$ rank($\mathbf{x}$)

$\mathbf{t}$ $\longleftarrow$ rank($\mathbf{y}$)

iniBin $\longleftarrow$ makeBin(s = $\mathbf{s}$, t = $\mathbf{t}$, bounds =
list(s = ($0,n$), t = ($0,n$)), e = $n$, depth = $0$) $\triangleright$
Construct initial bin

binList $\longleftarrow$ init(iniBin) $\triangleright$ Initialize bin list

stopStat $\longleftarrow$ stopper(binList) $\triangleright$ Initial stop check

15: while any stopStat are FALSE do $\triangleright$ Continue as long as bins
can be split

oldBins $\longleftarrow$ binList[stopStat] $\triangleright$ Stopped bins

oldStop $\longleftarrow$ stopStat[stopStat] $\triangleright$ All TRUE

newBins $\longleftarrow$ {} $\triangleright$ Variable to store splitting
results

for $bin$ in binList[!stopStat] do

20: append splitter($bin$) to newBins $\triangleright$ Add split results

end for

newStop $\longleftarrow$ stopper(newBins) $\triangleright$ Check stop criteria
on new bins

binList $\longleftarrow$ append newBins to oldBins

stopStat $\longleftarrow$ append newStop to oldStop

25: end while

return binList

end function

Of course, the splitting and stopping logic described in Sections 5.1 and 5.2
are at the heart of this algorithm. Algorithm 2 presents pseudo-code for the
splitting logic. The stopping logic is not given a pseudo-code version, as it
only involves computing and checking numerous properties of each bin against a
pre-determined set of thresholds. Section 7 discusses a particular
implementation of this in the R programming language.121212In this
implementation, there is the ability to specify a custom initial split
function with the argument init. Choosing init to randomly halve a bin along
one margin coincides with the splitting logic described earlier, as every
split produces a score of zero for the initial uniform data.

Algorithm 2 Score maximizing splitter

Input

bin \- the bin object to be split with elements s, t, bounds, e, and depth

scorer \- function accepting a vector of coordinates in increasing order and
the expected number of points in bin and returning a vector of scores
corresponding to splits at each internal coordinate

5:function maxScoreSplit(bin, scorer)

get s, t, s.bounds, t.bounds, e, depth from bin

sort s, t in increasing order to give $sSrt,tSrt$

$\mathbf{c}_{s}$ $\longleftarrow$ append ( $sSrt$[1]-1, $sSrt$ )

$\mathbf{c}_{t}$ $\longleftarrow$ append ( $tSrt$[1]-1, $tSrt$ )
$\triangleright$ Add a split coordinate below all points

10: $\mathbf{d}_{s}$ $\longleftarrow$ scorer( append ( s.bounds[1],
$\mathbf{c}_{s}$, s.bounds[2] ), e )

$\mathbf{d}_{t}$ $\longleftarrow$ scorer( append ( t.bounds[1],
$\mathbf{c}_{t}$, t.bounds[2] ), e )

$s_{Max}$, $t_{Max}$ $\longleftarrow$ the indices of the maxima of
$\mathbf{d}_{s}$, $\mathbf{d}_{t}$

if all $\mathbf{d}_{s}$ = $\mathbf{d}_{s}$[$s_{Max}$] AND all $\mathbf{d}_{t}$
= $\mathbf{d}_{t}$[$t_{Max}$] then

if $\mathbf{d}_{s}$[$s_{Max}$] $>$ $\mathbf{d}_{t}$[$t_{Max}$] then

15: split bin on $s$ at ceiling(mean(s.bounds))

else if $\mathbf{d}_{s}$[$s_{Max}$] $<$ $\mathbf{d}_{t}$[$t_{Max}$] then

split bin on $t$ at ceiling(mean(t.bounds))

else

split bin randomly on $s$ or $t$ at ceiling(mean(s.bounds)) or
ceiling(mean(t.bounds))

20: end if

else if $\mathbf{d}_{s}$[$s_{Max}$] $\geq$ $\mathbf{d}_{t}$[$t_{Max}$] then
$\triangleright$ Ties go to $s$

split bin on $s$ at $\mathbf{c}_{s}$[$s_{Max}$]

else

split bin on $t$ at $\mathbf{c}_{t}$[$t_{Max}$]

25: end if

return upper and lower split of bin

end function

Finally, Algorithm 3 gives the pseudocode for a chi scoring function with
limited minimal bin size. As written, this function would be provided as the
score in the maxScoreSplit function of Algorithm 2.

Algorithm 3 Marginal $\chi^{2}$ scores

Input

$\mathbf{x}$ \- a numeric vector of potential split coordinates in increasing
order

$e$ \- a numeric value giving the expected number points in a bin

$e_{min}$ \- the minimum expected number allowed for a split

5:

function chiScoring($\mathbf{x}$, $e$, $e_{min}$)

$n\longleftarrow$ length($\mathbf{x}$)

$cumulative\longleftarrow\mathbf{x}$[2] $-\mathbf{x}$[1] $\triangleright$
Initialize cumulative length

$density\longleftarrow e/($max$(\mathbf{x})-\mathbf{x}$[1]$)$ $\triangleright$
Density under uniformity

10: scores $\longleftarrow$ {} $\triangleright$ Initialize storage

for $i=2$ to $n-1$ do

$e_{i}\longleftarrow cumulative*density$ $\triangleright$ Expected count below
$i$

$o_{i}\longleftarrow i-2$ $\triangleright$ Observed count ignores bounds,
pseudo-point

if $e_{i}\geq e_{min}$ AND $e-e_{i}\geq e_{min}$ then

15: append
$\frac{(o_{i}-e_{i})^{2}}{e_{i}}+\frac{(n-3-o_{i}-e+e_{i})^{2}}{e-e_{i}}$ to
scores $\triangleright$ chi score of candidate split

else

append $0$ to scores

end if

$cumulative\longleftarrow cumulative+\mathbf{x}[i+1]-\mathbf{x}[i]$
$\triangleright$ Update length

20: end for

return $scores$

end function

Note that this algorithm has been written for the chi score but provides a
framework for any scoring function. The specifics of line 15 can be replaced
with any objective function, for example the mi score of Equation (28) with
$\frac{o_{i}}{n-3}\log\frac{o_{i}}{e_{i}}+\frac{n-3-o_{i}}{n-3}\log\frac{n-3-o}{e-e_{i}}$
to maximize the marginal Kullback-Liebler divergence from uniformity.

This framework also supports random binning. Rather than computing a function
that compares the observed distribution to a uniform one about a split, the
scores can be replaced with independent and identically distributed $U[0,1]$
realizations. The coordinate of the maximum score value will then be uniformly
distributed along each potential split coordinate and each margin, creating a
fully random recursive binning procedure. Additionally, as these random
realizations are greater than or equal to zero, adopting the indicator
multiplication of Equation (30) controls the bin size under this form of
random splitting.

##  7 The AssocBin package

Algorithms 1, 2, and 3 are implemented in the R package AssocBin for the mi
score, chi score, and random splitting score. The core functions are:

makeCriteria:

a function which captures its arguments and appends them into a single logical
expression

stopper:

a function which accepts a list of bins and a logical expression and evaluates
the expression within each bin using R’s lexical scoping

binner:

the wrapper function described in Algorithm 1 which accepts integer vectors x
and y; a stopper function which accepts a list of bins and returns a logical
vector; a splitter function which accepts a single bin and returns a pair of
bins partitioning the input bin; an init function which splits the initial bin
containing all points; and (optionally) additional arguments to pass to
internal function calls

chiScores, miScores, randScores:

functions which implement Algorithm 3 with line 15 replaced by the
corresponding score

maxScoreSplit:

the splitting function described in Algorithm 2 which accepts a bin and scorer
function in addition to ties and pickMax which allow for custom specification
of how to handle score ties and choose a maximum

splitX, splitY:

functions which accept a bin to be split, a numeric bd giving the coordinate
of the split, and the indices of values above and below bd and return two bins
resulting from a split of bin at bd along the corresponding margin

halfCutTie:

the tie-breaking logic described in Algorithm 2

By writing binner with fully modular components, the score function, splitting
logic, stop criteria, and ties can be modified to suit the preference of a
user. As these are all imagined as single argument functions within binner,
this requires that certain functions be defined in closures before use. To use
stopper, for example, the stop criteria returned by makeCriteria must be
passed as an argument to stopper within another function which returns a
single-argument function. Similarly, maxScoreSplit must have its scoring
function and minimum expected count set in a closure that then accepts a
single argument. Though this requires some extra set up, it limits the
arguments of binner and forces the user to consider these choices
intentionally in advance. The demos included in the package provide numerous
examples of this set up.

Additional helpers visualize and summarize the results. The binChi function
computes the $\chi^{2}$ statistic over a list of bins returned by binner and
the plotBinning function plots a bin list and scaterplot with optional fill.
Two additional functions, depthFill and residualFill, create gradient fills to
communicate depth or residual magnitude based on colour range and residual
function arguments. Such visualizations not only give insight into what region
of the data the algorithm deems most important, but are also useful as a
summarized display of the data. Rather than a scatterplot with potentially
thousands of points, these visualization display a handful of coloured regions
which can be read at a glance. The full package can be found on the author’s
GitHub.

##  8 Using the algorithm

Some results gained from the implementation of this method in R are presented
here. First, a series of visualizations of the algorithm on independent data
and strongly associated data are given to demonstrate how it works step-by-
step. Then, as both Heller et al. (2016) and Reshef et al. (2011) utilize
simulated independent $X$ and $Y$ to generate the $p$-values of their methods,
the null distribution of the recursive binning measure under different
splitting rules is explored. This is followed by the application of the method
to simulated data patterns from Newton (2009) and Liu et al. (2018). Finally,
a real data set is examined.

###  8.1 Simple examples

To illustrate the behaviour of the maximum splitting algorithm from Section 6,
consider applying it to two extremes: random independent data and data in
perfect agreement. Scatterplots of both of these simulated data sets with
$n=1,000$ are shown in Figure 2.

|
---|---
(a) | (b)
Figure 2: The (a) simulated random data and (b) perfect rank agreement data
used to illustrate the flow of the algorithm.

As the ranks of both data sets have no gaps, the initial step of the algorithm
will not find a marginal maximal split. Indeed, for any split at a point the
number of points on either side will match expectation exactly. Therefore, the
algorithm begins by splitting the initial bin in half by adding a vertical
edge at 500 in both. This gives the bins seen in Figure 3.

|
---|---
(a) | (b)
Figure 3: The (a) simulated random data and (b) perfect rank agreement data
with the first split indicated.

Once so halved, the ranks are no longer necessarily uniform within the bins.
Marginal gaps are introduced in both due to the split. Therefore, the next
split is not made at the halfway point, but is instead made to optimize the
deviation of counts from uniformity as measured by the score function. While
identifying the location of these splits on either side in the random uniform
data is difficult, the expected split in the line data would be another
halving of the bins, as this produces empty corner bins which should each
contain a quarter of the observations. Indeed, this is the result seen in
Figure 4.

|
---|---
(a) | (b)
Figure 4: The (a) simulated random data and (b) perfect rank agreement data
with the first two splits indicated, dividing both data sets into four bins.

Such a step-by-step demonstration could be continued, to the tedium of the
reader, but an easier visualization can summarize such stepwise inspection.
For each bin, the implementation described in Section 7 keeps track of the bin
depth as well as its other features. More information about the algorithm’s
path can be gleaned by simply running the algorithm and shading the bins
according to their depth. Setting a maximum depth of six and shading in this
way produces Figure 5.

|
---|---
(a) | (b)
Figure 5: The (a) simulated random data and (b) perfect rank agreement data
split to a maximum depth of six, with bins shaded darker according to their
depth. The algorithm splits to the maximum depth along the line and not
elsewhere while the in the random data continues splitting in every location.

There are stark differences between the depth patterns for these two data
sets, as might be expected. While the uniform random data continues splitting
every bin somewhat haphazardly with the exception of small early splits that
meet the stop criteria (the expected number of observations being 10 or fewer
or no observations within the bin), the line data causes a very particular
pattern of depths to emerge. It seems to chase the linear pattern by
introducing many more splits along its length, in a sense introducing a
greater density of bins to the regions with greater density of points.

By chasing these regions, the algorithm produces a striking pattern in the
residuals of its final bins. This pattern occurs in, for example, the Pearson
residuals of bins with expected counts $e_{i}$ and observed counts $o_{i}$,

$r_{i}\text{sign}(o_{i}-e_{i})\sqrt{\frac{(o_{i}-e_{i})^{2}}{e_{i}}}.$

Large positive residuals occur all along the line and large negative ones
occur away from it. Figure 6 plots these residuals using hue to convey their
sign, with the convention that blue represents negative values and red
represents positive ones, and saturation to convey their magnitude, with
darker saturation indicating a larger magnitude.

|
---|---
(a) | (b)
Figure 6: (a) Simulated random data and (b) simulated data with perfect rank
agreement split to a maximum depth of six, with bins shaded according to their
Pearson residuals. Negative residuals are shaded blue while positive residuals
are shaded red, with darker shading indicating a larger residual magnitude.
Large positive residuals occur along the line and large negative ones occur
elsewhere.

While bins in the random data do not display pronounced shading, indicating
small residuals, the line pattern has large positive residuals all along its
length and large negative residuals elsewhere. In particular, note the deep
blue shading in upper left and bottom right quadrants. Even without the
points, the structure of the data is easily discerned from the shading of
these bins alone.

In this example, taking a sum of the squared Pearson residuals to get the
$\chi^{2}$ statistic gives values of 87.8 for the random data and 7000 for the
line data. Of course, as the bins produced here do not follow a regular grid
pattern and are generated by maximizing a score analogous to the $\chi^{2}$,
the distribution of the statistic warrants an investigation under
independence. As these bins are produced using optimization at each step, it
is not clear that $n_{bin}$ bins would have a $\chi^{2}$ distribution with
$n_{bin}-1$ degrees of freedom as is usually the case, even with the
controlled minimum bin size.

###  8.2 The null distribution

Even presuming the $\chi^{2}$ distribution applied to the $\chi^{2}$ statistic
on these bins, the degrees of freedom would be given by $n_{bin}-1$. As
$n_{bin}$ is a result of the stop criteria and splitting logic, the null
distribution will depend on both of these. The nature of its dependence on the
splitting logic is explored here. Three different splitting rules are
considered: random splits, chi score maximizing splits, and mi score
maximizing splits. All will lead to different final bins when applied to the
same data, and therefore different statistic values. To give the
$\chi^{2}_{n_{bin}}$ distribution a better chance of fitting, all were
restricted so that bins with expected counts less than 5 were not created.

Consequently, the null distribution is investigated by generating 10,000
independent bivariate data sets $\mathbf{x},\mathbf{y}\in\mathbb{R}^{10,000}$
by randomly and independently shuffling the vector of ranks
$\\{1,2,\dots,10000\\}$ on each margin. Every $\mathbf{x},\mathbf{y}$ pair is
recursively binned to every maximum depth between two and ten using random
binning, binning maximizing the chi score, and binning maximizing the mi score
at each step. Splitting is stopped when the expected number of observations
within a bin falls below ten or if a bin contains no observations. For each
binning, the $\chi^{2}$ statistic is computed for the final bins based on the
observed and expected counts, and this statistic is recorded along with the
number of final bins.

####  8.2.1 The impact of depth

To start, consider a plot of the $\chi^{2}$ statistic against the number of
bins for these simulated data when bin splits are selected randomly and
uniformly among observed point coordinates within a parent bin. Figure 7
displays this plot with both margins $\log_{10}$ transformed to make the
horizontal width of clusters for small depths as wide as those for large
depths. A dashed line is plotted on top of the points at the 99% critical
value of $\chi^{2}_{n_{bin}-1}$ for reference.

Figure 7: The $\chi^{2}$ statistic plotted by the final number of bins for
10,000 simulated independent data sets for each of 9 depth settings with
random binnings. The dashed line displays the 99% critical value of the
$\chi^{2}$ distribution with one degree of freedom less than the number of
bins, which is conservative for all depths.

Due to interference by the expected and minimum count stop criteria, the
maximum depth stop criterion does not produce the same number of bins every
time. If every bin were split at every step, a maximum depth of $d$ would
always produce $2^{d}$ bins but small bins are produced by chance which are
not further split even when the maximum depth is increased. Nonetheless, a
larger maximum depth will produce more bins as long as some can be split. This
produces the clusters of bin counts at each depth in Figure 7.

A more interesting result comes from a comparison between the 99% critical
value of the $\chi^{2}_{n_{bin}-1}$ distribution and statistic values for all
depths. Given that 10,000 null cases are generated for each depth, roughly 100
points are expected above this critical value if the 0.99 quantile of the
$\chi^{2}_{n_{bin}-1}$ corresponds with the 0.99 quantile of the null
distribution. Instead, far fewer points are observed above the line. This
suggests the $\chi^{2}_{n_{bin}-1}$ 99% critical value provides a conservative
approximation of this critical value for the null distribution under random
splitting. That is, choosing quantiles from this distribution to test for
dependence at level 0.01 implies a test with a level of at most 0.01 for the
null distribution. Indeed, quantile regression performed later suggests the
$\chi^{2}_{n_{bin}}$ distribution is generally a conservative approximation
for the null distribution of the $\chi^{2}$ statistic under random splitting.

The same result is not expected for splits chosen to maximize the chi or mi
score due to this deliberate maximization. The $\chi^{2}$ statistic for the
simulated data split to optimize the chi score is shown plotted against the
final number of bins in Figures 8 and the same plot for splits optimizing the
mi score are shown in 9. Points are, again, coloured according to the maximal
depth allowed, and the mean statistic values and number of bins for each of
these depth settings are denoted with corresponding coloured squares.
Additionally, the 99% critical value for the $\chi^{2}$ distribution with
$n_{bin}-1$ degrees of freedom is plotted against $n_{bin}$.

Figure 8: A comparison of the number of final bins and $\chi^{2}$ statistic
for 10,000 simulated independent data sets split to maximize the chi score
over each of 9 depth settings. The dashed line displays the 99% critical value
of the $\chi^{2}$ distribution with one degree of freedom less than the number
of bins. For all depths, many more points lie above this critical value than
would be expected if the data followed a $\chi^{2}_{n_{bin}-1}$ distribution.

Two changes occur in Figure 8 compared to random binning. First, the statistic
values tend to be larger than the null case for the same value of $n_{bin}$.
Indeed, the majority of score-maximizing binned statistics are above the
$\chi_{n_{bin}-1}^{2}$ 99% critical line when the maximal depth is only 4, all
observations are above the line when the depth is 6, and the gap between the
critical value curve and the true distribution grows larger as the number of
bins increases.

Second, the number of bins no longer displays clear separation by the
restriction on maximal depth. Whereas Figure 7 has clusters along the
horizontal margin corresponding to each maximal depth, the maximum binning bin
counts are not strongly grouped. This may be a result of the maximized method
selecting smaller bins on average than random binning, resulting in fewer bins
due to the minimum size and expected count restrictions that stop splitting of
these small bins early. This dynamic ‘smears’ the number of bins for each
maximum depth, leading to large overlap between groups by maximum depth.

Figure 9: A plot of the $\chi^{2}$ statistic by the final number of bins for a
binning which maximized the mi score at each step. The plot looks almost
identical to that generated by splitting the chi score.

Figure 9 shows a similar pattern for the mutual information statistic. The
$\chi^{2}$ statistic value rapidly increases in the number of bins, resulting
in a distribution which is well above the $\chi^{2}_{n_{bin}-1}$ 0.99
quantile. Indeed, both seem to produce similar null distributions for the
$\chi^{2}$ statistic when applied to the same data.

To better compare the $\chi^{2}_{n_{bin}-1}$ quantiles to the null
distributions across all numbers of bins, quantile regression of the statistic
value on the bin depth was carried out using simulated data and the quantreg
package in R (Koenker, 2023). Specifically, regression of the $0.95$, $0.99$,
and $0.999$ quantiles of the $\chi^{2}$ statistic was carried out for binning
under all three splitting methods. The results are shown in Figure 10. Just as
indicated by the earlier plots, maximized splitting with either score leads to
_inflated_ statistic values, that is statistic which is stochastically greater
than a $\chi^{2}_{n_{bin}-1}$ random variable. Both maximized methods produce
almost identical upper quantiles. In contrast, random splitting produces a
$\chi^{2}$ statistic which is stochastically less than a
$\chi^{2}_{n_{bin}-1}$ random variable, suggesting the $\chi^{2}_{n_{bin}}$
distribution could be used to generate conservative $p$-values (at least in
the tails).

Figure 10: Fit quantile regression lines for several quantiles for all
different splitting methods. Both maximized splitting methods lead to similar
larger quantiles that grow faster with depth and are further apart, though the
chi score leads to slightly larger $\chi^{2}$ statistic values.

These lines not only make the difference in the distributions under the
different splitting methods clear, but also could be used in practice to
determine rejection or acceptance of the null hypothesis of independent data.
Of course, the quantiles depend on the number of observations, so this plot is
only demonstrative unless the sample size is 10,000.

In short, both depth and splitting method are highly relevant to the
distribution of statistics computed on the final bins. The naive
$\chi_{n_{bin}-1}^{2}$ distribution on the bins seems to provide a
conservative distributional approximation for random splitting, if any
maximization is used it cannot be used for accurate rejection or to
approximate $p$-values. Indeed, the critical value of the statistic in this
case generated must take the number of bins into account. Most simply, a large
simulation as performed here can be used as a reference and the empirical
distributions could be used to generate critical values or empirical
$p$-values for each number of bins. These could be smoothed or interpolated to
cover gaps in the data due to certain bin counts occurring less frequently by
chance.

###  8.3 Depth and sample size

To investigate other sample sizes, 10,000 samples of both size $n=100$ and $n$
= 1,000 were generated and split to a maximum depth of 10. Random binning and
binning maximizing the chi score were applied to the collection of samples for
both sample sizes with the same stop criteria as the null investigation
immediately preceding. Plots of the resulting $\chi^{2}$ statistic and number
of bins are compared to the plot for $n$ = 10,000 in Figure 11 for the case of
random splitting and in Figure 12 for the case of spitting which maximizes the
chi score at each step. Splits maximizing the mi score were not investigated
here due to the very similar null distribution between the null under these
splits and when maximizing the chi score.

|  |
---|---|---
(a) | (b) | (c)
Figure 11: $\chi^{2}$ statistic and number of bins for 10,000 samples of simulated random data split by random binning to a maximum depth of ten with the minimum expected bin count restricted to 5 for (a) $n=100$, (b) $n$ = 1,000, and (c) $n$ = 10,000. As the sample size increases, the maximum number of bins and separation between the depths increases. |  |
---|---|---
(a) | (b) | (c)
Figure 12: $\chi^{2}$ statistic and number of bins for 10,000 samples of
random data split maximizing the chi score with a minimum expected bin count
of 5 to a maximum depth of ten for (a) $n=100$, (b) $n$ = 1,000, and (c) $n$ =
10,000. As before, the sample size limits the maximum number of bins produced
by the algorithm.

In both cases, the primary impact of increasing the sample size is to increase
the maximum number of bins produced by the algorithm. For smaller sample
sizes, the minimum bin size stop criteria and bin split restrictions end
splitting before many bins are created. This could have been anticipated. With
a sample size of 100, for example, splitting the data into 10 bins implies an
expected count of 10 per bin, at which point the stop criteria would cease
further splitting. Irregular bin shapes complicate this, of course, but the
basic idea stands. Fewer points means smaller expectations to compare to the
constant stop criteria.

Importantly, this does not seem to impact the conservative approximation of
the $\chi^{2}$ statistic under random binning by the $\chi^{2}_{n_{bin}-1}$
distribution. For all three sample sizes, the observed statistic values have
similar location and spread given the number of bins. As the sample size
increases, the statistics merely spread along the curve dictated by the
$\chi^{2}_{n_{bin}-1}$ distribution.

A more impactful pattern is seen for splitting maximizing the $\chi^{2}$
score. As the sample size increases, the whole distribution of the $\chi^{2}$
statistic moves upward for a given number of bins. This is particularly
obvious if the mean points are compared to the line providing the 99% quantile
of the $chi^{2}_{n_{bin}-1}$ distribution. In Figure 12(a), none of the mean
points lie above this line for $n_{bin}>10$ while in Figure 12(c) all of them
lie well above this line.

These different behaviours are probably a result of the way the chi score
maximizing algorithm chases local patterns, as demonstrated on the perfect
line data. Random binning splits agnostically, but chi maximizing binning
searches through all splits for the one leading to the largest $\chi^{2}$
statistic. Large $\chi^{2}$ values therefore occur by chance for random
binning but are actively sought by maximized binning. With a larger sample, a
greater variety of local patterns will be included in any realization which
can inflate the $\chi^{2}$ statistic. The maximized algorithm actively
searches for such patterns, stochastically increasing the resulting
distribution.

###  8.4 Simulated data patterns

Besides the null distribution, an important aspect to the use of recursive
binning to measure association is the ability of the alogrithm to detect non-
null patterns in the data. Newton (2009) provides supplemental code to
generate several interesting data patterns, which have since been used by Liu
et al. (2018) and Heller et al. (2016) to test their methods and also appear
on the Wikipedia article for Pearson’s correlation. This code is adapted and
used here to repeatedly generate observations following each pattern as a test
of recursive binning. Samples of 1,000 observations from each pattern are
shown in Figure 13, and the marginal ranks of these samples are shown in
Figure 14. Both figures have the margins dropped for brevity.

Figure 13: Point patterns used in Newton (2009) and Liu et al. (2018) to
demonstrate different dependence structures. In the last pattern, the data are
independent despite the obvious two-dimensional structure. Figure 14: The
patterns of Figure 13 converted to marginal ranks. Refer to these patterns as
the wave, rotated square, circle, valley, cross, ring, and uniform noise
respectively.

All but the final pattern contain non-linear dependence between the vertical
and horizontal variables, and still display strong patterns when these
variables are converted to ranks. The final sample, in contrast, consists of
four clusters of equal size which can be constructed by the product of
independent bimodal marginal densities. Consequently, in the ranks the sample
appears uniform. The final pattern therefore acts as a control to ensure that
structure without dependence is not detected by recursive binning applied to
the ranks.

For each pattern, both the chi score maximizing binning and random binning
algorithms were applied with maximum depths ranging from 1 to 10 and the same
stop criteria as previously (stop splitting for bins with an expected count of
ten or less and empty bins). This gives a sequence of $n_{bin}$, $\chi^{2}$
coordinates for each pattern as the depth restriction is increased, just as in
Figures 7 and 8. Rather than visualize these as points, the progression of
bins and statistics with increasing maximum depth can be emphasized by joining
the points for sucessive depth restrictions to give paths for each pattern in
the $n_{bin},\chi^{2}$ space. The same can be done for the 10,000 simulated
null samples of 1,000 points from Section 8.2 to give 10,000 null lines for
comparison. Figure 15 displays paths coloured by pattern to match Figure 14
alongside the null paths coloured in gray. A dashed line at the 0.95 quantile
of the $\chi^{2}_{n_{bin}}$ distribution is added to simulate a rejection
boundary in Figure 15(b).

|
---|---
(a) | (b)
Figure 15: Paths for every pattern in $n_{bin}$, $\chi^{2}$ statistic space
across depth restrictions for (a) maximized chi score splitting and (b) random
splitting. While the maximized splitting is deterministic for a sample,
leading to smooth curves, randomized binning leads to rough and erratic paths.
Despite this, both display roughly the same ordering of the patterns by
colour, and all patterns have paths far above any of the 10,000 simulated null
cases.

Figure 15(a) shows smooth curves increasing in both the number of bins and the
$\chi^{2}$ statistic as the maximum depth is increased while Figure 15(b)
shows erratic paths which occasionally decrease in the statistic value as
depth increases. Holding this major difference aside for a moment, all the
curves for patterns with dependence lie well above the null paths, while the
uniform pattern sits within the null curves for both splitting methods. This
suggests empirical $p$-values less than 0.0001 for every pattern when the
depth restriction is greater than 3. Additionally, the order of the curves is
similar between the two, with the wave and the cross giving the largest
$\chi^{2}$ statistic values for a given number of bins and the rotated square
and circle giving the smallest statistic values.

Returning to the difference in path smoothness, the erratic nature of the
random binning and smoothness of the maximized binning are natural
consequences of the different split methods. The maximized splitting algorithm
chooses the maximum split at each step and so behaves deterministically for a
given sample. In contrast, the random splitting algorithm proceeds non-
deterministically and will generate different bin counts and statistic values
for a given depth restriction every time it is run, even for constant data. To
better see the difference this makes for the separation of the patterns and to
evaluate the performance of both based on more than a single exemplar, samples
of 1,000 points from each pattern are generated independently and subsequently
split either by chi-maximizing or random binning. The resulting paths for each
of the 100 independent samples are displayed in Figure 16, with the median
path plotted using a thicker grey-accented line.

|
---|---
(a) | (b)
Figure 16: Paths of (a) maximized and (b) random binning applied to 100
independent realizations of the seven simulated data patterns compared to the
null paths with median paths plotted with thicker lines accented by grey. Both
spitting regimes create clear separation between the null paths and paths of
patterns with dependence and both display the same ordering of patterns. The
random splitting, however, creates more erratic and variable paths.

The clouds of paths show roughly the same ordering under both random binning
and chi-maximized binning. Both place the wave (in dark green) and the cross
(in light green) above the others, followed closely by the valley (in pink)
and the ring (in yellow). Halfway between these these patterns (all with
locally linear sections) and the null patterns sit the rotated square (in
orange) and circle (in blue). Though these latter two result in smaller
statistic values than the others, they remain easily distinguished from the
null paths in both Figures 16(a) and 16(b). Indeed, under random splitting the
$\chi^{2}_{n_{bin}-1}$ 0.95 quantile plotted with a dashed line neatly
separates the null paths from these two patterns and consequently all
patterns.

The close agreement of their ordering suggests that randomized binning is as
powerful as maximized binning at detecting these patterns. The clear
separation of the null curves from the curves of all patterns over all
repetitions implies empirical $p$-values $<0.0001$ for all patterns when the
algorithm is allowed to reach its maximal depth using only bin size stop
criteria. Any of these patterns would therefore be regarded as significant by
both regimes for commonly used rejection levels. Moreover, as the order of
curves by pattern is the same for both, the greatest difference between the
maximized binning and random binning seems to be the increased the statistic
value alone. If a rejection level is chosen on the $\chi^{2}$ statistic
proportional to the null quantile, we expect that both regimes would show the
same tendency for acceptance or rejection of every pattern. Given that random
binning requires less computation and has a conservative $p$-value provided by
the $\chi^{2}_{n_{bin}-1}$ distribution, the similar ordering of these curves
provides a strong argument for the use of random recursive binning of the
ranks to detect and quantify dependence between variables.

####  8.4.1 Progression of bins in the maximizing algorithm

The median lines are somewhat more interesting, as they cross each other
frequently. This suggests the ordering of these different patterns, and
potentially the power of the algorithm under each splitting method to detect
them, is dependent on the maximum depth. Under chi-maximized binning, the
median line for the rotated square is above all others for very small numbers
of bins, indicating there is a pattern in the data which can be captured
quickly with few splits. In contrast, the median line for the ring pattern
starts below most others for both splitting methods, indicating that a certain
number of bins are required before the pattern is detected. It seems each
pattern has a natural _resolution_ : a number of bins required to identify the
dependence.

Directly related are the different slopes of these curves for each pattern.
Patterns which generate large residuals within a few splits will grow rapidly
at first, and then may slow if the following residuals are relatively small.
The median path of the valley pattern, for example, grows quicker than most
other patterns in the number of bins for the first few depths before its rate
of growth slows. In contrast, the $\chi^{2}$ statistic of the ring pattern is
unexceptional before a a certain depth, at which point the statistic values
jump abruptly. More insight into both the slope of these lines and the
resolution of each pattern can be gleaned by plotting the bins at each depth
for each pattern shaded by residual as in Figure 17.

---
Figure 17: Bins for chi score maximizing splits at increasing depths. By row,
depths 2, 4, 6, 8, and 10 are displayed in order. The final bins reflect and
summarize the pattern of points given to the algorithm.

As before, the hue of the shading is determined by the sign of the bin’s
Pearson residual, blue for negative and red for positive, and the saturation
is determined by its magnitude. To ensure fair comparison of the residuals,
all saturations are determined relative to the maximum residual observed
across all depths. The rapid early growth of the rotated square can be spied
immediately by its shading in the first row. After only two splits, the empty
top left corner and dense top right are detected by the maximized binning
algorithm and contribute large residuals, leading to a large early $\chi^{2}$
statistic. Once these regions of low and high density around the margins are
identified everywhere, however, the algorithm is left splitting the nearly
uniform interior of the rotated square, and so grows parallel to the null
paths. One interpretation of this “elbow” point is the point at which the
recursive binning has adequately captured the pattern in the data. Additional
splits beyond this point are ineffective at increasing the score because they
occur essentially randomly based on previous splits.131313Another
interpretation of this behaviour is that the algorithm is effectively
splitting the noise. Once large areas of relatively high and low density have
been identified, the main aspect of the data determining whether further
splits are productive is the level of noise. So, for example, the perfect line
of Figure 6(b) always benefits from further splits, but the noisier pattern of
the valley is adequately captured at some point.

In contrast, the pattern observed for the ring requires a certain depth to
detect. The radial symmetry of this pattern means that most early splits fail
to identify the less dense regions in the corners and the centre and the more
dense region along the circumference of the ring. It is only as the depth
exceeds 6 that strong positive and negative residuals are found in any bins,
before that point the residuals are not as large as for the other patterns.

Comparing the final row of Figure 17 to randomized binning at a depth of ten
in Figure 18, a clear advantage of the maximized binning is better
representation of the the underlying pattern. Especially for the cross, wave,
and valley, random binning produces a larger proportion of thin bins which
obscure the pattern of residuals and points. In contrast, maximized binning
tends to chase local patterns, leading to many small rectangular bins around
high density areas that provide a sense of the pattern they summarize. Should
a visual summary of the data be desired, maximized binning will give more
consistent and clear results than random binning.

Figure 18: Bins for random splitting at a depth of 10

##  9 A real data example: S&P 500 returns

For a real data example, consider the S&P 500 constituent data from Hofert and
Oldford (2018). The raw data contains a time series of 505 stock prices from
the first day of 2007 to the last day of 2009 for stocks included in the S&P
500 index. The goal is to evaluate the pairwise dependence present within the
negative log-returns of the 461 stocks with complete records over this period
using recursive binning over all ${461\choose 2}=$ 106,030 pairs. The negative
log-return of a stock is the negative logarithm of the ratio in its end-of-day
price over two consecutive days, explicitly

$-\log\frac{S_{t}}{S_{t-1}}$

for a stock with value $S_{t}$ at time $t$. As there are 756 days recorded in
the data set, there are 755 log-returns for each stock.

To remove further time dependencies between stock returns day-to-day, the
negative log-returns for each stock are fit with an ARMA(1,1)-GARCH(1,1) model
and the residuals are taken as a new set of independent pseudo-observations,
see details in Hofert and Oldford (2018). The raw data are taken from the
qrmdata package (Hofert et al., 2022) and processed by code adapted from the
SP500 demo from the zenplots package (Hofert and Oldford, 2020) in order to
compute the log-returns and generate the pseudo-observations for recursive
binning.

Recursive binning as in Section 8.4 is applied to each of the 106,030 pairs of
755 pseudo-observations to evaluate dependence. Specifically, splitting
maximizes the chi score under the constraint that no bins be created with
expected counts less than 5 and splitting is stopped when the expected count
in a bin is less than 10, the bin contains no observations, or the depth of
the bin is 6. After binning, the $\chi^{2}$ statistic is computed over the
bins to measure the departure of the observed distribution from that expected
under uniformity. Figure 19 displays the resulting statistic values and
$n_{bin}$ for every pair compared to the null quantiles for uniform pairs with
1,000 points from Section 8.2 estimated by quantile regression141414This
choice of null should be somewhat conservative for the case of 755 points. As
shown in Figure 12, smaller sample sizes lead to fewer bins and smaller
$\chi^{2}$ statistic values for the same maximum depth. Therefore, taking the
$\chi^{2}$ statistic over simulated uniform samples of a slightly larger
sample size provides an approximately correct, but slightly conservative, null
distribution.. Hues from blue to red encode the empirical $p$-value of an
observed statistic for an S&P 500 pair at a given number of bins in the null
data. Partial transparency (alpha-blending) is also used in this plot, but the
sheer number of points makes it ineffective at representing the density, so
marginal histograms have been added to give a better sense of the location
where values are most concentrated.

Figure 19: The distribution of $n_{bin}$ and $\chi^{2}$ statistics for the
106,030 S&P 500 pairs split to a maximum depth of 6 compared to some upper
quantile estimates of the $\chi^{2}$ statistic from $n_{bin}$ under the null
distribution of independence. Most pairs in the S&P 500 data appear to be
highly significant.

Most striking in this plot is the significance of nearly every pair. Only 334
of the 106,030 pairs (0.3%) have empirical $p$-values less than 0.95, and the
$\chi^{2}$ statistic values are centred well above the fit null quantiles
despite their slightly conservative nature. The pairs generally lie in a
single large cluster in $\chi^{2}$ statistic and $n_{bin}$ values and the
marginal histograms indicate that within this cluster most points are
concentrated in a small region at its center. Only a few dozen pairs lie
outside this cluster, including the two very large $\chi^{2}$ statistic values
for relatively small values of $n_{bin}$ and the smallest $\chi^{2}$ statistic
values with the smallest values of $n_{bin}$. This suggests that almost every
pair in this data set contains some level of dependence, inviting further
exploration.

Figure 20: The S&P 500 pairs with the largest $\chi^{2}$ statistic values
after recursive binning under the splitting and stopping logic described
earlier. All pairs show strong linear dependence, especially in the upper and
lower tails.

First, consider the exceptionally large statistic values. Figure 20 displays a
matrix of the 36 pairs with the largest $\chi^{2}$ statistic values in
decreasing statistic order from top left to bottom right. For each pair, a
scatterplots of the marginal ranks is augmented by a plot of the final
binnings coloured by the sign of the Pearson residual (red for positive and
blue for negative) and shaded by magnitude. The range of hues is kept constant
through all subplots to support direct comparisons between any two binnings.

Immediately apparent in every plot is a strong positive linear relationship
between the pairs. Bins along the diagonal line from the bottom left to the
top right have more points than expected, while those in the top left and
bottom right corners have fewer. Most notably, the bins in the top right and
bottom left corners tend to have far more observations than expected. This
suggests particularly strong _tail dependence_ in these pairs, loosely the
probability that large or small values of two random variables occur
simultaneously. Precisely, upper tail dependence at $p\in[0,1]$ is defined as

$P\big{(}Z_{1}>Q_{1}(p)~{}\vline~{}Z_{2}>Q_{2}(p)\big{)}$

and lower tail dependence as

$P\big{(}Z_{1}\leq Q_{1}(p)~{}\vline~{}Z_{2}\leq Q_{2}(p)\big{)}$

for random variables $Z_{1},Z_{2}$ with respective quantile functions
$Q_{1}(\cdot),Q_{2}(\cdot)$ (Hofert and Oldford, 2018). As the upper and lower
bins approximate these conditional probabilities, albeit for separate quantile
values, the shaded residuals in the tails of these plots communicate how much
larger the conditional tail probabilities are than would be expected under
independence.

Though detecting this tail dependence is not the goal of recursive binning,
the most interesting pairs it identifies overlap considerably with the pairs
that have the largest upper tail dependence. Indeed, 8 the top 10 pairs ranked
by upper tail dependence in Hofert and Oldford (2018) can be found in Figure
20. The top two relationships, in particular, correspond to different classes
of stock in the same company. This explains not only their strong dependence
(especially in the tails), but also their outlying positions in Figure 19 with
similarly large $\chi^{2}$ statistics for relatively few bins.

At the other end of the spectrum, Figure 21 displays the pairs with the
smallest $\chi^{2}$ statistic values in decreasing statistic order from top
left to bottom right. Residuals are shaded according to the same hues as
Figure 20, so a comparison of bins between the two plots is possible. In
contrast to the pairs with strong association, the points, bins, and residuals
show no obvious patterns for any of these pairs. The magnitude of residuals is
generally small and the points show essentially random scatter. It is
unsurprising that these pairs are deemed unexceptional by the recursive
binning algorithm.

Figure 21: The S&P 500 pairs with the smallest $\chi^{2}$ statistic values.
These pairs look like realizations of random uniform data, with no strong
patterns to the bins, points, or residuals.

Finally consider the pairs with moderate statistic values. Figure 22 displays
the 36 pairs with values closest to the median $\chi^{2}$ statistic across all
106,030 pairs. Again, the shading is consistent with all previous plots to
allow comparisons. These ‘middling’ pairs show relatively weak positive,
linear relationships without the strong tail dependence of the pairs with the
largest $\chi^{2}$ statistic values. Nonetheless, there is a concentration of
bins which are shaded faintly red about the main diagonal and others shaded
faintly blue in the top left and bottom right corners distinctive of a
positive linear relationship. This observation, along with clear separation of
statistic values for the simulated data patterns of Section 8.4 from the null
distribution, supports the conclusion that the large statistic values across
the majority of pairs in the S&P 500 data are not spurious, but instead
reflect the power of recursive binning to detect even weak dependence.

Figure 22: S&P 500 pairs with $\chi^{2}$ statistic values nearest to the
median value. Though the dependence is weak compared to Figure 20, there is
still a concentration of bins with more points than expected along the
diagonal and bins with fewer points than expected in the off-diagonal corners,
suggesting weak positive linear relationships between these pairs.

##  10 Conclusions

This investigation of recursive binning as a method to measure association
garners several key observations. First, and perhaps most important, the power
of recursive binning to detect association does not appear to depend greatly
on whether chi maximizing or randomly placed edges are used to split bins.
Both the maximizing splits and the random splits showed pronounced separation
of non-null patterns from the null distribution in the simulation study of
Section 8.4. Using random splits gives a number of advantages: it is
computationally faster, conceptually simpler, and produces a null distribution
which is conservatively approximated by the $\chi^{2}_{n_{bin}-1}$
distribution if the spliting logic and stop criteria maintain the rule of
thumb that the expected count of every bin be $\geq 5$. If detecting
dependence is the only goal of an investigation, random splits may therefore
be preferred to maximizing ones despite, or perhaps because of, their
simplicity.

In contrast, maximized score binning produces bins which serve as a superior
visual summary of a pairwise relationship. Splits which maximize the chi
score, for example, will split off empty sections of a bin and separate
regions of particularly high density. The end result is a collection of bins
which represent the underlying pattern of points much better than the bins
which result from random splits. This comes at a cost of inflated statistic
values, requiring modeling or simulation at every application to estimate the
significance of an observed statistic value.

In either case, the pattern present in the data impacts the path of the
algorithm. For simple linear patterns, dependence can be detected with
relatively few splits and bins. More complex patterns may take many splits to
detect. For many patterns tested here, it seems a natural depth or resolution
is present. Splits below this natural depth may increase the $\chi^{2}$
statistic drastically, but after it is reached the statistic grows more slowly
at a rate comparable to the null over successive splits.

Recursive binning is a promising method to measure association. It displays
high power in the detection of non-linear relationships, can be used to
generate a summary visualization of data pairs, and seems to naturally
highlight local dependence such as tail dependence in real data. Though the
simulations here are not comprehensive, they are highly suggestive of a
practical and powerful tool for sorting and summarizing pairwise relationships
in large data sets.

## References

  * Agresti (1981) Alan Agresti.  Measures of nominal-ordinal association.  _Journal of the American Statistical Association_ , 76(375):524–529, 1981.
  * Anscombe (1973) Francis J. Anscombe.  Graphs in statistical analysis.  _The American Statistician_ , 27(1):17–21, 1973\.
  * Breiman (2001) Leo Breiman.  Statistical modeling: The two cultures.  _Statistical Science_ , 16(3):199–231, 2001.
  * Breiman and Friedman (1985) Leo Breiman and Jerome H. Friedman.  Estimating optimal transformations for multiple regression and correlation.  _Journal of the American Statistical Association_ , 80(391):580–598, 1985.
  * Cao et al. (2021) Dan Cao, Yuan Chen, Jin Chen, Hongyan Zhang, and Zheming Yuan.  An improved algorithm for the maximal information coefficient and its application.  _Royal Society Open Science_ , 8(2):201424, 2021\.
  * Carr et al. (1987) Daniel B. Carr, Richard J. Littlefield, W.L. Nicholson, and J.S. Littlefield.  Scatterplot matrix techniques for large N.  _Journal of the American Statistical Association_ , 82(398):424–436, 1987.
  * Chen et al. (2016) Yuan Chen, Ying Zeng, Feng Luo, and Zheming Yuan.  A new algorithm to optimize maximal information coefficient.  _PLoS One_ , 11(6):e0157567, 2016.
  * Choi et al. (2010) Seung-Seok Choi, Sung-Hyuk Cha, and Charles C. Tappert.  A survey of binary similarity and distance measures.  _Journal of Systemics, Cybernetics and Informatics_ , 8(1):43–48, 2010.
  * Cochran (1952) William G. Cochran.  The $\chi^{2}$ test of goodness of fit.  _The Annals of Mathematical Statistics_ , pages 315–345, 1952.
  * Cohen (1960) Jacob Cohen.  A coefficient of agreement for nominal scales.  _Educational and Psychological Measurement_ , 20(1):37–46, 1960.
  * Cramér (1924) Harald Cramér.  Remarks on correlation.  _Scandinavian Actuarial Journal_ , 1924(1):220–240, 1924.
  * Ding et al. (2017) A. Adam Ding, Jennifer G. Dy, Yi Li, and Yale Chang.  A robust-equitable measure for feature ranking and selection.  _The Journal of Machine Learning Research_ , 18(1):2394–2439, 2017.
  * Dümcke et al. (2014) Sebastian Dümcke, Ulrich Mansmann, and Achim Tresch.  A novel test for independence derived from an exact distribution of $i$th nearest neighbours.  _PloS One_ , 9(10), 2014.
  * Embrechts et al. (2001) Paul Embrechts, Filip Lindskog, and Alexander McNeil.  Modelling dependence with copulas.  Technical report, ETH Zurich, 2001.
  * Fairfield Smith (1957) H. Fairfield Smith.  On comparing contingency tables.  _The Philippine Statistician_ , 6:71–81, 1957.
  * Friedman and Stuetzle (2002) Jerome H. Friedman and Werner Stuetzle.  John W. Tukey’s work on interactive graphics.  _Annals of Statistics_ , pages 1629–1639, 2002.
  * Garcia et al. (2012) Salvador Garcia, Julian Luengo, José Antonio Sáez, Victoria Lopez, and Francisco Herrera.  A survey of discretization techniques: Taxonomy and empirical analysis in supervised learning.  _IEEE transactions on Knowledge and Data Engineering_ , 25(4):734–750, 2012.
  * Genest and Rémillard (2004) Christian Genest and Bruno Rémillard.  Test of independence and randomness based on the empirical copula process.  _Test_ , 13(2):335–369, 2004.
  * Goodman and Kruskal (1979) Leo A. Goodman and William H. Kruskal.  _Measures of Association for Cross Classifications_.  Springer, 1979.
  * Gorfine et al. (2012) Malka Gorfine, Ruth Heller, and Yair Heller.  Comment on ”Detecting Novel Associations In Large Data Sets”.  _Science_ , pages 1–6, 2012.
  * Gretton et al. (2007) Arthur Gretton, Kenji Fukumizu, Choon Teo, Le Song, Bernhard Schölkopf, and Alex Smola.  A kernel statistical test of independence.  _Advances in Neural Information Processing Systems_ , 20, 2007.
  * Hald (1998) Anders Hald.  _A History of Mathematical Statistics from 1750 to 1930_.  Wiley, 1998.
  * Heller et al. (2013) Ruth Heller, Yair Heller, and Malka Gorfine.  A consistent multivariate test of association based on ranks of distances.  _Biometrika_ , 100(2):503–510, 2013.
  * Heller et al. (2016) Ruth Heller, Yair Heller, Shachar Kaufman, Barak Brill, and Malka Gorfine.  Consistent distribution-free k-sample and independence tests for univariate random variables.  _The Journal of Machine Learning Research_ , 17(1):978–1031, 2016.
  * Hoeffding (1948) Wassily Hoeffding.  A non-parametric test of independence.  _The Annals of Mathematical Statistics_ , pages 546–557, 1948.
  * Hofert and Oldford (2018) Marius Hofert and R. Wayne Oldford.  Visualizing dependence in high-dimensional data: An application to S&P 500 constituent data.  _Econometrics and Statistics_ , 8:161–183, 2018.
  * Hofert and Oldford (2020) Marius Hofert and R. Wayne Oldford.  Zigzag expanded navigation plots in R: The R package zenplots.  _Journal of Statistical Software_ , 95(1):1–44, 2020.
  * Hofert et al. (2022) Marius Hofert, Kurt Hornik, and Alexander J. McNeil.  _qrmdata: Data Sets for Quantitative Risk Management Practice_ , 2022\.  URL https://CRAN.R-project.org/package=qrmdata.  R package version 2022-05-31-1.
  * Jiang and Liu (2013) Bo Jiang and Jun S. Liu.  Sliced inverse regression with variable selection and interaction detection.  _arXiv preprint arXiv:1304.4056_ , 652:17, 2013.
  * Jiang et al. (2015) Bo Jiang, Chao Ye, and Jun S. Liu.  Nonparametric k-sample tests via dynamic slicing.  _Journal of the American Statistical Association_ , 110(510):642–653, 2015.
  * Khamis (2008) Harry Khamis.  Measures of association: How to choose?  _Journal of Diagnostic Medical Sonography_ , 24(3):155–162, 2008.
  * Kinney and Atwal (2014a) Justin B. Kinney and Gurinder S. Atwal.  Equitability, mutual information, and the maximal information coefficient.  _Proceedings of the National Academy of Sciences_ , 111(9):3354–3359, 2014a.
  * Kinney and Atwal (2014b) Justin B. Kinney and Gurinder S. Atwal.  Reply to Reshef et al.: Falsifiability or bust.  _Proceedings of the National Academy of Sciences_ , 111(33):E3364–E3364, 2014b.
  * Koenker (2023) Roger Koenker.  _quantreg: Quantile Regression_ , 2023.  URL https://cran.r-project.org/web/package=quantreg.  R package version 5.97.
  * Lee and Huh (2003) Seung-Chun Lee and Moon Yul Huh.  A measure of association for complex data.  _Computational statistics & data analysis_, 44(1-2):211–222, 2003.
  * Liebetrau (1983) Albert M. Liebetrau.  _Measures of association_.  Sage, 1983.
  * Liu et al. (2018) Ying Liu, Victor de la Pena, and Tian Zheng.  Kernel-based measures of association.  _WIREs Computational Statistics_ , 10(2), 2018.
  * Lopez-Paz et al. (2013) David Lopez-Paz, Philipp Hennig, and Bernhard Schölkopf.  The randomized dependence coefficient.  _arXiv preprint arXiv:1304.7717_ , 2013.
  * Móri and Székely (2019) Tamás F. Móri and Gábor J. Székely.  Four simple axioms of dependence measures.  _Metrika_ , 82(1):1–16, 2019.
  * Newton (2009) Michael A Newton.  Introducing the discussion paper by Székely and Rizzo.  _The Annals of Applied Statistics_ , 3(4):1233–1235, 2009.
  * Pearson (1900) Karl Pearson.  On the criterion that a given system of deviations from the probable in the case of a correlated system of variables is such that it can be reasonably supposed to have arisen from random sampling.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 50(302):157–175, 1900.
  * Plackett (1983) Robin L. Plackett.  Karl Pearson and the chi-squared test.  _International Statistical Review_ , pages 59–72, 1983.
  * Rahman (2018) Adam Rahman.  _Preserving measures structure during generation and reduction of multivariate point configurations_.  PhD thesis, University of Waterloo, 2018.  URL http://hdl.handle.net/10012/13365.
  * Reimherr and Nicolae (2013) Matthew Reimherr and Dan L. Nicolae.  On quantifying dependence: a framework for developing interpretable measures.  _Statistical Science_ , 28(1):116–130, 2013.
  * Rényi (1959) Alfréd Rényi.  On measures of dependence.  _Acta Mathematica Academiae Scientiarum Hungarica_ , 10(3-4):441–451, 1959.
  * Reshef et al. (2011) David N. Reshef, Yakir A. Reshef, Hilary K. Finucane, Sharon R. Grossman, Gilean McVean, Peter J. Turnbaugh, Eric S. Lander, Michael Mitzenmacher, and Pardis C. Sabeti.  Detecting novel associations in large data sets.  _Science_ , 334(6062):1518–1524, 2011.
  * Reshef et al. (2014) David N. Reshef, Yakir A. Reshef, Michael Mitzenmacher, and Pardis C. Sabeti.  Cleaning up the record on the maximal information coefficient and equitability.  _Proceedings of the National Academy of Sciences_ , 111(33):E3362–E3363, 2014.
  * Reshef et al. (2018) David N. Reshef, Yakir A. Reshef, Pardis C. Sabeti, and Michael Mitzenmacher.  An empirical study of the maximal and total information coefficients and leading measures of dependence.  _The Annals of Applied Statistics_ , 12(1):123–155, 2018.
  * Schweizer and Wolff (1981) Berthold Schweizer and Edward F. Wolff.  On nonparametric measures of dependence for random variables.  _Annals of Statistics_ , 9(4):879–885, 1981.
  * Scott (1988) David W. Scott.  A note on choice of bivariate histogram bin shape.  _Journal of Official Statistics_ , 4(1):47, 1988\.
  * Shannon (1948) Claude E. Shannon.  A mathematical theory of communication.  _The Bell System Technical Journal_ , 27(3):379–423, 1948.
  * Siburg and Stoimenov (2010) Karl Friedrich Siburg and Pavel A. Stoimenov.  A measure of mutual complete dependence.  _Metrika_ , 71(2):239–251, 2010.
  * Silvey (1964) S.D. Silvey.  On a measure of association.  _The Annals of Mathematical Statistics_ , pages 1157–1166, 1964.
  * Simon and Tibshirani (2014) Noah Simon and Robert Tibshirani.  Comment on ”Detecting Novel Associations In Large Data Sets” by Reshef et al., Science dec 16, 2011.  _arXiv preprint arXiv:1401.7645_ , 2014.
  * Sklar (1996) Abe Sklar.  Random variables, distribution functions, and copulas: a personal look backward and forward.  _Institute of Mathematical Statistics Lecture Notes – Monograph Series_ , pages 1–14, 1996.
  * Stephens (1974) Michael A. Stephens.  Edf statistics for goodness of fit and some comparisons.  _Journal of the American Statistical Association_ , 69(347):730–737, 1974.
  * Stigler (1989) Stephen M Stigler.  Francis Galton’s account of the invention of correlation.  _Statistical Science_ , pages 73–79, 1989.
  * Székely and Rizzo (2009) Gábor J Székely and Maria L. Rizzo.  Brownian distance covariance.  _The Annals of Applied Statistics_ , pages 1236–1265, 2009.
  * Thas and Ottoy (2004) Olivier Thas and Jean-Pierre Ottoy.  A nonparametric test for independence based on sample space partitions.  _Communications in Statistics-Simulation and Computation_ , 33(3):711–728, 2004.
  * Tjøstheim et al. (2022) Dag Tjøstheim, Håkon Otneim, and Bård Støve.  Statistical dependence: Beyond pearson’s $\rho$.  _Statistical science_ , 37(1):90–109, 2022.
  * Wilkinson and Wills (2008) Leland Wilkinson and Graham Wills.  Scagnostics distributions.  _Journal of Computational and Graphical Statistics_ , 17(2):473–491, 2008.
  * Wilkinson et al. (2005) Leland Wilkinson, Anushka Anand, and Robert Grossman.  Graph-theoretic scagnostics.  In _IEEE Symposium on Information Visualization, 2005. INFOVIS 2005._ , pages 157–164. IEEE, 2005.
  * Zheng et al. (2021) Wenjun Zheng, Dejian Lai, and K. Lance Gould.  A simulation study of a class of nonparametric test statistics: a close look of empirical distribution function-based tests.  _Communications in Statistics-Simulation and Computation_ , pages 1–17, 2021.