Bayesian Dependence Tests for Continuous, Binary and Mixed Continuous-Binary Variables

Abstract

Tests for dependence of continuous, discrete and mixed continuous-discrete variables are ubiquitous in science. The goal of this paper is to derive Bayesian alternatives to frequentist null hypothesis significance tests for dependence. In particular, we will present three Bayesian tests for dependence of binary, continuous and mixed variables. These tests are nonparametric and based on the Dirichlet Process, which allows us to use the same prior model for all of them. Therefore, the tests are “consistent” among each other, in the sense that the probabilities that variables are dependent computed with these tests are commensurable across the different types of variables being tested. By means of simulations with artificial data, we show the effectiveness of the new tests.

1. Introduction

Tests for dependence of continuous, discrete and mixed continuous-discrete variables are fundamental in science. The standard way to statistically assess if two (or more) variables are dependent is by using null-hypothesis significance tests (NHST), such as ${\chi }^2$-test, Kendall’s τ, etc. However, these tests are affected by the drawbacks which characterize NHST [1,2,3]. An NHST computes the probability of getting the observed (or a larger) value of the statistics under the assumption that the null hypothesis of independence is true, which is obviously not the same as the probability of variables being dependent on each other, given the observed data. Another common problem is that the claimed statistical significance might have no practical impact. Indeed, the usage of NHST often relies on the wrong assumptions that p-values are a reasonable proxy to the probability of the null hypothesis and that statistical significance implies practical significance.

In this paper, we propose a collection of Bayesian dependence tests. The questions we are actually interested in—for example, Is variable Y dependent on Z? or Based on the experiments, how probable is Y dependent on Z?—are actually questions about posterior probabilities. Answers to these questions are naturally provided by Bayesian methods. The core of this paper is thus to derive Bayesian alternatives to frequentist NHST and to discuss their inference and results. In particular, we present three Bayesian tests for dependence of binary, continuous and mixed variables. All of these tests are nonparametric and based on the Dirichlet Process. This allows us to use the same prior model for all the tests we develop. Therefore, they are “consistent” in the sense that the probabilities of dependence we compute are commensurable across the tests. This is another main difference about such an approach and the use of p-values, since the latter usually cannot be compared across different types of tests.

To address the issue of how to choose the prior parameters in case of lack of information, we propose the use of the Imprecise Dirichlet Process (IDP) [4]. It consists of a family of Dirichlet processes with fixed prior strength and and prior probability measure free to span the set of all distributions. In this way, we obtain as a byproduct a measure of sensitivity of inferences to the choice of the prior parameters.

Nonparametric tests based on the Dirichlet Process and on similar ideas to those presented in this paper have also been proposed in [4] to develop a Bayesian rank test, in [5] for a Bayesian signed-rank test, in [6] for a Bayesian Friedman test and in [7] for a Bayesian test that accounts for censored data.

Several alternative Bayesian methods are available for testing of independence. The test of linear dependence between two continuous univariate random variables can be achieved by fitting a linear model and inspecting the posterior distribution of the correlation coefficient. A more sophisticated test based on a Dirichlet Process Mixture prior is instead presented in [8] to deal with linear and nonlinear dependences. Other methods were proposed for testing of independence based on a contingency table [9,10,11]. The main difference between these works and the work presented in this paper is that we provide tests for continuous, categorical (binary) and mixed variables using the same approach. This allows us to derive a very general framework to test independence/dependence (these tests could be used for instance for feature selection in machine learning [12,13,14,15]).

By means of simulations on artificial data, we use our test to decide if two variables are dependent. We show that our Bayesian test achieves equal or better results than the frequentist tests. We moreover show that the IDP test is more robust, in the sense that it acknowledges when the decision is prior-dependent. In other words, the IDP test suspends the judgment and becomes indeterminate when the decision becomes prior dependent. Since IDP has all the positive features of a Bayesian test and it is more reliable than the frequentist tests, we propose IDP as a new test for testing dependence.

2. Dirichlet Process

The Dirichlet Process was developed by Ferguson [16] as a probability distribution on the space of probability distributions. Let $\mathbb{X}$ be a standard Borel space with Borel σ-field ${\mathcal{B}}_{\mathbb{X}}$ and $\mathbb{P}$ be the space of probability measures on $(\mathbb{X},{\mathcal{B}}_{\mathbb{X}})$ equipped with the weak topology and the corresponding Borel σ-field ${\mathcal{B}}_{\mathbb{P}}$. Let $\mathbb{M}$ be the class of all probability measures on $(\mathbb{P},{\mathcal{B}}_{\mathbb{P}})$. We call the elements $\mu \in \mathbb{M}$ nonparametric priors.

An element of $\mathbb{M}$ is called a Dirichlet Process distribution $\mathcal{D}(\alpha )$ with base measure α if for every finite measurable partition $B_1,\cdots ,B_m$ of $\mathbb{X}$, the vector $(P(B_1),\cdots ,P(B_m))$ has a Dirichlet distribution with parameters $(\alpha (B_1),\cdots ,\alpha (B_m))$, where $\alpha (·)$ is a finite positive Borel measure on $\mathbb{X}$. Consider the partition $B_1=A$ and $B_2=A^c=\mathbb{X}\setminus A$ for some measurable set $A\in \mathbb{X}$, then if $P\sim \mathcal{D}(\alpha )$ from the definition of the DP we have that $(P(A),P(A^c))\sim Dir(\alpha (A),\alpha (\mathbb{X})-\alpha (A))$, which is a β distribution. From the moments of the β distribution, we can thus derive that:

$$\begin{array}{cc}\hfill \mathcal{E}[P(A)]={\displaystyle \frac{\alpha (A)}{\alpha (\mathbb{X})}},& \mathcal{E}[{(P(A)-\mathcal{E}[P(A)])}^2]={\displaystyle \frac{\alpha (A)(\alpha (\mathbb{X})-\alpha (A))}{(\alpha {(\mathbb{X})}^2(\alpha (\mathbb{X})+1))}},\hfill \end{array}$$

(1)

where we have used the calligraphic letter $\mathcal{E}$ to denote expectation with respect to the Dirichlet process. This shows that the normalized measure $\alpha (·)/\alpha (\mathbb{X})$ of the DP reflects the prior expectation of P, while the scaling parameter $\alpha (\mathbb{X})$ controls how much P is allowed to deviate from its mean $\alpha (·)/\alpha (\mathbb{X})$. Let $s=\alpha (\mathbb{X})$ stand for the total mass of $\alpha (·)$ and ${\alpha }^{*}(·)=\alpha (·)/s$ stand for the probability measure obtained by normalizing $\alpha (·)$. If $P\sim \mathcal{D}(\alpha )$, we shall also describe this by saying $P\sim Dp(s,{\alpha }^{*})$ or, if $\mathbb{X}=\mathbb{R}$, $P\sim Dp(s,G_0)$, where $G_0$ stands for the cumulative distribution function of ${\alpha }^{*}$.

