Matthieu Marbac

Mohammed Sedki

A Family of Blockwise One-Factor Distributions for Modelling High-Dimensional Binary Data

Keywords: Binary data, EM algorithm, High-dimensional data, IFM procedure, Model selection, One-factor copulas

We introduce a new family of one factor distributions for high-dimensional binary data. The model provides an explicit probability for each event, thus avoiding the numeric approximations often made by existing methods. Model interpretation is easy since each variable is described by two continuous parameters (corresponding to its marginal probability and to its strength of dependency with the other variables) and by one binary parameter (defining if the dependencies are positive or negative). An extension of this new model is proposed by assuming that the variables are split into independent blocks which follow the new one factor distribution. Parameter estimation is performed by the inference margin procedure where the second step is achieved by an expectation-maximization algorithm. Model selection is carried out by a deterministic approach which strongly reduces the number of competing models. This approach uses a hierarchical ascendant classification of the variables based on the empirical version of Cramer's V for selecting a narrow subset of models. The consistency of such procedure is shown. The new model is evaluated on numerical experiments and on a real data set. The procedure is implemented in the R package MvBinary available on CRAN.

Introduction

Binary data are increasingly emerging in various research fields, particularly in economics, psychometrics or in life sciences [START_REF] Cox | Analysis of binary data[END_REF][START_REF] Collett | Modelling binary data[END_REF]. To carry out statistical inference, it is important to have at hand flexible distributions for such data. However, there is a shortage of multivariate distributions for binary data [START_REF] Genest | A primer on copulas for count data[END_REF]. Indeed, many approaches have been developed by considering that the binary variables are responses of several explanatory variables [START_REF] Glonek | Multivariate logistic models[END_REF][START_REF] Nikoloulopoulos | Multivariate logit copula model with an application to dental data[END_REF][START_REF] Genest | Predicting dependent binary outcomes through logistic regressions and meta-elliptical copulas[END_REF]. However, these models cannot manage data composed with only binary variables.

Since binary variables are easily accessible and poorly discriminative, the binary data sets are often composed of many variables. Thus, the modelling of high-dimensional binary data is an important issue. Moreover, classical models suffer from the curse of dimensionality since they involve too many parameters [START_REF] Bellman | Dynamic Programming[END_REF]. Therefore, specific distributions should be introduced to manage such data.

Many authors have been interested in defining the properties of a multivariate distribution which permit an easy interpretation and inference [START_REF] Nikoloulopoulos | Finite normal mixture copulas for multivariate discrete data modeling[END_REF][START_REF] Panagiotelis | Pair copula constructions for multivariate discrete data[END_REF]. Thus, [START_REF] Nikoloulopoulos | Copula-based models for multivariate discrete response data[END_REF] lists the five following features that define a distribution with good properties: (F1) Wide range of dependence, allowing both positive and negative dependence; (F2) Flexible dependence, meaning that the number of bivariate marginals is (approximately) equal to the number of dependence parameters; (F3) Computationally feasible cumulative distribution function for likelihood estimation; (F4) Closure property under marginalization, meaning that lower-order marginals belong to the same parametric family; (F5) No joint constraints for the dependence parameters, meaning that the use of covariate functions for the dependence parameters is straightforward.

The modelling by dependency trees [START_REF] Chow | Approximating discrete probability distributions with dependence trees[END_REF]) is a pioneer approach for assessing the distribution of binary variables. A strength of this method is the easy maximization of likelihood function by the Kruskal algorithm [START_REF] Kruskal | On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF], which estimates the tree of minimal length. Although the tree structure leads to benefits (estimation, visualisation and interpretation), it is limited to simple dependency relations. Moreover, it does not provide parameters for measuring the strength of the dependencies between two variables.

A naïve approach for modelling binary variables is to use a product of Bernoulli distributions. However, in spite of the parsimony induced by the independence assumption, this approach leads to severe biases when variables are dependent. Thus, a mixture model with conditional independence assumption can capture the main dependencies [START_REF] Goodman | Exploratory latent structure analysis using both identifiable and unidentifiable models[END_REF]. [START_REF] Celeux | Clustering criteria for discrete data and latent class models[END_REF] propose different parsimonious models to deal with high-dimensional data. However, this mixture-based method suffers primarily from a lack of interpretation of dependencies. Indeed, there is no parameters for directly reflecting the strength of the dependencies between variables.

The quadratic exponential binary distribution [START_REF] Cox | The analysis of multivariate binary data[END_REF] is considered as the binary version of the multivariate Gaussian distribution. However, this model does not retain its exact form under marginalization, but closure under marginalization can be achieved approximately [START_REF] Cox | A note on the quadratic exponential binary distribution[END_REF]. This model is not really suitable for highdimensional data since it implies a quadratic number of parameters.

The modelling of spatial binary data can be achieved by latent Gaussian Markov Random Fields [START_REF] Pettitt | A conditional autoregressive Gaussian process for irregularly spaced multivariate data with application to modelling large sets of binary data[END_REF][START_REF] Weir | Binary probability maps using a hidden conditional autoregressive Gaussian process with an application to Finnish common toad data[END_REF] or lattice-based Markov Random Fields like the Ising model [START_REF] Gaetan | Spatial statistics and modeling[END_REF]. These approaches can deal with high-dimensional data since they have the Markov properties. However, their using for non-spatial data is not really doable except when the model at hand is known. Indeed, this approach requires to define the neighbourhood of each site. The combinatorial issue of model selection also prevents the use of such approaches when data are non-spatial.