Let $P\sim Dp(s,{\alpha }^{*})$ and f be a real-valued bounded function defined on $(\mathbb{X},\mathcal{B})$. Then the expectation with respect to the Dirichlet Process of $E[f]$ is

$$\mathcal{E}\left[E(f)\right]=\mathcal{E}\left[\int fdP\right]=\int fd\mathcal{E}[P]=\int fd{\alpha }^{*}.$$

(2)

One of the most remarkable properties of the DP priors is that the posterior distribution of P is again a DP. Let $X_1,\cdots ,X_n$ be an independent and identically distributed sample from P and $P\sim Dp(s,{\alpha }^{*})$, then the posterior distribution of P given the observations, denoted as $P_{X|X^n}$, is

$$P_{X|X^n}\sim Dp\left(s+n,\frac{s}{s+n}{\alpha }^{*}+\frac{1}{s+n}\sum _{i=1}^n{\delta }_{X_i}\right),$$

(3)

where ${\delta }_{X_i}$ is an atomic probability measure centered at $X_i$ and $X^n=\{X_1,\cdots ,X_n\}$. This means that the Dirichlet Process satisfies a property of conjugacy, in the sense that the posterior for P is again a Dirichlet Process with updated unnormalized base measure $\alpha +{\sum }_{i=1}^n{\delta }_{X_i}$. From Equations (1)–(3), we can easily derive the posterior mean and variance of $P(A)$ and, respectively, posterior expectation of f. Hereafter we list some useful properties of the DP that will be used in the sequel (see Chapter 3 in [17]).

(a)

In case $\mathbb{X}=\mathbb{R}$, since P is completely defined by its cumulative distribution function F, a-priori we say $F\sim Dp(s,G_0)$ and a posteriori we can rewrite (3) as follows:

$$F_{X|X^n}\sim Dp\left(s+n,\frac{s}{s+n}{G}_0+\frac{n}{s+n}\frac{1}{n}\sum _{i=1}^n{I}_{[X_i,\infty )}\right),$$

(4)

where I is the indicator function and $\frac{1}{n}{\sum }_{i=1}^n{I}_{[X_i,\infty )}$ is the empirical cumulative distribution.

(b)

Consider an element $\mu \in \mathbb{M}$ which puts all its mass at the probability measure $P={\delta }_x$ for some $x\in \mathbb{X}$. This can also be modeled as $Dp(s,{\delta }_x)$ for each $s>0$.

(c)

Assume that $P_1\sim Dp(s_1,{\alpha }_1^{*})$, $P_2\sim Dp(s_2,{\alpha }_2^{*})$, $({\omega }_1,{\omega }_2)\sim Dir(s_1,s_2)$ and $P_1$, $P_2$, $({\omega }_1,{\omega }_2)$ are independent, then Section 3.1.1. in [17]:

$${\omega }_1{P}_1+{\omega }_2{P}_2\sim Dp\left(s_1+s_2,\frac{s_1}{s_1+s_2}{\alpha }_1^{*}+\frac{s_2}{s_1+s_2}{\alpha }_2^{*}\right).$$

(5)

(d)

Let $P_{X|X^n}$ have distribution $Dp(s+n,\frac{s}{s+n}{\alpha }^{*}+\frac{1}{s+n}{\sum }_{i=1}^n{\delta }_{X_i})$. We can write

$$P_{X|X^n}={\omega }_{0}P+\sum _{i=1}^n{\omega }_i{\delta }_{X_i},$$

(6)

where ${\sum }_{i=0}^n{\omega }_i=1$, $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$ and $P\sim Dp(s,{\alpha }^{*})$. This follows from (b)–(c).

An issue in the use of the DP as prior measure on P is how to choose the infinite dimensional parameter $G_0$ in case of lack of prior information. There are two avenues that we can follow. The first assumes that prior ignorance can be modelled satisfactorily by a so-called noninformative prior. In the DP setting, the only noninformative prior that has been proposed so far is the limiting DP obtained for $s\to 0$, which has been introduced by [16] and discussed by [18]. The second approach suggests that lack of prior information should be expressed in terms of a set of probability distributions. This approach known as Imprecise Probability [19,20,21] is connected to Bayesian robustness [22,23,24] and it has been extensively applied to model prior (near-)ignorance in parametric models. In this paper, we implement a prior (near-)ignorance model by considering a set of DPs obtained by fixing s to a strictly positive value and letting $G_0$ span the set of all distributions. This model has been introduced in [4] with the name of Imprecise Dirichlet Process (IDP).

3. Bayesian Independence Tests

Let us denote by X the vector of variables ${[Y,Z]}^T$ so that the n observations of X can be rewritten as

$$X^n=\{X_1,\cdots ,X_n\},$$

(7)

that is, a set of n vector-valued i.i.d. observations of X. We also consider an auxiliary variable $X^{\prime }$ together with X. We assume that $X,X^{\prime }$ are independent variables from the same unknown distribution and that ${X^{\prime }}^n=X^n$, that is, we have the same observations of X and $X^{\prime }$.

Let P be the unknown distribution of $X,X^{\prime }$ and assume that the prior distribution of P is $Dp(s,{\alpha }^{*})$. Our goal is to compute the posterior of P. The posterior of P is given in (3) and, by exploiting (6), we know that

$$P_{X|X^n}\sim {\omega }_{0}P+\sum _{i=1}^n{\omega }_i{\delta }_{X_i},$$

(8)

with $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$ and $P\sim Dp(s,{\alpha }^{*})$. The distribution $P_{X^{\prime }|{X^{\prime }}^n}$ of $X^{\prime }$ is similarly defined.

The questions we pose in a statistical analysis can all be answered by querying this posterior distribution in different ways. We adopt this posterior distribution to devise Bayesian counterparts of the independence hypothesis tests.

3.1. Bayesian Bivariate Independence Test for Binary Variables

Let us assume that the variables $Y,Z\in \{0,1\}$ (that is, they are binary). Our aim is to devise a Bayesian independence test for binary variables based on the DP. We will also show that our test is a Bayesian generalisation of the frequentist ${\chi }^2$-test for independence applied to binary variables. We start by defining the following quantities:

$$E\left(I_{\{0,0\}}(X)I_{\{1,1\}}(X^{\prime })|X^n,{X^{\prime }}^n\right)=\int \int I_{\{0,0\}}(X)I_{\{1,1\}}(X^{\prime })dF(X|X^n)dF(X^{\prime }|{X^{\prime }}^n),$$

where we have exploited the independence of $X,X^{\prime }$ and here $F(X|X^n)$ denotes the posterior cumulative distribution of $P_{X|X^n}$ defined in (8). From (8), it can easily be verified that

$$E\left(I_{\{0,0\}}(X)I_{\{1,1\}}(X^{\prime })|X^n,{X^{\prime }}^n\right)={\omega }_{00}{\omega }_{11},$$

where

$${\omega }_{00}={\omega }_0\int I_{\{0,0\}}(X)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{0,0\}}(Y_i,Z_i),$$

(9)

and

$${\omega }_{11}={\omega }_0\int I_{\{1,1\}}(X)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{1,1\}}(Y_i,Z_i),$$

(10)