The general approach to model multivariate distributions is to use copulas. Indeed, copulas [START_REF] Nelsen | An Introduction to Copulas[END_REF][START_REF] Joe | Multivariate models and multivariate dependence concepts[END_REF] can be used to build a multivariate model by defining, on the one hand, the one-dimensional marginal distributions, and, on the other, the dependency structure. Among the copulas, the Gaussian and the Student ones are very popular since they model the pairwise dependencies. However, their likelihood has not a closed form when the variables are discrete. It can be approached by numerical integrations which is not doable for high-dimensional data. Moreover, they require a quadratic number of parameters which leads to the curse of dimensionality for high-dimensional data. Alternatively the Archimedean copulas are relevant to reduce the number of parameters since they use a single parameter to model the dependencies between all the variables. Thus, this parameter characterizes a general dependency over the whole variables but it also limits the interpretation. For instance, positive and negative dependencies cannot be modelled simultaneously. Moreover, the evaluation of the likelihood requires the evaluation of an exponential number of terms, so it is not doable for high-dimensional data. Finally, vine copulas [START_REF] Kurowicka | Dependence modeling: vine copula handbook[END_REF] are a powerful alternative since they allow the specification of a joint distribution on d variables with given margins by specifying d 2 bivariate copula and conditional copula. Note that the vine copulas generalize and increase the flexibility with respect to the dependencies trees.

The one factor copulas [START_REF] Knott | Latent variable models and factor analysis[END_REF] enable to reduce the number of parameters and thus to deal with high-dimensional data. This approach assumes that the dependencies between the observed variables are explained by a continuous latent variable. This approach can be used for modelling continuous data variables [START_REF] Krupskii | Structured factor copula models: theory, inference and computation[END_REF], extreme-value continuous data [START_REF] Mazo | A flexible and tractable class of one-factor copulas[END_REF] or ordinal data [START_REF] Nikoloulopoulos | Factor copula models for item response data[END_REF].

In this work, we introduce a new family of one factor distributions that can be written as a specific one factor copula. For modelling more complex dependency structures, we extend this family by allowing a partition of the set of observed variables into independent blocks, where each block follows the new one factor distribution. The resulting family respects the five features listed by [START_REF] Nikoloulopoulos | Copula-based models for multivariate discrete response data[END_REF]. According to this specific distribution, each variable is described by three parameters: a continuous parameter indicating its marginal probability, a continuous parameter indicating the strength of the dependency with the rest of variables of the block (through the latent variable) and a discrete parameter indicating if the dependency is positive or negative.

Since the proposed distribution is a specific copula for discrete data, parameter inference is achieved by a two step procedure named Inference Function for Margin (IFM, see [START_REF] Joe | Multivariate models and multivariate dependence concepts[END_REF][START_REF] Joe | Asymptotic efficiency of the two-stage estimation method for copula-based models[END_REF]). Model selection consists in finding the best partition of the variables into blocks according to the Bayesian Information Criterion (BIC; [START_REF] Schwarz | Estimating the dimension of a model[END_REF]; [START_REF] Neath | The Bayesian information criterion: background, derivation, and applications[END_REF]). Although this information criterion is defined with the maximum likelihood estimates, an extension has been proposed with the parameter estimates resulting from the IFM [START_REF] Gao | Composite likelihood bayesian information criteria for model selection in high-dimensional data[END_REF]. For high-dimensional data, an exhaustive approach computing the BIC for each possible model is not doable. Therefore, we propose a deterministic two step procedure for the model selection. First, a small subset of models is extracted from the whole competing models by a deterministic procedure based on a Hierarchical Ascendant Classification (HAC) of the variables by using their empirical Cramer's V. Second, the BIC is computed for the models belonging to this subset and the model maximizing this criteria is returned. We show that this approach is asymptotically consistent. Indeed, Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo statistical methods[END_REF] is used for detecting the model maximizing the BIC criterion. Alternatively, a Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo statistical methods[END_REF] can also be used for detecting the model maximizing the BIC criterion. However, we numerically show that the deterministic procedure obtains similar results, in a strongly reduced computing time, as the stochastic one. Therefore, we advise to use the deterministic procedure.

The paper is organised as follows. Section 2 introduces the new family of the specific one factor distributions per independent blocks. Section 3 presents the parameter inference with the IFM procedure. The model selection issue is detailed in Section 4. Section 5 numerically compares both model selection procedures. Section 6 illustrates the approach on a real data set. Section 7 concludes this work. All the mathematical proofs are in appendix. The R package MvBinary implements the proposed method and contains the real data set. It is available on CRAN and the url http://mvbinary.r-forge.r-project.org/ proposes a tutorial for reproducing the application described in Section 6.

Multivariate distribution of binary variables 2.1 Blocks of independent variables

The aim is to model the distribution of the d-variate binary vector X = (X 1 , . . . , X d ).

Variables are grouped into b independent blocks for modelling different kinds of dependencies. Thus, the vector ω = (ω 1 , . . . , ω d ) determines the block of each variables since ω j = b indicates that X j is assigned to block b with 1 ≤ b ≤ b. Therefore, independence between blocks implies

∀1 ≤ j ≤ j ≤ d : ω j = ω j =⇒ X j ⊥ X j . (1)
Vector ω defines a model which is unknown and which has to be infered from the data. Variables affiliated to block b are mutually dependent and are denoted by X {b} = (X j : ω j = b). Obviously, this approach allows to model dependencies between all the variables (i.e. b = 1 then ω j = 1 for all 1 ≤ j ≤ d) or independence between all the variables (i.e. b = d then ω j = j for all 1 ≤ j ≤ d). The probability mass function (pmf) of the realisation x = (x 1 , . . . , x d ) is

p(x|ω, θ) = b b=1 p(x {b} |θ b ), (2)
where θ = (θ b ; b = 1, . . . , b) groups the model parameters, where θ b groups the parameters of the variables of block b. Finally, p(.|θ b ) is the pmf of variables affiliated to block b and each block is assumed to follow the one-factor distribution described in the following.

One-factor distribution per blocks 2.2.1 Conditional block distribution

In block b, dependencies between variables are characterised through a random continuous variable U b which follows a uniform distribution on [0, 1]. More precisely, variables of block b are independent conditionally on U b . So, the pmf of variables affiliated to block b is

p(x {b} |u b , θ b ) = j∈Ω b p(x j |u b , θ j ), (3)
where θ b = (θ j ; j ∈ Ω b ), θ j denotes the parameters related to variable X j detailed below, and where Ω b = {j : ω j = b} is the set of the indices of the variables of block b. Therefore, the specific conditional distribution of x {b} is assumed to be a product of Bernoulli distributions whose parameters are defined according to the value of u b . Indeed,

for j ∈ Ω b p(x j |u b , θ j ) = p x j j (1 -p j ) 1-x j with p j = (1 -ε j )α j + ε j 1 δ j {u b <α j } 1 1-δ j {u b >1-α j } , (4)
where θ j = (α j , ε j , δ j ) groups the parameters related to variable X j where:

• the continuous parameter α j ∈ (0, 1) corresponds to the marginal probability that X j = 1 since one can easily verify that for j ∈ Ω b ,

1 0 p(X j = 1|u b , θ j )du b = α j ,
• the continuous parameter ε j ∈ (0, 1) reflects the dependency strength between X j and the other variables of block j since the stronger the ε j , the more correlated are the variables of the block (see Proposition 2.3),

• the binary parameter δ j ∈ {0, 1} indicates the nature of the dependency, since δ j = 1 if the observed variable is positively dependent with the latent variable and δ j = 0 otherwise. Thus, two variables X j and X j affiliated to the same block (i.e. ω j = ω j ) are positively correlated if δ j = δ j and they are negatively correlated if

δ j = 1 -δ j .
Note that the model identifiability is discussed is the next section.

The parametrization of ( 4) is convenient for the model interpretation. However, we introduce the following new parametrization which simplifies the likelihood computation. Conditionally on u ω j , x j follows a Bernoulli distribution whose the parameters are only determined by a relation between u ω j and real β j = α δ j j (1 -α j ) 1-δ j which corresponds to the marginal probability that X j = δ j . Indeed, for u ω j ∈ [0, β j ), the conditional distribution X j |u ω j , θ j is a Bernoulli distribution B(λ j ) where λ j = (1 -ε j )α j + ε j δ j . Moreover, for u ω j ∈ [β j , 1], the conditional distribution X j |u b , θ j is a Bernoulli B(ν j ) where ν j = (1 -ε j )α j + ε j (1 -δ j ). Thus, (4) can be summarized as follows

p(x j |u b , θ j ) = λ x j j (1 -λ j ) 1-x j if 0 ≤ u b < β j ν x j j (1 -ν j ) 1-x j if β j ≤ u b < 1 . (5)

Marginal block distribution

Obviously, the realizations u b are not observed, but the distribution of the observed variables X {b} results from the marginal distribution of the pair (X {b} , U b ). So, the pmf of x {b} is defined by

p(x {b} |θ b ) = 1 0 p(x {b} |u b , θ b )du b . (6)
We now describe the properties of the block distribution. All proofs are given in Appendix A. For respecting the feature (F3) of [START_REF] Nikoloulopoulos | Copula-based models for multivariate discrete response data[END_REF] and for dealing with high-dimensional data, the block distribution needs to have a closed form. This explicit pmf is detailed in the following proposition.

f b (j 0 ; θ b ) = j 0 j=1 ν x (b,j) (b,j) (1 -ν (b,j) ) 1-x (b,j) d b j=j 0 +1 λ x (b,j) (b,j) (1 -λ (b,j) ) 1-x (b,j) , (8)
where x (b,j) := x σ b (j) denotes the j-th variable (according to permutation σ b ) assigned to block b, λ (b,j) := λ σ b (j) , ν (b,j) := ν σ b (j) and where j 0 j=j 0 +1 is one. The strength of the proposed model is its easy interpretation. The parameter interpretation is allowed by the property of identifiability now presented.

Proposition 2.2 (Model identifiability) The distribution defined by ( 7) is identifiable under the following constraints:

• δ (b,1) = 1 if d b > 2, δ (b,1) = 1; • ε (b,1) = ε (b,2) if d b = 2; • δ (b,1) = 1 and ε (b,1) = 0 if d b = 1; where δ (b,j) := δ σ b (j) , ε (b,j) := ε σ b (j) .
The proposed model allows a wide range of dependencies. The following proposition is related the model parameters and Cramer's V. Thus, we can see that the full dependency (respectively anti-dependency) can be modelled by putting ε j = ε j , α j = α j and δ j = δ j (respectively ε j = ε j , α j = 1 -α j and δ j = 1 -δ j ).

Proposition 2.3 (Dependency measures) The dependency between two binary variables is often measured with Cramer's V. For the distribution defined by (7), Cramer's V between X j and X j is zero. Moreover, for j and j and

β (b,j) ≤ β (b,j ) V (X j , X j ) = ε (b,j) ε (b,j ) β (b,j) (1 -β (b,j ) ) β (b,j ) (1 -β (b,j) ) . ( 9
)
3 Parameter inference

Inference function for Margins

We observed a sample x = (x 1 , . . . , x n ) assumed to be composed of n independent realizations of the proposed model. The likelihood related to model ω is defined by

p(x|ω, θ) = n i=1 b b=1 p(x i{b} |θ b ). ( 10
)
The log-likelihood function is defined by