where in the last equality we have exploited the fact that $X^{\prime }$ has the same distribution as X and also the same observations. The two quantities ${\omega }_{00},{\omega }_{11}$ include two terms. The first is the term due to the prior $dF\sim Dp(s,{\alpha }^{*})$ and the second term is due to the observations.

Similarly, we compute

$$E\left(I_{\{0,1\}}(X)I_{\{1,0\}}(X^{\prime })|X^n,{X^{\prime }}^n\right)={\omega }_{01}{\omega }_{10},$$

where

$${\omega }_{01}={\omega }_0\int I_{\{0,1\}}(X)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{0,1\}}(Y_i,Z_i),$$

(11)

and

$${\omega }_{11}={\omega }_0\int I_{\{1,0\}}(X)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{1,0\}}(Y_i,Z_i).$$

(12)

Summing up, ${\omega }_{00},{\omega }_{1,0},{\omega }_{0,1},{\omega }_{11}$ represent the posterior probabilities of the events $(0,0)$ (that is, $Y=0$ and $Z=0$), $(1,0)$, $(0,1)$ and $(1,1)$, respectively, according to the posterior joint distribution $F(X|X^n)$.

Theorem 1. 

The variables Y and Z are said to be concordant (dependent) with posterior probability $(1-\gamma )$ provided that

$$\mathcal{P}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)>(1-\gamma ),$$

(13)

and they are said to be discordant provided that

$$\mathcal{P}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})<0|X^n)>(1-\gamma ),$$

(14)

where $\mathcal{P}$ is the probability computed with respect to $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$ and $dF\sim Dp(s,{\alpha }^{*})$. Finally, they are said to be simply dependent with posterior probability $(1-\gamma )$ provided that

$$0\notin (1-\gamma )HDI(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})|X^n),$$

(15)

where $HDI$ denotes the posterior Highest Density Interval of $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$.

Proof. 

We just derive the third statement. The other two statements are analogue. We first consider the indicator functions

$$I_{\{0\}}(Y),I_{\{1\}}(Y),I_{\{0\}}(Z),I_{\{1\}}(Z),$$

(16)

and same for the auxiliary variables $Y^{\prime },Z^{\prime }$. By computing the expectation of these functions, we can obtain the marginals of the variables $Y,Z$ with respect to the joint $P_X$:

$$\begin{array}{c}\hfill {\omega }_{0•}:=E(I_{\{0\}}(Y)|X^n)={\omega }_0\int I_{\{0\}}(Y)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{0\}}(Y_i),\end{array}$$

(17)

$$\begin{array}{c}\hfill {\omega }_{1•}:=E(I_{\{1\}}(Y)|X^n)={\omega }_0\int I_{\{1\}}(Y)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{1\}}(Y_i),\end{array}$$

(18)

$$\begin{array}{c}\hfill {\omega }_{•0}:=E(I_{\{0\}}(Z)|X^n)={\omega }_0\int I_{\{0\}}(Z)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{0\}}(Z_i),\end{array}$$

(19)

$$\begin{array}{c}\hfill {\omega }_{•1}:=E(I_{\{1\}}(Z)|X^n)={\omega }_0\int I_{\{1\}}(Z)dF(X)+\sum _{i=1}^n{\omega }_i{I}_{\{1\}}(Z_i),\end{array}$$

(20)

where ${\omega }_{0•}$ (resp. ${\omega }_{1•}$) denotes the marginal with respect to Z when $Y=0$ (resp. $Y=1$), while ${\omega }_{•0}$ (resp. ${\omega }_{•1}$) denotes the marginal with respect to Y when $Y=0$ (resp. $Y=1$).

Then, by exploiting independence between X and $X^{\prime }$, we derive

$$\begin{array}{c}\hfill E(I_{\{0\}}(Y)I_{\{0\}}(Z^{\prime })|X^n,{X^{\prime }}^n)={\omega }_{0•}{\omega }_{•0},E(I_{\{1\}}(Y)I_{\{0\}}(Z^{\prime })|X^n,{X^{\prime }}^n)={\omega }_{1•}{\omega }_{•0},\end{array}$$

(21)

$$\begin{array}{c}\hfill E(I_{\{0\}}(Y)I_{\{1\}}(Z^{\prime })|X^n,{X^{\prime }}^n)={\omega }_{0•}{\omega }_{•1},E(I_{\{1\}}(Y)I_{\{1\}}(Z^{\prime })|X^n,{X^{\prime }}^n)={\omega }_{1•}{\omega }_{•1}.\end{array}$$

(22)

We are now ready to define the independence test. If the two variables $Y,Z$ are independent, then the vector

$$v=({\omega }_{00},{\omega }_{10},{\omega }_{01},{\omega }_{11})-({\omega }_{0•}{\omega }_{•0},{\omega }_{1•}{\omega }_{•0},{\omega }_{0•}{\omega }_{•1},{\omega }_{1•}{\omega }_{•1}),$$

has zero mean. Note that the first component of the vector v is $E(I_{\{0,0\}}(X)-I_{\{0\}}(Y)I_{\{0\}}(Z^{\prime })|X^n,{X^{\prime }}^n)$ and thus is a well-defined quantity with respect to our probabilistic model (similarly for the other terms). Therefore, the independence test reduces to checking whether the $(1-\gamma )\%$ highest density credible region (HCR) of v includes the zero vector. It can be easily verified that $|v_l|=|v_m|$ for each $l,m$ component of v. In fact, we have

$${\omega }_{i•}={\omega }_{ij}+{\omega }_{i\overline{j}},{\omega }_{•j}={\omega }_{ij}+{\omega }_{\overline{i}j},$$

for $i,j\in \{0,1\}$ and $\overline{i}=1-i$, $\overline{j}=1-j$, and so

$$\begin{array}{ccc}\hfill {\omega }_{ij}-({\omega }_{ij}+{\omega }_{i\overline{j}})({\omega }_{ij}+{\omega }_{\overline{i}j})& =& {\omega }_{ij}({\omega }_{ij}+{\omega }_{\overline{i}j}+{\omega }_{i\overline{j}}+{\omega }_{\overline{i}\overline{j}})-({\omega }_{ij}+{\omega }_{i\overline{j}})({\omega }_{ij}+{\omega }_{\overline{i}j})\hfill \\ & =& {\omega }_{ij}{\omega }_{\overline{i}\overline{j}}-{\omega }_{i\overline{j}}{\omega }_{\overline{i}j}.\hfill \end{array}$$

Therefore, it is enough to check whether

$$0\notin (1-\gamma )\%HDI(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})).$$

If this is the case, then we can declare that the two variables are dependent with probability $(1-\gamma )$. Here, the multiplier 2 in $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$ is only a scaling factor so that $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$ varies in $[-0.5,0.5]$. ☐

From the proof of Theorem 1 it is evident the similarity of the test with the frequentist ${\chi }^2$-test for independence. Both tests use the difference between the joint and the product of the marginals as a measure of dependence. The advantage of the Bayesian approach is that we compute posterior probabilities for the hypothesis in which we are interested and not the probability of getting the observed (or a larger) difference under the assumption that the null hypothesis of independence is true.

The probabilities computed in Theorem 1 depend on the prior information $F\sim Dp(s,{\alpha }^{*})$. In this paper we adopt IDP as prior model. We can then perform a Bayesian nonparametric test that is based on extremely weak prior assumptions, and easy to elicit, since it requires only the choice of the strength s of the DP instead of its infinite-dimensional parameter ${\alpha }^{*}$. The infinite-dimensional parameter ${\alpha }^{*}$ is free to vary in the set of all distributions.

Let us consider for instance (13). Each one of these priors gives a posterior probability $\mathcal{P}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$. We can characterize this set of posteriors by computing the lower and upper bounds $\underline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$ and $\overline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$. Inferences with IDP can be computed by verifying if

$$\underline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)>(1-\gamma ),\overline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)>(1-\gamma ),$$

(23)

and then by taking the following decisions:

• if both the inequalities are satisfied, then we declare that the two variables are dependent with probability larger than $1-\gamma$;
• if only one of the inequalities is satisfied (which has necessarily to be the one for the upper), we are in an indeterminate situation, that is, we cannot decide;
• if both are not satisfied, then we declare that the probability that the two variables are dependent is lower than the desired probability of $1-\gamma$.

When IDP returns an indeterminate decision, it means that the evidence from the observations is not enough to declare that the probability of the hypothesis being true is either larger or smaller than the desired value $1-\gamma$; more observations are necessary to reach a reliable decision.

Theorem 2. 

The upper probability $\overline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}=m{\delta }_{(0,0)}+(1-m){\delta }_{(1,1)}$ with

$$m=\left\{\begin{array}{cc}0& if{e}_1+{\omega }_0<e_0,\hfill \\ 1& if{e}_0<e_1-{\omega }_0,\hfill \\ {\displaystyle \frac{{\omega }_0+e_1-e_0}{2{\omega }_0}}& other,\hfill \end{array}\right.$$

(24)

where $e_0={\sum }_{i=1}^n{\omega }_i{I}_{\{0,0\}}(Y_i,Z_i)$ and $e_1={\sum }_{i=1}^n{\omega }_i{I}_{\{1,1\}}(Y_i,Z_i)$. The lower probability $\underline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}=m{\delta }_{(1,0)}+(1-m){\delta }_{(0,1)}$ with the same m as before but $e_0={\sum }_{i=1}^n{\omega }_i{I}_{\{1,0\}}(Y_i,Z_i)$ and $e_1={\sum }_{i=1}^n{\omega }_i{I}_{\{0,1\}}(Y_i,Z_i)$.

Proof. 

We are interested in the quantity $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$. It is clear that in order to maximize the probability that ${\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}>0$ we must put all the prior mass on ${\omega }_{00}{\omega }_{11}$. Let us call $m=I_{\{0,0\}}(X)dF(X)$. Then $\int I_{\{1,1\}}(X)dF(X)=1-I_{\{0,0\}}(X)dF(X)=1-m$. From (9)–(12), we have that ${\omega }_{00}={\omega }_{0}m+e_0$ and ${\omega }_{11}={\omega }_0(1-m)+e_1$ with $m\in [0,1]$. By computing the derivative with respect to m we have

$$\begin{array}{c}{\displaystyle \frac{d}{dm}}{\omega }_{00}{\omega }_{11}={\omega }_0\left({\omega }_0(1-m)+e_1\right)+\left({\omega }_{0}m+e_0\right)\left(-{\omega }_0\right),\hfill \end{array}$$

whose zero is $m={\displaystyle \frac{{\omega }_0+e_1-e_0}{2{\omega }_0}}$, which is also a maximum. Hence, the maximum can be either on $m={\displaystyle \frac{{\omega }_0+e_1-e_0}{2{\omega }_0}}$ or on the extremes $m=0$ or $m=1$. This can be easily verified by checking when ${\omega }_0+e_1-e_0<0$ (so $m=0$) or ${\omega }_0+e_1-e_0>2{\omega }_0$ (so $m=1$). The lower probability $\underline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$ can be determined using a similar reasoning. ☐

Since $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$, the computation of $\underline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$, $\overline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)$ can be obtained by Monte Carlo sampling. The following pseudo-code describes how to compute the upper (the lower can be computed in a similar way).

• Initialize the counter $P_c$ to 0 and the array V to empty;
• For $i=1,\cdots ,N_{mc}$

  (a)

  sample $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$;

  (b)

  compute ${\omega }_{00},{\omega }_{01},{\omega }_{10},{\omega }_{11}$ as in (9)–(12) by choosing $dF(X)=m{\delta }_{(0,0)}(X)+(1-m){\delta }_{(1,1)}(X)$ with m defined in Theorem 2;

  (c)

  compute $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$ and store the result in V;

  (d)

  if $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0$ then $P_c=P_c+1$ else $P_c=P_c+0$.

• compute the histogram of the elements in V (this gives us the plot of the posterior of $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$)
• compute the posterior upper probability that $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$ is greater than zero as $\overline{\mathcal{P}}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})>0|X^n)\approx P_c/N_{mc}$.

The number of Monte Carlo samples $N_{mc}$ is equal to 100 thousand in the next examples and figures.

The lower and upper $HDI$ intervals in Theorem 1 can also be obtained as in Theorem 2 and computed via Monte Carlo sampling ($HDI$ can be computed using the values stored in V, see pseudo-code). Hereafter we will denote the two intervals corresponding to the lower and upper distributions as $\underline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$ and $\overline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$, respectively.

The only prior parameter that must be selected with IDP is the prior strength s. The value of s determines how quickly the posteriors corresponding to the lower and upper probabilities converge as the number of observations increases. We select $s=0.5$—this means that we need at least 4 concordant binary observations to take a decision with $1-\gamma =0.95$. In other words, for $s=0.5$ we need two observations of type $Y=0,Z=0$ and two of type $Y=1,Z=1$ to guarantee that both $1-\gamma =0.95\%$ $HDI$ intervals, i.e., $\underline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$ and $\overline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$, do not include the zero. For any number of (and configuration of) observations less than four, the test is always indeterminate (i.e., no decision can be taken). Thus, four is the minimum number of observations that is required to take a decision. This choice is arbitrary and subjective, but is our measure of cautiousness. We make clearer the meaning of determinate and indeterminate in the following example.

Example 1. 

Let us consider the following three matrices of 10 paired binary i.i.d. observations

$$\begin{array}{c}{X}_a^{10}={\left[\begin{array}{cccccccccc}0& 0& 0& 0& 0& 1& 1& 1& 1& 1\\ 0& 0& 0& 0& 0& 1& 1& 1& 1& 0\end{array}\right]}^T,\\ X_b^{10}={\left[\begin{array}{cccccccccc}0& 0& 0& 0& 0& 1& 1& 1& 1& 1\\ 0& 0& 0& 0& 0& 1& 1& 0& 0& 0\end{array}\right]}^T,\\ X_c^{10}={\left[\begin{array}{cccccccccc}0& 1& 1& 0& 0& 1& 1& 1& 0& 1\\ 0& 0& 0& 1& 0& 1& 1& 1& 1& 1\end{array}\right]}^T.\end{array}$$