L(α, δ, ε; x, ω) = n i=1 b b=1 ln p(x i{b} |θ b ), (11)
where α = (α j ; j = 1, . . . , d), δ = (δ j ; j = 1, . . . , d) and ε = (ε j ; j = 1, . . . , d). The proposed distribution is a multivariate copula-based model since each multivariate parametric distribution can be defined as a copula. When the model at hand is a copula with discrete margins, the maximization of the likelihood is quite difficult. Therefore, we use the Inference Function for Margins (IFM) procedure [START_REF] Joe | Multivariate models and multivariate dependence concepts[END_REF]. This estimation procedure is based on two optimization steps. The first step maximizes the likelihood of univariate margins. The second step maximizes the dependency parameters with the univariate parameters hold fixed from the first step. [START_REF] Joe | Asymptotic efficiency of the two-stage estimation method for copula-based models[END_REF] shows the asymptotical efficiency of such a procedure. Thus, the parameters θ = ( α, δ, ε) are estimated by the two following steps: Margin step: for j ∈ {1, . . . , d}

αj = 1 n n i=1
x ij ,

Dependency step: ( δ, ε) = arg max (δ,ε) L( α, δ, ε; x, ω).
The margin step is easily performed, but the search of ( δ, ε) at the dependency step requires solving equations having no analytical solution (except when d b = 2). This step is also achieved by using the latent structure of the data when d b > 2 (details are given in Section 3.2). When

d b = 2, for j ∈ Ω b : δ(b,2) = 1 if n 11 ≥ α(b,1) α(b,2) 0 if n 11 > α(b,1) α(b,2) and ε(b,1) = ε(b,2) = |n 11 -α(b,1) α(b,2) | β(b,1) (1 -β(b,2) ) , (12)
where

n 11 = 1 n n i=1 x ij 1 x ij 2 with j 1 ∈ Ω b , j 2 ∈ Ω b and j 1 = j 2 .

An EM algorithm for the dependency step

Since the blocks of the one-factor distributions imply latent variables, it is natural to perform the dependency steps of the IFM procedure with an Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF][START_REF] Mclachlan | The EM algorithm and extensions[END_REF]

when d b > 2.
The complete-data likelihood is defined by

L(θ; x, u, ω) = n j=1 L(θ j ; x, u, ω) (13)
where

L(α j , δ j , ε j ; x, u, ω) = n i=1 z ij (x ij ln λ j + (1 -x ij ) ln(1 -λ j )) (14) + (1 -z ij ) (x ij ln ν j + (1 -x ij ) ln(1 -ν j )) ,
where

z ij = 1 if 0 ≤ u iω j < β j and z ij = 0 if β j ≤ u iω j ≤ 1.
The EM algorithm is an iterative algorithm which alternates between two steps: the computation of conditional expectation of the complete-data log-likelihood (E step) and its maximization (M step) on (δ, ε). Note that the estimate α is not modified by the algorithm. Its iteration [r] is written as: E step: Computation of the complete-data log-likelihood, for j ∈ {1, . . . , d}

t ij (θ [r] ) = E[Z ij |x i , ω, θ [r] j ] = λ [r] j β [r] j λ [r] j β [r] j + ν [r] j (1 -β [r] j ) . ( 15
)
M step: Maximization over (δ j , ε j ), for j ∈ {1, . . . , d} δ

[r+1] j = 1 {maxε j L( αj ,δ j =1,ε j ;x,t [r] ,ω)>maxε j L( αj ,δ j ,ε j ;x,t [r] ,ω)} , (16)
ε [r+1] j = arg max ε j L(α j , δ [r+1] j , ε j ; x, t [r] , ω), (17)
where θ

[r]