(25)

They correspond to different degrees of dependence. Figure 1 shows the lower and upper distributions of $2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10})$ and the relative $95\%$ HDI, i.e., $\underline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$ and $\overline{HDI}(2({\omega }_{00}{\omega }_{11}-{\omega }_{01}{\omega }_{10}))$, for the three cases $a,b,c$ (the filled in areas). In case (a), the two variables are dependent (concordant) with probability greater than $0.95$, since all the mass of the lower and upper distributions are in the interval $[0,0.5]$. In the second case, we are in an indeterminate situation, that is, the lower and upper are in disagreement, which means that the inference is prior dependent. In the third case, we can only say that they are not dependent at $95\%$ since both the $HDI$ intervals include the zero.

3.2. Bayesian Bivariate Independence Test for Continuous Variables

Let us assume that variables $Y,Z\in \mathbb{R}$, that is, they are real continuous variables. Our aim is to devise a Bayesian independence test for continuous variables based on the DP. We will also show that our test is a Bayesian generalisation of Kendall-τ test for independence. This test uses results from [25] that derived a Bayesian Kendall’s τ statistic using DP. As before, we introduce auxiliary variables $Y^{\prime },Z^{\prime }$. We start by defining the following quantities:

$$\begin{array}{ccc}\hfill T_1& =& \{(Y,Z,Y^{\prime },Z^{\prime }):(Y-Y^{\prime })(Z-Z^{\prime })>0\},\hfill \\ \hfill T_2& =& \{(Y,Z,Y^{\prime },Z^{\prime }):(Y-Y^{\prime })(Z-Z^{\prime })<0\}.\hfill \end{array}$$

$T_1$ and $T_2$ are concordance measures. We can then compute

$$E[I_{T_1}-I_{T_2}]=\int \int (I_{T_1}(X,X^{\prime })-I_{T_2}(X,X^{\prime }))dF(X|X^n)dF(X^{\prime }|{X^{\prime }}^n),$$

(26)

where we have exploited the independence of $X,X^{\prime }$ and here $F(X|X^n)$ denotes the posterior cumulative distribution of $P_{X|X^n}$. This quantity is equal to

$$\begin{array}{cc}\hfill E[I_{T_1}-I_{T_2}]& ={\omega }_0^2\int \int (I_{T_1}(X,X^{\prime })-I_{T_2}(X,X^{\prime }))dF(X)dF(X^{\prime })\hfill \\ \hfill & +2\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{T_1}(X_i,X^{\prime })-I_{T_2}(X_i,X^{\prime }))dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{T_1}(X_i,X_j)-I_{T_2}(X_i,X_j)),\hfill \end{array}$$

where we have exploited the fact that $X^{\prime }$ has the same distribution as X and the same observations. Given $({\omega }_0,\cdots ,{\omega }_n)$, it can be seen that the first two terms depend on the prior distribution $F\sim Dp(s,{\alpha }^{*})$ and the last term is only due to the observations.

Theorem 3. 

The variables Y and Z are said to be concordant (dependent) with posterior probability $(1-\gamma )$ provided that

$$\mathcal{P}(E[I_{T_1}-I_{T_2}]/2>0|X^n)>(1-\gamma ),$$

(27)

and they are said to be discordant provided that

$$\mathcal{P}(E[I_{T_1}-I_{T_2}]/2<0|X^n)>(1-\gamma ),$$

(28)

where $\mathcal{P}$ is the probability computed with respect to $({\omega }_0,{\omega }_1,\cdots ,{\omega }_n)\sim Dir(s,1,\cdots ,1)$ and $dF\sim Dp(s,{\alpha }^{*})$. Finally, they are said to be simply dependent with posterior probability $(1-\gamma )$ provided that

$$0\notin (1-\gamma )HDI(E[I_{T_1}-I_{T_2}]/2|X^n),$$

(29)

where $HDI$ denotes the posterior Highest Density Interval of $E[I_{T_1}-I_{T_2}]/2$.

The divisor 2 in $E[I_{T_1}-I_{T_2}]/2$ is only a scaling factor so that the expectation lies in $[-0.5,0.5]$. The theorem simply follows from the fact that $E[I_{T_1}-I_{T_2}]$ is the same measure of dependence used in Kendall’s τ test. In this respect, it is worth to highlight the connection with Kendall’s τ. By exploiting the properties of DP, we have that the posterior mean of $E[I_{T_1}-I_{T_2}]$ for large n is approximately equal to.

$$\mathcal{E}(E[I_{T_1}-I_{T_2}]|X^n)\approx {\displaystyle \frac{1}{(n+1)n}}\sum _{i=1}^n\sum _{j=1}^n(I_{T_1}(X_i,X_j)-I_{T_2}(X_i,X_j))$$

(30)

and this is exactly Kendall’s sample τ coefficient. In fact, Kendall’s sample τ coefficient is defined as:

$$T=2\sum _{1\le i<j\le n}{\displaystyle \frac{A_{ij}}{n(n-1)}}=2\sum _{i=1}^{n-1}\sum _{j=i+1}^n{\displaystyle \frac{A_{ij}}{n(n-1)}},$$

(31)

with

$$A_{ij}=\left\{\begin{array}{cc}1,\hfill & \mathit{if}(Y_i-Y_j)(Z_i-Z_j)>0,\hfill \\ -1,\hfill & \mathit{if}(Y_i-Y_j)(Z_i-Z_j)<0.\hfill \end{array}\right.$$

Observe that T can also be rewritten as:

$$T=\sum _{i=1}^n\sum _{j=1}^n{\displaystyle \frac{A_{ij}}{n(n-1)}},$$

(32)

in terms of all the $A_{ij}$ pairs, which is proportional to (30) for large n. This clarifies the connection between our Bayesian test of dependence for continuous variables based on $E[I_{T_1}-I_{T_2}]/2$ and Kendall’s τ test.

As for the dependence test for binary variables, we will make inferences using IDP. Inferences with IDP can computed by verifying if

$$\underline{\mathcal{P}}(E[I_{T_1}-I_{T_2}]/2>0|X^n)>(1-\gamma ),\overline{\mathcal{P}}(E[I_{T_1}-I_{T_2}]/2>0|X^n)>(1-\gamma ).$$

(33)

Theorem 4. 

The upper probability $\overline{\mathcal{P}}(E[I_{T_1}-I_{T_2}]/2>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}=0.5{\delta }_{X_0^a}+0.5{\delta }_{X_0^b}$ with $X_0^a>X_0^b>X_i$ for $i=1,\cdots ,n$. The lower probability $\underline{\mathcal{P}}(E[I_{T_1}-I_{T_2}]/2>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}=0.5{\delta }_{X_0^a}+0.5{\delta }_{X_0^b}$ with $X_0=(Y_0,Z_0)$ and $Y_0^a>Y_0^b>Y_i$ and $Z_0^a<Z_0^b<Z_i$ for $i=1,\cdots ,n$.

Proof. 

We have that

$$\begin{array}{cc}\hfill E[I_{T_1}-I_{T_2}]& ={\omega }_0^2\int \int (I_{T_1}(X,X^{\prime })-I_{T_2}(X,X^{\prime }))dF(X)dF(X^{\prime })\hfill \\ \hfill & +2\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{T_1}(X_i,X^{\prime })-I_{T_2}(X_i,X^{\prime }))dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{T_1}(X_i,X_j)-I_{T_2}(X_i,X_j)).\hfill \end{array}$$