j = (α j , δ [r] j , ε [r] j ), λ [r] j = (1 -ε [r] j )α j + ε [r]δ [r] j j , ν [r] j = (1 -ε [r] j )α j + ε [r]1-δ [r] j j
. Thus, the M step involves the search of the maximum over ε j ∈]0, 1[ of L(α j , δ j , ε j ; x, t [r] , ω). This maximization is easily performed since it only leads to solve a quadratic equation as shown by Appendix B.

Model selection

Information criterion

Model selection is obviously necessary when we are faced with model-based statistical inference. When the model pmf is given by ( 2), selecting a model means identifying the repartition of the variables into independent blocks. The challenge also consists of finding the best model according to the data among a set of competing models. In a Bayesian framework, the best model is defined by the model having the highest posterior probability. By assuming that uniformity holds for the prior distribution of ω, the best model also maximizes the integrated likelihood p(x|ω) where

p(ω|x) ∝ p(x|ω) with p(x|ω) = Θ p(x|ω, θ)p(θ|ω)dθ, (18)
and p(θ|ω) corresponds to the prior distribution of the parameters of model ω. However, this integral has not a closed form. In thus case, the BIC [START_REF] Schwarz | Estimating the dimension of a model[END_REF] is used for approaching the logarithm of the integrated likelihood by using a Laplace approximation.

It is defined by

BIC(ω) = L( θω ; x, ω) - ν ω 2 ln(n), (19)
where ν ω corresponds to the number of continuous parameters involved in model ω and where θω is the MLE of model ω. As shown by [START_REF] Gao | Composite likelihood bayesian information criteria for model selection in high-dimensional data[END_REF], the MLE can be replaced in ( 19) by the estimate provided by the IFM procedure. Thus, we want to obtain model ω which maximizes the BIC criterion among all the competing models.

The number of competing models is too huge for applying an exhaustive approach. Therefore, Section 4.2 presents a deterministic procedure for model selection. This procedure applies a filter among the competing models and only selects d models. Then, the BIC criterion is computed for each of the selected models. We show that this procedure returns the correct model ω asymptotically with probability one. Moreover, Section 4.3 presents a stochastic algorithm which finds ω . Section 5 shows that both procedures have the same behaviour for detecting the true model, but that the deterministic procedure is drastically faster than the stochastic procedure. Both procedures are implemented in the R package MvBinary, but we advise to only use the deterministic procedure for computing reasons.

Deterministic approach for model selection

This deterministic procedure has two steps. First, the reduction step reduces the number of competing models to only d competing models. Second, the comparison step computes the BIC criterion for each of the d resulting models and the model maximizing the BIC criterion is returned.

The reduction step decreases the number of competing models by using the empirical dependencies between the variables. Indeed, it performs the Hierarchical Ascendant Classification (HAC) of the variables based on the empirical Cramer's V. This procedure proposes d partitions corresponding to the d competing models on which the BIC criterion will be computed. Each model proposed by the HAC is relevant since it models the strongest empirical dependencies. Moreover, the HAC provides embedded partitions of variables and then reduces the calls to the EM algorithm.

The deterministic procedure based on HAC performs the model selection with the two following steps: Reduction step performs the HAC based on the empiric Cramer's V to defined the d partitions of the variables. Comparison step computes BIC ω [k] for k = 1, . . . , d, where ω [k] is such that each block b is composed by the variables affiliated to class b by the partition into k classes of the HAC. The procedure returns arg max k=1,...,d BIC ω [k] . Proposition 4.1 (Consistency of the HAC-based procedure) The HAC-based procedure is asymptotically consistent ( i.e. it returns the true model with probability one when n grows to infinity).

Proof is given in Appendix B.

Stochastic approach for model selection

Model ω can be determined through a Metropolis-Hastings algorithm [START_REF] Robert | Monte Carlo statistical methods[END_REF]. This algorithm performs a random walk over the competing models and its unique invariant distribution is proportional to exp BIC(ω) . Therefore, ω is the mode of its stationary distribution. It is also sampled with probability one by the algorithm when the number of iterations R grows to infinity.

At iteration [r], the algorithm samples a model candidate ω from the distribution q(.|ω [r] ) where ω [r] corresponds to the current model. More precisely, candidate ω is equal to the current model ω [r] except for variable j [r] randomly sampled which is affiliated into block b [r] randomly sampled in {1, . . . , max(ω [r] ) + 1}. This candidate is accepted with a probability equal to

ρ [r] = exp BIC( ω) q(ω [r] | ω) exp BIC(ω [r] ) q( ω|ω [r] ) . ( 20
)
This algorithm performs R iterations and returns the model maximizing the BIC criterion. In practice, there may be almost absorbing states, so different initialisations of this algorithm ensure to visit ω . Thus, starting from ω [0] , uniformly sampled among the competing models, the algorithm performs R iterations and returns arg max r=1,...,R BIC(ω [r] ).

Its iteration [r] performs the two following steps: Candidate step: ω is sampled from q(.|ω [r] ).

Acceptance/reject step: defined ω [r] with r-1] otherwise .