We want to maximize $I_{T_1}(X,X^{\prime })$. Since $\int \int I_{T_1}(X,X^{\prime }){\delta }_{X_0^a}(X){\delta }_{X_0^a}(X^{\prime })dXd{X}^{\prime }=0$, we need at least two Dirac’s deltas. Hence, we consider the mixture $dF=m{\delta }_{X_0^a}+(1-m){\delta }_{X_0^b}$ with $X_0^a>X_0^b>X_i$ for $i=1,\cdots ,n$. Then we have that

$$\begin{array}{cc}\hfill E[I_{T_1}-I_{T_2}]& =m(1-m){\omega }_0^2+2\sum _{i=1}^n{\omega }_0{\omega }_i+\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{T_1}(X_i,X_j)-I_{T_2}(X_i,X_j)),\hfill \end{array}$$

and so we have maximized the second term. For the first term depending on $m(1-m)$, the maximum is obtained at $m=1/2$. For the lower probability, the proof is similar. ☐

The lower and upper $HDI$ intervals can also be obtained as in Theorem 4. Again in this case, the value of s determines how quickly lower and upper posteriors converge as the number of observations increases. We choose $s=0.5$ as for the binary test.

Example 2. 

Also in this case we consider three matrices of 10 paired continuous i.i.d. observations

$$\begin{array}{c}{X}_a^{10}={\left[\begin{array}{cccccccccc}-0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& 0.3& 0.2& 0.1\\ -0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& -0.3& -0.2& -0.1\end{array}\right]}^T,\\ X_b^{10}={\left[\begin{array}{cccccccccc}-0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& 0.3& 0.2& 0.1\\ -0.1& -0.2& -0.3& -0.4& -0.5& 0.5& -0.4& -0.3& -0.2& -0.1\end{array}\right]}^T,\\ X_c^{10}={\left[\begin{array}{cccccccccc}-0.1& 0.2& 0.3& -0.4& -0.5& -0.5& 0.4& 0.3& 0.2& 0.1\\ 0.1& -0.2& -0.3& -0.4& 0.5& 0.5& 0.4& -0.3& -0.2& -0.1\end{array}\right]}^T.\end{array}$$

(34)

They correspond to different degrees of dependence. Figure 2 shows the lower and upper posteriors for the three cases $a,b,c$ and the relative $HDI$ intervals at $95\%$ probability (the filled in areas). In case (a), the two variables are dependent (concordant) with probability greater than $0.95$, since all the mass of the lower and upper distributions are in the interval $[0,0.5]$. In the second case, we are in an indeterminate situation, that is, the lower and upper are in disagreement, which means that the inference is prior dependent. In the third case, we can only say that they are not dependent at $95\%$ since both the $HDI$ intervals include the zero.

3.3. Bayesian Bivariate Independence Test for Mixed Continuous-Binary Variables

Let us assume that the variables $Y\in \mathbb{R}$ and $Z\in \{0,1\}$. Our aim is to devise a Bayesian independence test based on the DP. We introduce the auxiliary variable $X^{\prime }$ as done before. To derive our test, we start by defining the following indicator:

$$I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime }).$$

This indicator is one if $X=(Y,0)$ and $X^{\prime }=(Y^{\prime },1)$, with $Y>Y^{\prime }$ and zero otherwise. We can compute

$$E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]=\int \int (I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime }))dF(X|X^n)dF(X^{\prime }|{X^{\prime }}^n),$$

(35)

where we have exploited the independence of $X,X^{\prime }$. $F(X|X^n)$ denotes the posterior cumulative distribution of $P_{X|X^n}$. This quantity is equal to

$$\begin{array}{cc}\hfill E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]& ={\omega }_0^2\int \int (I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime }))dF(X)dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{(Y^{\prime },\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z^{\prime }))dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{(Y_i,\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z_i))dF(X)\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{(Y_j,\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z_j)).\hfill \end{array}$$

For large n, we have that

$$\mathcal{E}(E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]|X^n)\approx \frac{1}{(n+1)n}\sum _{i=1}^n\sum _{j=1}^n{I}_{(Y_j,\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z_j),$$

(36)

which is equal to the rank of Y in the observations $(Y,0)$ with respect to the observations $(Y,1)$. Therefore, our dependence test is rank-based. It is clear that, in the case of independence of the variables Y and Z, the mean rank is equal to $0.125$. Hence, we can formulate an independence test for mixed variables.

Theorem 5. 

The variables Y and Z are dependent with posterior probability $(1-\gamma )$ provided that

$$0\notin (1-\gamma )HDI(4E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]-0.5|X^n),$$

(37)

where $HDI$ denotes the posterior Highest Density Interval of $4E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]-0.5$.

The theorem follows from the fact that in case of independence between variables Y and Z the mean rank (36) scaled by 4 and shifted of $-0.5$ is equal to 0. Also in this case, we make inferences using IDP.

Theorem 6. 

The upper probability $\overline{\mathcal{P}}(4E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]-0.5>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}=m{\delta }_{X_0^a}+(1-m){\delta }_{X_0^b}$ with $X_0^a$ equal to $(-\infty ,1)$ and $X_0^b$ equal to $(\infty ,0)$ and

$$m=\left\{\begin{array}{cc}0& if{\omega }_0+e_0<e_1,\hfill \\ 1& if{e}_1<e_0-{\omega }_0,\hfill \\ {\displaystyle \frac{{\omega }_0+e_0-e_1}{2{\omega }_0}}& other,\hfill \end{array}\right.$$

with $e_0={\sum }_{i=1}^n{\omega }_i{I}_{\{0\}}(Z_i)$ and $e_1={\sum }_{i=1}^n{\omega }_i{I}_{\{1\}}(Z_i)$. The lower probability $\underline{\mathcal{P}}(4E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]-0.5>0|X^n)$ is obtained by a prior measure ${\alpha }^{*}={\delta }_{X_0}$ with $X_0$ equal to $(-\infty ,0)$.

Proof. 

Consider the quantity