ω [r] = ω with probability ρ [r] ω [
5 Numerical experiments

Suitability of the HAC-based procedure

This simulation shows the relevance of competing models provided by the reduction step of the HAC-based procedure. Data are simulated from the proposed model with the following parameters d = 10, δ j = 1, α j = 0.4 and

ω j = 1 if 1 ≤ j ≤ 5 2 if 6 ≤ j ≤ 10 . ( 21
)
For different sizes of sample n and strengths of dependencies ε j , we check if the true model belongs to the list of models returned by reduction step of the HAC-based procedure.

Table 1 shows the results obtained on 50 samples for different values of (n, ε j ). Thus, whatever the strength of dependencies, the procedure asymptotically proposes the true model. Obviously, for a fixed sample size, results are better when the dependencies are strong since the number of times where the true model belongs to the list of models is increasing with the dependency strength.

n|ε

Comparison of model selection procedures

Both procedures of model selection are compared on data sampled from the proposed model with the parameters defined in (21). To compare the quality of the estimates returned by both procedures, we use the Kullback-Leibler divergence. As shown by Table 2, both procedures are consistent since the Kullback-Leibler divergence asymptotically vanishes. Moreover, the estimates have the same quality (equal value of the Kullback-Leibler divergence). However, the HAC-based procedure is considerably faster than the Metropolis-Hastings procedure as shown by Table 3

Model selection for high-dimensional data

This section shows the behaviour of the HAC-based procedure in high dimension. Data are generated from a model with blocks of five dependent variables (d b = 1), with equal marginal probabilities (α j = 0.4) and equal dependency strength (ε j = 0.4 and δ j = 1).

For different sizes of sample and numbers of variables, 50 samples are generated. Table 4 shows the relevance of the deterministic procedure by using the Adjusted Rand Index (ARI) to compare the true partition of the variables into blocks and its estimated. Indeed, whatever the number of variables, the procedure provides the true model with probability one when n grows to infinity. However, for small samples the procedure can provide a model slightly different to the true model. n|d 10 20 50 100 200 50 0.11 ( 0.18 ) 0.07 ( 0.03 ) 0.10 ( 0.06 ) 0.10 ( 0.04 ) 0.06 ( 0.02 ) 100 0.35 ( 0.33 ) 0.35 ( 0.24 ) 0.24 ( 0.14 ) 0.22 ( 0.07 ) 0.15 ( 0.03 ) 200 0.85 ( 0.27 ) 0.78 ( 0.20 ) 0.67 ( 0.11 ) 0.56 ( 0.07 ) 0.43 ( 0.05 ) 400 0.97 ( 0.09 ) 0.95 ( 0.07 ) 0.95 ( 0.05 ) 0.91 ( 0.05 ) 0.86 ( 0.05 ) 800 1.00 ( 0.00 ) 1.00 ( 0.01 ) 1.00 ( 0.01 ) 1.00 ( 0.01 ) 1.00 ( 0.01 ) Table 4: Mean (in bold) and standard deviation (in parenthesis) of the ARI between ω 0 and ω .

6 Application to plant distribution in the USA Dataset Data has been extracted from the USA plants database, July 29, 2015. It describes 35583 plants by indicating if they occur in 69 states (USA, Canada, Puerto Rico, Virgin Islands, Greenland and St Pierre and Miquelon). By modelling the data distribution, the flora variety of each states could be characterized. Moreover, one can expect bring out geographic dependencies between the variables. The data are available in the R package MvBinary which implements the proposed method.

Experiment conditions

The model selection is achieved by the deterministic algorithm (see Section 4.2) where the Ward criterion is used for the HAC. The EM algorithm is randomly initialized 40 times and it is stopped when two successive iterations increase the log-likelihood less than 0.01.

Model coherence

Figure 1 shows the relevance of the dependencies detected by the estimated model. Indeed, Figure 1a shows the correspondence between Cramer's V computed with the model parameter and the empirical Cramer's V, for each pair of variables claimed to be dependent by the estimated model. Moreover, Figure 1b shows that the estimated model well represents the main dependencies.

The estimated model is composed of 10 blocks of dependent variables. Figure 2 shows that this block repartition has a geographic meaning.

Model interpretation

Parameters permit an easy interpretation of the whole distribution. The mean per block of the values of αj and εj are summarized by Figure 3. Note that the dependencies q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q 0.0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q FALSE TRUE Each block is composed of highly dependent variables (high values of parameters εj and δj = 1). Therefore, the knowledge of one variables of a block provides strong information about the other variables affiliated into this block. For instance, the most dependent block is Block 10 (composed by Prince Edward Island, Nova Scotia, New Brunswick, New Hampshire, Vermont, Maine, Québec and Ontario). Thus, a plant occurs in Ontario with probability αOntario = 0.14 while it occurs with a probability 0.83 if this plant occurs in Québec. The least dependent block is composed of tropical states (Virgin Islands, Puerto Rico and Hawaii). These weaker dependencies can be explained by large geographic distance. Finally, parameters α j allow to described the region by their amount of plants. Cold regions (Blocks 2, 3 and 10) obtains small values of αj while the "sun-belt" obtains large values of this parameter.

Conclusion

In this paper, we have introduced a new family of distributions for large binary datasets. This family implies that the variables are grouped into independent blocks and that each block follows a specific one factor distribution. This new family has many good properties. Indeed, it verifies the five features required by [START_REF] Nikoloulopoulos | Copula-based models for multivariate discrete response data[END_REF] for a "good" distribution. Moreover, it permits an easy interpretation of the whole distribution. The variable repartition puts the light on the main dependencies. Moreover, each variable is summarized by its marginal probability (parameter α j ) and by its strength (parameter ε j ) and its kind (parameter δ j ) of dependency with the other block variables. Finally, this model is suitable for modelling large binary data since its number of parameters is linear in d.

We have proposed to circumvent the combinatorial problem of model selection with a deterministic procedure which reduces the number of competing models by using the empirical dependencies. Although this procedure does not ensure the maximization of the BIC, its consistency has been demonstrated. Numerical experiments have shown that this approach provides estimates having the same quality as a stochastic (and optimal) procedure, but it strongly reduces the computing time. The R package MvBinary implements both procedures of inference and contains the data set used in the application.

Many extension of this work can be envisaged. Indeed, parsimony extensions could be introduced by imposing equality constraints between the block parameters (e.g ∀j ∈ Ω b , ε j = c b where c b ∈]0, 1]). Moreover, more complex dependencies could be modelled by considering more than one factor and by keeping the same kind of parametrization. However, the parameter estimation and the likelihood computation would be more complex. Indeed, the pmf of block b would be defined as a sum of (d b + 1) 2 terms, while it is currently a sum of d b + 1 terms.

Finally, this model could be an answer to difficult task of the binary data clustering with intra-component dependencies. Indeed, the clustering aim could be achieved by considering a finite mixture of the proposed distribution. However, the challenge of model selection would be a complex issue. Moreover, the model identifiability should be carefully studied.

A Proofs of the model properties Proof of Proposition 2.1 It suffices to remark that (6) can be decomposed into d b + 1 integrals whose bounds are given by the coefficients β (b,j) . By using the conditional independence between variables given in (3) and the conditional distribution of x j given by ( 5

b ) = β (b,1) 0 p(x {b} |u b , θ b )du + d b -1 j=1 β (b,j+1) β (b,j) p(x {b} |u b , θ b )du + 1 β (b,d b ) p(x {b} |u b , θ b )du = β (b,1) f b (0; θ b ) + d b -1 j=1 (β (b,j+1) -β (b,j) )f b (j; θ b ) + (1 -β (d) )f b (d b ; θ b ). Proof of Proposition 2.2 We define that the distribution is identifiable if for two vec- tors of parameters θ b = (α j , ε j , δ j ; j ∈ Ω b ) and θ b = (α j , ε j , δ j ; j ∈ Ω b ) such that ∀x {b} , p(x {b} |θ b ) = p(x {b} |θ b ) then θ b = θ b . (22)
Without loss of generality, we assume that α j ≤ α j+1 . The equality α j = α j is directly obtained since ∀j ∈ Ω b , α j = p(x j = 1|θ b ) = p(x j = 1|θ b ) = α j . The proof distinguishes three cases: one variable in the block (i.e. d b = 1) with the constraints δ (b,1) = 1 and ε (b,1) = 0; two variables in the block (i.e. d b = 2) with the constraints δ (b,1) = 1 and ε (b,1) = ε (b,2) ; more than two variables in the block (i.e. d b > 2) with the constraint δ (b,1) = 1. Proofs use the following probability: ∀

(j 1 , j 2 ) ∈ Ω b , p(x j 1 = 1, x j 2 = 1|θ b ) =    α j 1 α j 2 + ε j 1 ε j 2 α j 1 (1 -α j 2 ) if δ j 2 = 1 α j 1 α j 2 -ε j 1 ε j 2 α j 1 α j 2 if δ j 2 = 0 and α j 1 + α j 2 < 1 α j 1 α j 2 -ε j 1 ε j 2 (1 -α j 1 )(1 -α j 2 ) if δ j 2 = 0 and α j 1 + α j 2 ≥ 1 (23) If δ (b,j) = δ (
b,j) then without loss of generality we assume that δ (b,j) = 1 and δ (b,j) = 0. From (23), p(x 1 = 1, x j = 1|θ b ) > α (b,1) α (b,j) = α (b,1) α (b,j) > p(x 1 = 1, x j = 1|θ b ) but this is in contradiction to (22), so ∀j ∈ Ω b , δ (b,j) = δ (b,j) . Therefore, we have to prove the equality ε (b,j) = ε (b,j) . Case 1 (d b = 1 with δ (b,1) = 1 and ε (b,1) = 0). Then parametrization assumes that only parameter α is free. Equality α j = α j implies θ b = θ b . Case 2 (d b = 2 with δ (b,1) = 1 and ε (b,1) = ε (b,2) ). By using constraints ε (b,1) = ε (b,2) and ε (b,1) = ε (b,2) and by using ( 23), then ε 2 (b,1) = ε 2 (b,1) . Thus, θ b = θ b . Case 3 (d b > 2 with δ (b,1) = 1). ( 23) is verified by θ and θ . Moreover, we know that α j = α j and δ j = δ j , for j = 1, . . . , d. So, the following system S arises from ( 23

If ε (b,1) = ε (b,1) then ∃t = 1 such that ε (b,1) = tε (b,1) . Then, the first d b lines of (S) imply that ∀j = 2, . . . , d b , ε (b,j) = tε (b,j) . Since the last line of (S) implies that ε (b,2) ε (b,3) = ε (b,2) ε (b,3) /t 2 , positivity of ε (b,j) permits to conclude that ε (b,1) = ε (b,1) , so ∀j = 2, . . . , d b , ε (b,j) = ε (b,j) . Thus, θ b = θ b .

Proof of Proposition 2.3 We denote p hh = P (X σ b (j) = h, X σ b (j ) |ω, θ) with j < j and h ∈ {0, 1} and h ∈ {0, 1}. Then p 11 = α (b,j) α (b,j ) + r p 01 = (1 -α (b,j) )α (b,j ) -r p 10 = α (b,j) (1 -α (b,j ) ) -r p 00 = (1 -α (b,j) )(1 -α (b,j ) ) + r where r = ε (b,j) ε (b,j ) β (b,j) (1 -β (b,j ) ). Thus, ( 9) is obtained by applying the definition of Cramer's V.

B Details about the M-step of the EM algorithm

By using the definition of αj , αj = n 10 + n 01 where n 10 = 1 n n i=1 x ij (1 -t ij ) δ j (t ij ) 1-δ j and n 11 = 1 n n i=1 x ij (t ij ) δ j (1 -t ij ) δ j . Moreover, the expectation of the complete-data likelihood related to variable j is written as L(α j , δ j , ε j ; x, t, ω) = n 10 ln((1 -ε j )(n 11 + n 10 )) + n 11 ln((1 -ε j )(n 11 + n 10 ) + ε j ) (25)

+ n 00 ln(1 -(1 -ε j )(n 11 + n 10 )) + n 01 ln(1 -(1 -ε j )(n 11 + n 10 ) -ε j ),

where n 00 = 1 n n i=1 (1 -x ij )(1 -t ij ) δ j (t ij ) 1-δ j and n 01 = 1 n n i=1 (1 -x ij )(t ij ) δ j (1 -t ij ) δ j . For a fixed value of δ j , the argmax over ε j of L(α j , δ j , ε j ; x, t, ω) is denoted by ε j|δ j . The estimation of ε j|δ j is obtained by setting to zero the derivative of L(α j , δ j , ε j ; x, t, ω) over ε j . So, remarking that n 01 = 1 -n 11 -n 10 -n 00 ,

s 1 = -B - √ ∆ 2A
and

s 2 = -B + √ ∆ 2A , (28)
where ∆ = B 2 -4AC. By noting that ε j ∈]0, 1[, and that s 1 = (n 11 +n 10 )n 10 +(1-n 11 -n 10 )n 01 -2(n 11 +n 10 )(1-n 11 -n 10 ) < 0, we conclude that ε j|δ j = max(0, s 2 ).

Proposition 2. 1 (

 1 Explicit distribution) Let σ b be the permutation of Ω b such that for 1 ≤ j < j ≤ d b the following inequality holds β (b,j) ≤ β (b,j ) , where β (b,j) := β σ b (j) and where d b = card(Ω b ) is the number of variables assigned to block b. The integral defined by (6) has the following closed form p(x {b} |θ b ) = d b j=0 (β (b,j+1) -β (b,j) )f b (j; θ b ), (7) where we define β (b,0) = 0 and β (b,d b +1) = 1. Finally function f b (.) is defined by

  Empiric Cramer's V and Cramer's V estimated by the model.

  Boxplot of the Empiric Cramer's V for the modelled and not modelled dependencies

Figure 1 :

 1 Figure 1: Visualisation of the dependencies taken into account by the model.

Figure 2 :Figure 3 :

 23 Figure 2: Geographic coherence of the blocks of states (color indicates the block assignment)

  ), function p(x {b} |u b , θ b ) is a piecewise constant function of u b . Thus, for u b ∈ [β (b,j) , β (b,j+1) [, p(x {b} |u b , θ b ) is constant and equal to f b (j) defined by (8). Then, p(x {b} |u b , θ

  ) for (j 1 , j 2 ) = {(1, 2), (1, 3), . . . , (1, d b ), (2, 3,1) ε (b,2) = ε (b,1) ε (b,2) ε (b,1) ε (b,3) = ε (b,1) ε (b,3) ε (b,d b ) = ε (b,1) ε (b,d b ) ε (b,2) ε (b,3) = ε (b,2) ε (b,3)

n

  11 + n 00 -1 1 -ε j|δ j + n 11 (1 -n 11 -n 10 ) n 11 + n 10 + ε j|δ j (1 -n 11 -n 10 ) + n 01 (n 11 + n 10 ) (n 11 + n 10 )ε j|δ j + (1 -n 11 -n 10 ) = 0.(26) This equation is equivalent to the following quadratic equationε 2 A + εB + C = 0, (27)whereA = -(n 11 + n 10 )(1 -n 11 -n 10 ), B = n 11 (n 11 + n ) + n 00 (1 -n 11 -n 10 ) -(n 11 + n 10 ) 2 -(1 -n 11 -n 10 ) 2 and where C = n 11 (1 -n 11 -n 10 ) + n 00 (n 11 + n 10 ) + A. Let s 1 and s 2 be the two solutions of (27):

Table 1 :

 1 Number of times where the true model belongs to the list of models returned by the reduction step of the HAC-based procedure on 50 samples.

	j 0.2 0.3 0.4 0.5 0.6
	50	0	1	4	28 37
	100	1	1	20 41 49
	200	0	9	40 48 20
	400	2	29 47 50 50
	800	3	44 50 50 50
	1600 21 49 50 50 50
	3200 40 50 50 50 50

Table 2 :

 2 . So, we recommend to use the HAC-based procedure to perform the model selection in high dimension. Kullback-Leibler divergence obtained with the estimates provided by both procedure of model selection.

	n|ε j	0.2		0.3			0.4	0.5
		HAC MH HAC MH HAC MH HAC MH
	50	0.15 0.16 0.24 0.27 0.36 0.36 0.39 0.46
	100 0.08 0.09 0.13 0.14 0.22 0.21 0.12 0.17
	200 0.04 0.04 0.10 0.10 0.08 0.11 0.05 0.06
	400 0.03 0.03 0.06 0.06 0.03 0.03 0.03 0.03
	n|ε j	0.2		0.3			0.4	0.5
		HAC MH HAC MH HAC MH HAC MH
	50	11	217	10	250	9	278	8	381
	100	12	241	11	250	10	354	8	633
	200	14	276	13	308	11	662	9	912
	400	16	296	15	509	12	1218	9	933

Table 3 :

 3 Computing time in seconds required by the two procedures of model selection.

Consistency of the HAC-based procedure

The proof of Proposition 4.1 is done in three steps. First, we show that the HAC-based procedure applied to the theoretical Cramer's matrix is consistent. Second, we show that this result holds in a neighbourhood of the theoretical Cramer's matrix. Third, we conclude by using the convergence in probability of the empiric Cramer's matrix to the theoretical one.

Let M 0 ∈ [0, 1] d×d be the dissimilarity matrix based on Cramer's V computed with the true distribution defined by model ω 0 and its parameters θ 0 . So, for 1

with V 0 (X j , X j ) is the theoretical Cramer's V between X j and X j defined by

Since the true model ω 0 involves independence between blocks of variables, for 1 ≤ j, j ≤ d with ω 0 j = ω 0 j , M 0 (j, j ) = 1. We denote by µ 0 the greatest value of M 0 when the variables belong to the same block for the true model ω 0

Note that µ 0 < 1 since the variables affiliated into the same block are dependent. Finally,

b ; b = 1, . . . , d) denotes the partition provided by the HAC at its iteration [r], where Ω

[r] b is the set of the indices of the variables affiliated to block b at iteration [r]. We consider that the HAC is used with a classical dissimilarity measure D(., .) (min, max, mean or Ward).

Proof At iteration [0], each variable is affiliated into its own block, so

2 ) = arg min (j 1 ,j 2 ) M 0 (j 1 , j 2 ), then

The Ω [1] verifies (32).

At iteration [r], by definition ∀b, ∀(j, j

b 2 ). There are two cases to be considered, for all j 1 ∈ Ω

Since µ 0 < 1, (32) is verified.

Corollary B.2 (Consistency with theoretical matrix)

The HAC based on the dissimilarity matrix M provides the true model at its iteration db 0 where b 0 is the number of blocks defined by ω 0 .

Proof It is the only partition of b [0] classes which respects Proposition B.1.

Corollary B.3 (Consistency in a neighbourhood of the theoretical matrix)

The HAC based on dissimilarity matrixM belonging to a neighbourhood of M 0 , denoted by N (M 0 ), provides the true model at its iteration db 0 where

Proof Proof is based on the same reasoning as the proof of Proposition B.1, since we have

Proof of Proposition 4.1 The Law of Large numbers implies that the observed probability of each couple (j, j ) converges in probability to its true value: for any h ∈ {0, 1} and h ∈ {0, 1} phh

where

where ph• = ph0 + ph1 and p•h = p0h + p1h . Thus, the Mapping theorem (see for instance Theorem 2.7 page 21 of Billingsley ( 2013)) implies that V converges in probability to V 0 . So,

Thus, by applying Corollary B.3, the probability that ω 0 belongs to the model subset provided by the HAC procedure is equal to one. The consistency of the BIC criterion concludes the proof.