$$\begin{array}{cc}\hfill E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]& ={\omega }_0^2\int \int (I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime }))dF(X)dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{(Y^{\prime },\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z^{\prime }))dF(X^{\prime })\hfill \\ \hfill & +\sum _{i=1}^n{\omega }_0{\omega }_i\int (I_{(Y_i,\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z_i))dF(X)\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{(Y_j,\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z_j)),\hfill \end{array}$$

and ${\alpha }^{*}=m{\delta }_{X_0^a}+(1-m){\delta }_{X_0^b}$ with $X_0^a$ equal to $(-\infty ,1)$ and $X_0^b$ equal to $(\infty ,0)$. Thus

$$\begin{array}{cc}\hfill E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]& =m(1-m){\omega }_0^2+m\sum _{i=1}^n{\omega }_0{\omega }_i{I}_{\{0\}}(Z_i)+(1-m)\sum _{i=1}^n{\omega }_0{\omega }_i{I}_{\{1\}}(Z_i))\hfill \\ \hfill & +\sum _{i=1}^n\sum _{j=1}^n{\omega }_i{\omega }_j(I_{(Y_j,\infty )}(Y_i)I_{\{0\}}(Z_i)I_{\{1\}}(Z_j).\hfill \end{array}$$

By computing the derivative

$${\displaystyle \frac{d}{dm}}E[I_{(Y^{\prime },\infty )}(Y)I_{\{0\}}(Z)I_{\{1\}}(Z^{\prime })]={\omega }_0-2{\omega }_{0}m+e_0-e_1=0,$$

we have that $m=\frac{{\omega }_0+e_0-e_1}{2{\omega }_0}$. The result is obtained by exploiting the fact that $m\in [0,1]$. For the lower probability, the computation is straightforward. ☐

The lower and upper $HDI$ intervals can also be obtained as in Theorem 4. We choose $s=0.5$ as for the previous tests.

Example 3. 

We consider three matrices of 10 paired binary-continuous i.i.d. observations

$$\begin{array}{c}{X}_a^{10}={\left[\begin{array}{cccccccccc}-0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& 0.3& 0.2& 0.1\\ 1& 1& 1& 1& 1& 0& 0& 0& 0& 0\end{array}\right]}^T,\\ X_b^{10}={\left[\begin{array}{cccccccccc}-0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& 0.3& 0.2& 0.1\\ 1& 1& 1& 1& 1& 0& 0& 0& 0& 1\end{array}\right]}^T,\\ X_c^{10}={\left[\begin{array}{cccccccccc}-0.1& -0.2& -0.3& -0.4& -0.5& 0.5& 0.4& 0.3& 0.2& 0.1\\ 1& 0& 0& 1& 1& 0& 0& 1& 0& 1\end{array}\right]}^T.\end{array}$$

(38)

Again, they correspond to different degrees of dependence. Figure 3 shows the lower and upper posteriors for the three cases $a,b,c$ and the relative $HDI$ intervals at $95\%$ probability (the filled in areas). In case (a), the two variables are dependent (concordant) with probability greater than $0.95$, since all the mass of the lower and upper distributions are in the interval $[0,0.5]$. In the second case, we are in an indeterminate situation, that is, the lower and upper are in disagreement. In the third case, we can only say that they are not dependent at $95\%$ since both the $HDI$ intervals include the zero.

4. Experiments

We compare our Bayesian testing approach in the three discussed main scenarios where both variables are binary, both are continuous and one is binary and the other is continuous. The goal is to decide whether the two variables are dependent or independent. We generate n samples ($n=20$ and 50) using the distributions defined in

Table 1. Data generation setup. In order to generate independent data, β is set to zero. Larger values of β increase their dependency.

Variable 1	Variable 2	Distribution
Binary	Binary	Multinomial distr. with $[P(00),P(01),P(10),P(11)]\propto [3,3+\beta ,3+\beta ,3]$.
Continuous	Continuous	Bivariate Gaussian with means 0 and covariance matrix $\left[\begin{array}{cc}10& \beta \\ \beta & 3\end{array}\right]$.
Binary	Continuous	Half of the samples have the binary variable set to zero and half to one. When that variable is zero, then for the continuous use $\mathrm{\Gamma }(10,2)$, otherwise $\mathrm{\Gamma }(10+\beta ,2+\beta )$.

. Ten thousand repetitions are used by forcing the variables to be independent (so $\beta =0$) and thousand repetitions where the variables are dependent, for each value of $\beta >0$. The value of β is varied as explained in the table. For each n, β and each of these twenty thousand samples (for which we know the correct result of the test), we run the new approach versus ${\chi }^2$ test, Kendall τ test and Kolmogorov–Smirnov test, respectively for the binary-binary, continuous-continuous and binary-continuous cases. For each run of each method, we record their p-values, while for the new approach we compute γ corresponding to the limiting credible region $1-\gamma$ wide where the decision changes between dependent and independent. Such value is related to the p-values of the other tests and can be used for decision making by comparing it against a threshold (just as it is done with the p-values). However, it should be observed that thresholds different from $0.05$ or $0.01$ are hardly used in practice in null hypothesis significance tests. Conversely, for a Bayesian tests $1-\gamma$ is a probability and, therefore, we can take decisions with probability $0.99$, $0.95$ but also $0.7$ or even $0.51$ depending on the application (and the loss function). However, instead of fixing a threshold (which is a subjective choice) to decide between the options dependent and non-dependent with probability $1-\gamma$, we use Receiver Operating Characteristic (ROC) curves. ROC curves give the quality of the approaches for all possible thresholds. The curves are calculated as usual by varying the threshold from 0 to 1 and computing the sensitivity (or true positive rate) and specificity (or one minus false positive rate) (this is slightly different from the common approach of drawing ROC curves as a function of the true positive rate and false positive rate [26,27,28]). ROC curves are always computed considering different degrees of dependence (different values for $\beta \ne 0$) against independence ($\beta =0$). We apply the same criterion to p-values for comparing the methods across a wider range of decision criteria. We have used the R package “pROC” to compute the ROC curves [29].

Figure 4, Figure 5 and Figure 6 present the comparison of the new approach (which we name as IBinary, ICont or IMixed to explicitly account for the types of variables been analyzed) using $s\approx 0$ against the appropriate competitor. With such choice of s, the new approach runs without indeterminacy and can be directly compared against usual methods. As we see in the figures, the new method performs very similar to each competitor, with the advantage of being compatible among different types of data (the p-values of the other methods, among different data types, cannot be compared to each other). This is useful when one works with multivariate models involving multiple data types. As expected, the quality of the methods increases with the increase of β and of the sample size.

Figure 7, Figure 8 and Figure 9 present the ROC curves for the methods ${\chi }^2$, Kendall τ and Kolmogorov–Smirnov, respectively. These curves are separated according to whether the instance is considered determinate or indeterminate by the new approach. In other words, for each one of the twenty thousand repetitions, we run the corresponding usual test and then we check whether the output of the new approach is determinate or indeterminate (applying $s=0.5$), and we split the instances accordingly (blue curves show the accuracy over instances that are considered easy (determinate cases) while green curves over instances that are hard (indeterminate cases)—we also present the overall accuracy of the method using red curves). As we see, such division is able to identify easy-to-classify and hard-to-classify cases, since the ROC curves for the cases deemed as indeterminate by the new approach suggest a performance not better than a random guess (green curves). using the new approach, This means that if we would devise another test (called “50/50 when indeterminate”) which returns the same response as IBinary, ICont or IMixed when they are determinate, and issues a random answer (with 50/50 chance) otherwise, then this “50/50 when indeterminate” test would have the same ROC curve as ${\chi }^2$, Kendall τ and Kolmogorov–Smirnov, respectively.

This suggests that the indeterminacy of IDP based tests is an additional useful information that our approach gives to the analyst. In these cases she/he knows that (i) her/his posterior decisions would depend on the choice of the prior DP measure; (ii) deciding between the two hypotheses under test is a difficult problem as shown by the comparison with the DP with $s=0$, ${\chi }^2$, Kendall τ and Kolmogorov–Smirnov. Based on this additional information, the analyst can for example decide to collect additional measurements to eliminate the indeterminacy (in fact when the number of observations goes to infinity the indeterminacy goes to zero).

This represents a second advantage of our IDP approach, once we have fixed the value of s (e.g., $s=0.5$) it can automatically identify the risky cases where a decision must be taken with additional caution. For this reason, we suggest to use the IDP based test for dependence and not $s=0$.

Finally,

Table 2. Area under the ROC curve (AUC) values for all the performed experiments using different values of s, β and n. IBinary shows the AUC for the new test applied to two binary variables and $s\approx 0$. The columns ${\chi }^2$ test, Det.cases, and Indet.cases show the AUC obtained by the ${\chi }^2$ test over all samples, only over samples considered determinate by IBinary (with the corresponding s) and finally only over samples considered indeterminate by IBinary.

s	n	β	$\mathit{IBinary}$	Chisq	Det.cases	Indet.cases
0.25	20	1	0.5562	0.5629	0.5653	0.4890
0.5	20	1	0.5544	0.5596	0.5645	0.5233
1	20	1	0.5491	0.5551	0.5642	0.5153
0.25	20	3	0.7341	0.7502	0.7567	0.4266
0.5	20	3	0.7388	0.7551	0.7686	0.4526
1	20	3	0.7330	0.7502	0.7717	0.4888
0.25	50	1	0.6372	0.6425	0.6449	0.5125
0.5	50	1	0.6319	0.6353	0.6393	0.4747
1	50	1	0.6366	0.6407	0.6492	0.4954
0.25	50	3	0.9145	0.9110	0.9127	0.5205
0.5	50	3	0.9130	0.9090	0.9115	0.4473
1	50	3	0.9134	0.9081	0.9123	0.5642

,

Table 3. Area under the ROC curve (AUC) values for all the performed experiments using different values of s, β and n. ICont shows the AUC for the new test applied to two continuous variables and $s\approx 0$. Kendall, Det.cases, and Indet.cases show the AUC obtained by Kendall’s test over all samples, only over samples considered determinate by ICont (with the corresponding s) and finally only over samples considered indeterminate by ICont.

s	n	β	$\mathit{ICont}$	Kendall	Det.cases	Indet.cases
0.25	20	1	0.5826	0.5858	0.5898	0.5101
0.5	20	1	0.5708	0.5729	0.5804	0.4987
1	20	1	0.5744	0.5742	0.5914	0.5004
0.25	20	2	0.7524	0.7506	0.7558	0.5037
0.5	20	2	0.7535	0.7502	0.7574	0.5203
1	20	2	0.7488	0.7407	0.7596	0.5447
0.25	50	1	0.6825	0.6888	0.6917	0.5051
0.5	50	1	0.6782	0.6869	0.6935	0.5633
1	50	1	0.6871	0.6960	0.7087	0.5204
0.25	50	2	0.9343	0.9191	0.9197	0.4933
0.5	50	2	0.9339	0.9208	0.9207	0.5487
1	50	2	0.9361	0.9205	0.9192	0.5499

and

Table 4. Area under the ROC curve (AUC) values for all the performed experiments using different values of s, β and n. IMixed shows the AUC for the new test applied to one binary and one continuous variables and $s\approx 0$. Kolmogorov–Smirnov (KS), Det.cases, and Indet.cases show the AUC obtained by KS test over all samples, only over samples considered determinate by IMixed (with the corresponding s) and finally only over samples considered indeterminate by IMixed.

s	n	β	$\mathit{IMixed}$	KS	Det.cases	Indet.cases
0.25	20	1	0.6159	0.6118	0.6139	0.5386
0.5	20	1	0.6150	0.5943	0.5989	0.5594
1	20	1	0.6132	0.6004	0.6104	0.5532
0.25	20	2	0.7176	0.7358	0.7392	0.5254
0.5	20	2	0.7202	0.7091	0.7159	0.4937
1	20	2	0.7163	0.7091	0.7233	0.4928
0.25	50	1	0.6997	0.7091	0.7109	0.4447
0.5	50	1	0.6966	0.7106	0.7149	0.4213
1	50	1	0.7076	0.7135	0.7224	0.4455
0.25	50	2	0.8526	0.8816	0.8832	0.3278
0.5	50	2	0.8497	0.8790	0.8818	0.3044
1	50	2	0.8562	0.8923	0.8986	0.2934

present the values for the Area under the curve (AUC) in Chaper 5 in [30] of the ROC curves discussed previously, as well as similar experimental setup but with different values of s: 0.25, 0.5 and 1. Table 2 has results for binary variable versus binary variable, Table 3 for continuous variable versus continuous variable, and Table 4 for continuous variable versus binary variable. Overall, results show that IBinary has similar performance as ${\chi }^2$ test, ICont has similar performance as Kendall’s τ test and IMixed is similar to Kolmogorov–Smirnov (KS) test. The most interesting outcome is the comparison, in each scenario, of the frequentist test over whole data, over only data samples that were considered determinate by the new test, and over only data samples that were considered indeterminate. We clearly see that the AUC values over the cases considered indeterminate are much inferior to the values over cases considered determinate, which indicates that the new test has a good ability to discriminate easy and hard cases. ROC curves for values of s other than 0.5 were omitted for clarity of exposition, but they are very similar to those obtained for $s=0.5$.

5. Conclusions

We have proposed three novel Bayesian methods for performing independence tests for binary, continuous and mixed binary-continuous variables. All of these tests are nonparametric and based on the Dirichlet Process. This has allowed us to use the same prior model for all the tests we have developed. Therefore, all the tests are “consistent”, in the sense that the probabilities of dependence we compute with these tests are commensurable across the tests.

We have presented two versions of these tests: one based on a noninformative prior and one based on a conservative model of prior ignorance (IDP). Experimental results show that the prior ignorance method is more reliable than both the frequentist test and the noninformative Bayesian one, being able to isolate instances in which these tests are almost guessing at random. For future work, we plan to extend this approach in two directions: (1) feature selection in classification; (2) learning the structure (graph) of Bayesian networks and Markov Random Fields. The idea is to use our dependence tests to replace the frequentist tests that are commonly used for that purpose and evaluate the gain in terms of performance. For instance in case (1), we then could compare the accuracy of a classifier whose features are selected using our tests with that of a classifier whose features are selected by using frequentist tests. Our new approach is suitable since it addresses two limitations of currently used tests: they are based on null-hypothesis significance tests, and they cannot be applied to categorical and continuous variables at the same time in a commensurable way.