Rémi Choquet
email: remi.choquet@cefe.cnrs.fr

Olivier Gimenez

Towards built-in capture-recapture mixed models in program E-SURGE

Keywords: Adaptive method, Cluster of individuals, Gauss-Hermite quadrature, GLMM, Hidden Markov chain, Normal random effects, State uncertainty

Towards built-in capture-recapture mixed models in program E-SURGE R. Choquet, O.

Introduction

With the improvement of Markov Chain Monte-Carlo (MCMC algorithms), the development of new software (e.g. Winbugs, Lunn et al. (2000)) and the availability of more and more powerful computers, mixed models are becoming more and more popular in ecological science [START_REF] Clark | Why environmental scientists are becoming bayesians[END_REF], and particularly in capture-recapture (CR) [START_REF] Wintrebert | Joint modelling of breeding and survival in the kittiwake using frailty models[END_REF][START_REF] Link | Modeling association among demographic parameters in analysis of open population capture-recapture data[END_REF][START_REF] Gimenez | Winbugs for population ecologists: Bayesian modeling using markov chain monte carlo methods[END_REF]. A non exhaustive list of advantages is:

-The possibility to explicitly decompose the variance into several components.

-The possibility to model dependence among parameters or random effects.

-That boundary estimates are avoided for small sample size data sets.

In the CR research area, mixed models remain under used despite their potential. One reason is that practitioners lack sufficient statistical training to code (or pseudocode in the case of WinBUGS) Bayesian implementations of mixed effects models. In non-Bayesian applications, [START_REF] Burnham | Model selection and Inference: A practical Information-theoretic approach[END_REF] and [START_REF] Royle | Random effects and shrinkage estimation in capturerecapture models[END_REF] considered the use of shrinkage estimator for time random effects. However, the mandatory assumption of non-boundary time fixed effect estimates are often violated, making this approach difficult to implement for a general class of model. Moreover, it is not applicable to the scale of individual effect as individual fixed effects estimates are not available. Another reason is with computation challenges in non-Bayesian applications such as numerical integration. This paper addresses the latter situation.

In this paper, we propose a first step toward implementing mixed CR models in E-SURGE. We will only consider a simple case with independent and identically distributed random effects. We use the independence property to reduce the dimension of the integral associated with the marginal likelihood. Thus because fitting problems become more and more complex, in the ideal case, we would like to consider an adaptive approach where the selected algorithm is suited to the model of interest. Developing adapted algorithms has proven in the past to be efficient in M-SURGE [START_REF] Choquet | M-SURGE: New software specifically designed for multistate recapture models[END_REF] to deal with models with complex age structures. The situation here is even more challenging and the road will be quite long before achieving this goal.

In line with [START_REF] Coull | The use of mixed logit models to reflect heterogeneity in capture-recapture studies[END_REF], [START_REF] Mcclintock | Estimating abundance using mark-resight when sampling is with replacement or the number of marked individuals is unknown[END_REF], [START_REF] Gimenez | Incorporating individual heterogeneity in capturerecapture models using numerical integration[END_REF], we consider a numerical approximation using Gauss-Hermite quadrature to calculate the marginal likelihood which has no analytical solution in presence of random effects. Gaussian quadrature is known to work efficiently for low dimension integrals over a large class of problems [START_REF] Lemenuel-Diot | Estimating heterogeneity in random effects models for longitudinal data[END_REF][START_REF] Heiss | Likelihood approximation by numerical integration on sparse grids[END_REF]Arndt et al., 2006). Numerical integration using Gauss-Hermite quadrature was used in [START_REF] Gimenez | Incorporating individual heterogeneity in capturerecapture models using numerical integration[END_REF] to fit the Cormack-Jolly-Seber (CJS) model with individual random effects. We implement this approach for the set of multievent model [START_REF] Pradel | Multievent: An extension of multistate capture-recapture models to uncertain states[END_REF], extensions of multistate CR models to handle state uncertainty, which are already implemented in E-SURGE [START_REF] Choquet | Program E-SURGE: A software application for fitting multievent models[END_REF]. In addition, we have considered an adaptive Gauss-Hermite quadrature to get efficient algorithms for group random effects. To get built-in random effects, we extend also the language in GEMACO [START_REF] Choquet | Automatic generation of multistate capture-recapture models[END_REF] to automatically build matrices of constraints: Mixed models can be specified with simple phrases.

We show how E-SURGE makes it is easy to fit mixed models using two applications. The first application illustrates a model accounting for a potential effect of the observation process with individual random effects applied to an European Dipper data set. The second application illustrates a model with group random effects applied to a Bank Vole data set.

Notations

K is the number of occasions, NI is the number of individuals, G is the number of groups, h is the data set of capture histories, I s is the identity matrix of size s, N (µ, σ ) is the normal distribution with mean µ and variance σ , Θ is the vector of biological parameters, Θ i the vector of biological parameters for individual i inside which the different biological parameters (survival, transition, capture, ...) are set in a row, see [START_REF] Choquet | Program E-SURGE: A software application for fitting multievent models[END_REF] for details.

Models under consideration

We consider the background of multievent models [START_REF] Pradel | Multievent: An extension of multistate capture-recapture models to uncertain states[END_REF] to handle state uncertainty in analyzing CR data, which are already implemented in E-SURGE. These models belong to the class of Hidden Markov Models (HMM) as the number of states is assume finite (e.g. [START_REF] Cappe | Inference in Hidden Markov models[END_REF]) and consequently to the class of State-Space Models (SSM) [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF]. Multievent models are a restriction of SSM in the sense that the number of states and observations are finite which is often the case for CR data.

Generalized linear mixed model

The class of mixed effects models that E-SURGE may consider can be expressed in the form of generalized linear mixed models (GLMM). We consider for both fixed and random effects, two general sets of effect:

-Set of effect 1: time, age, cohort and group effects.

-Set of effect 2: individual effect.

Considering the general form of GLMM f (Θ ) = Xβ + Zb with β the vector of fixed effects and b the vector of random effects is computationally demanding because of the dimension of the problem with so many potential effects. Thus, we have implemented the following restricted form of GLMM by constraining separately categories 1 and 2.

f (Θ i ) = X 0 β 0 + X i β 1 + L ∑ l=1 Z l b l + L+P ∑ l=L+1 Z l,i b l,i i = 1, . . . , NI (1)
where b l ∈ R s l and b l,i ∈ R are random effects given by

b l ∼ N (0, σ 2 l × I s l ), l = 1, . . . , L, b l,i ∼ N (0, σ 2 l ), l = L + 1, . . . , L + P.
(2) b l , l = 1, . . . , L are random effects associated to the set of effects 1, s l is the number of levels of the random effect l (s l = G for a group random effects), b l,i , l = L + 1, . . . , L + P are individual random effects assuming that individuals are independent. Matrices X i are individual-specific matrices of individual covariates. They are never stored in the computer because of the memory size needed but rather they are computed each time (see appendix A). In the same way, matrices Z l,i contain either 0, 1 or values of individual covariates and are never stored. Because we assume that individual are independent then matrices of covariance for each random effects are diagonal. We use this property to implement efficient algorithms.

General expression of the marginal likelihood

Assuming that individuals are independent, the likelihood for fixed effects for the entire set of capture histories is obtained as the product of the probability P(h i |β ) of the likelihood for each history h i

∏ h i ∈h Pr(h i |β )
, see [START_REF] Choquet | Program E-SURGE: A software application for fitting multievent models[END_REF], p 849 for details on multievent models.

For a GLMM like equation ( 1), we get the marginal likelihood:

L(β , σ ) = ∏ h i ∈h Pr(h i |β )ω(σ , x)dx. ( 3
)
where σ is the vector of variances associated to the random effects and ω(σ , x) is a product of normal densities associated to random effects b l or b l,i .

One might think that the dimension of the integral in equation ( 3) is of the size of the random effects and in most cases the evaluation of L(β , σ ) is numerically intractable. However, we illustrate below two situations for which the computation of equation ( 3) in E-SURGE can be made by reducing the dimension of the integral.

3.3 Independent and identically distributed random effect for individual E-SURGE can now handle mixed models with individual random effects only (L = 0):

f (Θ i ) = X 0 β 0 + X i β 1 + P ∑ l=1 Z l,i b l,i (4)
with b l,i in the form of equation (2). There is no limit for the number of random effects that we can build. However for P > 2 the fitting step may be time consuming.

In this case and similar to [START_REF] Lemenuel-Diot | Estimating heterogeneity in random effects models for longitudinal data[END_REF]; [START_REF] Gimenez | Incorporating individual heterogeneity in capturerecapture models using numerical integration[END_REF], the marginal likelihood (3) is the product of the probabilities of all individual encounter histories. So equation (3) can be rewritten:

L(β , σ ) = ∏ h∈h R P Pr(h|β , x)ω(σ , x)dx. (5)
Example 1: Survival varying with individual covariates and random effect. The following model has been used in [START_REF] Gimenez | Nonparametric estimation of natural selection on a quantitative trait using mark-recapture data[END_REF] with a constant survival across time but dependent on an individual covariates (m) and from an individual random effect (b i ).

logit(φ ik ) = β 0 + β 1 m i + b i , i = 1, . . . , I (6)
where b i ∼ N (0, σ 2 b ), (i.i.d). In equation ( 6) P = 1, the computation of the marginal likelihood according to equation (5) leads to the evaluation of several integrals with a single dimension instead of a potentially large dimension according to equation (3).

3.4 Independent and identically distributed random effect for group E-SURGE can also handle mixed models with group random effects only (P = 0) like:

f (Θ i ) = X 0 β 0 + X i β 1 + L ∑ l=1 Z l b l (7)
with b l in the form of equation (2). There is no limit for the number of random effects that we can build. However for L > 2 the fitting step may be time consuming."

In this case and similar to section 3.3, the marginal likelihood (3) is the product of the marginal probability for each group. We denote h g , the set of encounter histories inside a group g considered together. So equation (3) can be rewritten:

L(β , σ ) = G ∏ g=1 R L ∏ h∈h g Pr(h|β , x) ω(σ , x)dx. (8)
Example 2: We consider a basic model where recapture rates vary with a group random effect.

logit(p g ) = β 0 + b g , g = 1, . . . , G (9)
where b g ∼ N (0, σ 2 b ), (i.i.d). In this situation L = 1 so for evaluating the marginal likelihood (8) we have only G one dimensional integrals to evaluate.

Numerical integration

Integrals in (5) and (8) have no analytical expression so require numerical approximations. After reparameterization, those integrals have the form:

K(y) = R d k(y, x) exp -x t x dx 1 . . . dx d ( 10
)
where the dimension d depends on the number of random effects. For equation ( 5), d = P, see also [START_REF] Gimenez | Incorporating individual heterogeneity in capturerecapture models using numerical integration[END_REF] for the reparameterization in the case of the CJS model with an individual random effect. We will not consider in this paper model with dependence among individuals or groups, some of them are formally described in [START_REF] Cam | Contribution of capture-mark-recapture modeling to studies of evolution by natural selection[END_REF]. Dependences can be accomodated by normal random effects with non diagonal covariance matrices. The marginal likelihood associated to such models can also be expressed as equation ( 10) after reparameterization, see for example [START_REF] Heiss | Likelihood approximation by numerical integration on sparse grids[END_REF]. In the particular case of dependence of individuals and groups, the dimension d of the integral will depend of the structure of the dependence.

Order and accuracy

Let X N be a set of N quadrature nodes of R d , X N = {x (1) , . . . , x (N) } with x (n) ∈ R d for n = 1, . . . , N. Let Ω N be the set of quadrature weights associated with the nodes X N , Ω N = {ω 1 , . . . , ω N } with ω n ∈ R for n = 1, . . . , N. The integrals K(y) can be approximated by II N defined by

II N = N ∑ n=1 ω n × k(y, x n ) ( 11
)
II N is called a quadrature formula and K(y) -II N the error functional corresponding to K(y). For relevant sets X N and Ω N , the error functional decreases as N increases. Roughly, two methods of integration exist with different choices for Ω N and X N . The first method is based on the random generation of X, see [START_REF] Rosenberg | Bernstein polynomials and monte-carlo integration[END_REF]; [START_REF] Genz | A stochastic algorithm for high-dimensional integrals over unbounded regions with gaussian weight[END_REF] with or without equidistant properties [START_REF] Zaremba | Mathematical basis of monte carlo and quasi-monte carlo methods[END_REF]; [START_REF] Gonzalez | Numerical integration in logistic-normal models[END_REF]. Under weak conditions, the simulated value is unbiased and √ N-consistent by a law of large numbers and independent of the problem. The major drawback of this method relies on the difficulty in determining an appropriate value for N in order to control the error of integration K(y)-II N . The second method is based on the use of quadrature, see Arndt et al. (2006). The Laplace [START_REF] Liu | A note on gauss-hermite quadrature[END_REF] method is not considered precise enough for complex CR likelihoods. Consequently, higher orders of Gauss-Hermite quadrature will be considered to obtain better precision. In this case, N increases with d and the order of the quadrature r thus in some cases we can replace the index N by the two indices d and r. For d = 1, we directly use the Gauss-Hermite quadrature because it adapts easily to the normal density. The nodes of the quadrature X are related to the zeros of the r-th Hermite polynomial with N = 2 × r -1. The Gauss-Hermite quadrature with r nodes is exact for polynomials of order N. Here as stated in [START_REF] Gimenez | Incorporating individual heterogeneity in capturerecapture models using numerical integration[END_REF], k(y, x) is not a polynomial but rather a composite function of a polynomial following an inverse logit function, which can be approximated as closely as long as N is large enough.

Moreover under mild condition (if |K(y) -II d,r+1 | < κ|K(y) -II d,r |, κ < 1), we can obtain the error estimate associated with a polynomial of order r by considering two successive orders as:

|K(y) -II d,r | < 1 1 -κ |II d,r+1 -II d,r |.
One severe limitation of Gauss-Hermite quadrature is the effort growing with the dimension of the integral. Multivariate integration formulas are often constructed as a tensor product of quadrature formulas. Let II 1,r be a sequence of quadrature rules on R then

II d,r = ∑ x 1 ∈X . . . ∑ x d ∈X k(y, x 1 , . . . , x d ) d ∏ dd=1 ω dd (x dd )
For dimension d and order r, a regular grid involves (M = 2 × r -1) d nodes. The calculation of equation ( 10) very soon becomes intractable. Thus, we try to reduce the order as much as possible with an appropriate scaling.

Centering and Scaling

If the integrand in a numerical integration is not well behaved, algorithms can perform poorly. For one dimension, the order of integration can be set higher to overcome this problem. But it is sometimes still not enough, and the error of integration can decrease very slowly. The problem in multiple dimensions is even worse because the number of nodes increases exponentially. Thus for equation ( 8), we use centering and scaling in a similar way as [START_REF] Lemenuel-Diot | Estimating heterogeneity in random effects models for longitudinal data[END_REF] to improve the integration.

In order to base the approximation of the integral (10) on an appropriate range of values of x, two transformations are made: centering and scaling. If we assume that k(y, x) is close to a normal distribution of mean x * and variance S, k(y, x) × exp ((xx * ) t S(xx * )) will be close to a constant and the approximation of the integral exact. Thus, we take x * as the mode of k(y, x) and S as half of the inverse of the hessian matrix of the same function.

With the overall transformation x = U -1 × z + x * with S = U t U then equation ( 10) becomes

K(y) = 1 detU R d k(y,U -1 × z + x * ) exp -(U -1 ×z+x * ) t (U -1 ×z+x * ) dz 1 . . . dz d = 1 detU R d k 1 (y, z) exp -z t ×z dz 1 . . . dz d (12) with k 1 (y, z) = k(y,U -1 × z + x * ) exp -(U -1 ×z+x * ) t (U -1 ×z+x * )+z t z .
We compute U as in [START_REF] Lemenuel-Diot | Estimating heterogeneity in random effects models for longitudinal data[END_REF]. To get the inverse of U and when necessary (i.e., when d > 1), we use an appropriate algorithm [START_REF] Xu | A divide-and-conquer method for the takagi factorization[END_REF] for the Takagy Factorization (a singular value decomposition of the symmetric matrix S).

Avoiding underflow and overflow

One particularity is that the approximation II N is not straightforward to implement in the case of equation ( 8). In fact, the function k(y, x) is a polynomial of high degree, with parameters belonging to [0, 1]. Thus very soon, values of the function k will be not representable in finite numerical precision, i.e. k(y, x) < 2 -128 . This limitation is called underflow. Thus we reformulate equation ( 10) to avoid underflow.

One solution to deal with underflow in equation ( 10) is to consider log k(y, x) so that the sum of individual contributions is involved rather than the product.

K(y) = R d k(y, x) exp -x t x dx 1 . . . dx d = R d exp(log k(y, x)) exp -x t x dx 1 . . . dx d (13)
However, overflow also occurs with exp(log k(y, x)). One way to deal with this new problem is to scale the quantity log k(y, x) by an appropriate value called dev * .

K(y) = exp(0.5 × dev * ) × R d exp(log k(y, x) + 0.5 × dev * ) exp -x t x dx 1 . . . dx d = exp(0.5 × dev * ) × R d k 2 (y, x) exp -x t x dx 1 . . . dx d ( 14
)
where k 2 (y, x) = exp(log k(y, x) + 0.5 × dev * ).

5 Implementation in program E-SURGE

Details of the global algorithm

To get the MLE, we consider a Quasi-Newton algorithm with a first order difference scheme to approximate the gradient and a relative tolerance of 10 -6 . The error made by using a first or second degree finite difference scheme does not influence the result in most of the case (Appendix B). In this paper, the results are obtained with a second-order finite difference scheme although the use of the first order scheme led to the same results in both applications.

For a group random effect, we consider the adaptative scaling described previously. To ensure precision of the quadrature scheme, we update the matrix K every 15 iterations of the quasi-Newton algorithm.

Description of the random effects

We extend the model definition language of the tool GEMACO used in M-SURGE and E-SURGE. We refer to [START_REF] Choquet | Automatic generation of multistate capture-recapture models[END_REF] for the description of the language used for fixed effects.

Although individual covariates are permitted for a fixed effect, individual fixed effects are not allowed in E-SURGE. Therefore, we introduce a new built-in keyword factor denoted IND for individual random effects, and implement random effects for groups with the keyword RANDOM, which translate fixed effects into random effects. These additions fit naturally into E-SURGE's model specification syntax. However contrary to traditional effect like TIME, AGE, GROUP, direct addressing of levels of IND (one level corresponding to one individual) is not currently allowed. We extend also the operator + to concatenate fixed effect and random effect to generate mixed model of the form (4,7). Examples include:

The phrase "I+XIND(1)+IND" models equation ( 6).

The phrase "I+RANDOM(GROUP)" models equation (9).

More generally, two general forms of phrase are currently allowed "phrase1+phrase2.IND" for equation (4) and "phrase1+RANDOM(phrase2)" for equation ( 7), where phrase1 and phrase2 are any general phrases for fixed effects.

6 Applications

Application 1

We consider the European Dipper data used in [START_REF] Lebreton | Modeling survival and testing biological hypotheses using marked animals: A unified approach with case studies[END_REF]. Although affecting only the two last occasions, we found a trap-dependence effect using program U-CARE [START_REF] Choquet | U-CARE: Utilities for performing goodness of fit tests and manipulating CApture-REcapture data[END_REF]). This local effect on time may be a consequence on the observation being done on the nest: an individual seen one year on a nest have a greater probability to be detected the next year on the same nest. This effect may be not detected the first years because of the flood years: their livelihood is closely dependant on stream as their nest are close to the water. Therefore, we considered a modified version of the standard Cormack-Jolly-Seber model incorporating a trap-dependence effect to account for this lack of fit. More precisely, we considered two detection probabilities at time t p * t and p t depending on whether an individual was previously captured or not at the previous occasion t-1.

As in [START_REF] Royle | Modeling individual effects in the cormack-jolly-seber model: A statespace formulation[END_REF], we tested for individual effects in both survival (φ ) and capture (p * , p). Similar to [START_REF] Lebreton | Modeling survival and testing biological hypotheses using marked animals: A unified approach with case studies[END_REF], the fixed effect for survival is the flood effect (t = 1, 4, 5, 6 vs t = 2, 3). With E-SURGE, the sentence for the survival TIME(1 4 5 6,2 3)+IND builds a model where

logit(φ t=(1,4,5,6),i ) = β 1 + b 1,i , (15)
logit(φ t=(2,3),i ) = β 2 + b 1,i .
where b 1,i ∼ N (0, σ 2 1 ). Using the formulation for trap-dependence described in [START_REF] Gimenez | Parameter redundancy in multistate capture-recapture models[END_REF], the set of states used in E-SURGE is {'Alive and captured previously', 'Alive and not captured previously', 'Dead'}. The sentence for capture TIME(1:4)+FROM.TIME(5 6)+IND builds a model where

logit(p * t=(1,...,4),i ) = β 3 + b 2,i , logit(p t=(1,...,4),i ) = β 3 + b 2,i , logit(p * t=(5 6),i ) = β 4 + b 2,i , (16)
logit(p t=(5 6),i ) = β 5 + b 2,i .
where b 2,i ∼ N (0, σ 2 2 ). We fit this model φ ( f lood + ind), p(partial(m) + ind) using r = 15 for the quadrature without centering and scaling. In this example, centering and scaling failed to improve the precision of the Gauss-Hermite quadrature because h(x) is not close to a normal distribution. Table (1) shows that no individual effect is detected. For any of the two models φ ( f lood + ind), p(partial(m) + ind) and φ ( f lood + ind), p(m + ind) for the full trap-dependence, the deviance does not even change compared to φ ( f lood), p(partial(m)) and φ ( f lood), p(m). In each case, σ 1 and σ 2 are estimated close to zero. Thus, the individual effect detected in [START_REF] Royle | Modeling individual effects in the cormack-jolly-seber model: A statespace formulation[END_REF] can be explained by unmodeled trap-dependence. This is supported also by simulation where data with trap-dependent effect were generated and then analyzed with an individual effect, the model with individual effect was far better than the constant model (unpublished results). A step-by-step procedure for the implementation of the model in E-SURGE is given in Appendix C.

Application 2

Hantaviruses are the ethiological agents of several more or less severe diseases in humans. In Northern, Western and Central Europe, Puumala hantavirus causes Nephropathia epidemica, a mild form of haemorrhagic fever with renal syndrome [START_REF] Klein | Emergence and persistence of Hantaviruses, volume 315 of Wildlife and Emerging Zoonotic Diseases: The Biology, Circumstances and Consequences of Cross-Species Transmission[END_REF]. However, for the virus, humans act only as dead-ends and its reservoir host is a small mammal, where b g ∼ N (0, σ 2 ) and the first mathematical parameter(FIRSTE) is fixed to 1 in the biological scale (p age=1 = 1). Note that this approach is equivalent to conditioning at first capture.

We fit the model using r = 15 with centering and scaling. A step by step procedure for the implementation of the model in E-SURGE is given in Appendix D. Results are represented in Table 2. The LRT test between models φ (i), p(i + r(g)) and φ (i), p(i) was accepted/rejected.

Discussion

We have implemented effective algorithms in E-SURGE to fit models and get maximum likelihood parameter estimates with random effects in cases where independence of individuals or groups is assumed. Instead of a MCMC approach, we consider a Gauss-Hermite scheme to cope with high-dimensional integrals. This approach has several advantages over Monte-Carlo approaches. For example, we obtained a cheap measure of the error. A better measure could be achieved for both methods by computing the discrepancy [START_REF] Frank | Computing discrepancies of smolyak quadrature rules[END_REF], but this option comes at a price. Furthermore, the language inside GEMACO has been extended to allow built-in mixed models incorporating individual and group random effects with i.i.d. assumptions. Thus, practitioners without experience in MCMC modeling can access these modeling tools within the scope of a familiar, userfriendly software package.

However several challenges remain:

-First, built-in time, age and cohort random effects still need to be implemented in E-SURGE. Dimensions of the associated integrals are very high and additional algorithms must be implemented to fit these models in a reasonable amount of time.

Much progress have been done the past few years with MCMC methods or EMtype algorithms [START_REF] Chaubert-Pereira | Combinaisons Markoviennes et semi-Markoviennes de modeles de regression[END_REF]. Our ultimate goal is to implement the most appropriate algorithm for each model. -Second, in the present version of E-SURGE, the independence of groups needs to be specified in a menu of E-SURGE. Description or automatic detection of independence between groups remains to be implemented. -Third, models with non-diagonal variance-covariance matrices as in [START_REF] Pinheiro | Unconstrained parametrizations for variancecovariance matrices[END_REF] have to be implemented. The description of dependence between individuals (which would allow fitting animal models like in Meyer ( 2001)) or groups is the first step of this implementation. -Fourth, besides fitting the models, model selection remains a challenging exercise with AIC a subject of controversy in the context of mixed models. Other measures like CAIC [START_REF] Vaida | Conditional akaike information for mixed-effects models[END_REF]) should be investigated. -Fifth, the acute computation of the rank of a mixed model for CR remains to be done to investigate if redundant parameters are present in the models. In the current version of E-SURGE, the rank is set equal to the rank of the hessian of the likelihood. -Last but not least, goodness-of-fit procedures are missing for CR mixed models. We hope that further progress will be made with GOF as an inappropriate random effect can be selected even with marginal overdispersion. This is the case in application 1.

A Implementation of the structure of X i and Z l,i X i A × B i where B i is set of square matrices dependent on individuals. This structure allows us to consider model like: I+T.XIND. In this case, A is a time dependent matrix and the matrix B i is a (K -1) × (K -1) diagonal matrix which entries are the value of the covariates for individual i.

Z l,i = C l × D l,i where D 1,i is a square matrix depending of the individual and of the effect. This structure allows us to consider model like: I+XIND+RAND(XIND). In this case L = 0, P = 1, C 1 is a constant vector of 1 and D 1,i is the value of the covariates for individual i.

B Influence of the finite difference scheme

We check that the tolerance (Tol = 10 -6 ) used for the quasi-Newton algorithm on the gradient and the error (Err) made by the approximation of the gradient using a finite difference scheme and a numerical integration are consistant, i.e. of the same order. Let ε be the computer precision, we will demonstrate that the global error Err defined by:

Err = ∂ K(y) ∂ y | y 0 - N ∑ n=1 ω n × k(y 0 + √ ε, x n ) -k(y 0 , x n ) √ ε
can be decomposed as a sum of two sources of error; the error of the finite difference scheme applied to the gradient of g and the error made by the quadrature formulae applied to the gradient of g. By the Fubini theorem,

∇K(y 0 ) = R d ∇ y k(y 0 , x) exp -x t x dx 1 . . . dx d
The error made by approximating a derivatives by a first order finite-difference scheme is O( √ ε) (see [START_REF] Dennis | Numerical Methods for Unconstrained Optimization and Nonlinear Equations[END_REF]) so

Err = ∇K(y 0 ) - N ∑ n=1 ω n × ∇ y k(y 0 , x n ) + O( √ ε) Err = R d ∇ y k(y 0 , x) exp -x t x dx 1 . . . dx d - N ∑ n=1 ω n × ∇ y k(y 0 , x n ) + O( √ ε)
For Matlab on a 32 bits windows system, ε = 2.2204 × 10 -16 so as soon as the error made by the quadrature formulae is lower than Tol, the global error is lower than Tol. For a lower tolerance (Tol = 10 -6 ) used for the quasi-Newton algorithm on the gradient then a second-order finite-difference scheme should be used.

C Application 1 with E-SURGE

E-SURGE can accept capture-recapture data either in MARK or BIOMECO format. The two types of data file are not very different: each row corresponds to a particular capture history followed by the number of individuals with that history. In MARK format this count is followed by a semi-colon and in BIOMECO format the data are preceded by the number of different capture histories and the number. In this study, there is only one site. We will fit here the model φ ( f lood)p(partial(m) + ind) described in part 6.1.

C.1 Starting E-SURGE

From E-SURGE, start a new session named 'result.mod'. Read in the data file and tell to E-SURGE that there is no individuals covariates. Check that the numbers of capture occasions, groups and events are correct (in this case 7 capture occasions, 1 group and 2 events): E-SURGE makes assumptions about the number of states, but these need to be modified depending on the problem you have to treat. In this case, we change the number of states to three. We will also have to set the number of age classes to one as for the present, we won't consider any age effect.

C.2 Fitting the model

Press the Modify button, and change the settings to specify there a single age class, a single group and 3 states (see Fig. 1).

Fig. 1 E-SURGE:

The number of states is set to 3 and the number of age classes is set to one.

The main menu should show the changes. Fitting the models in E-SURGE involves four steps:

1. The Gepat step: specifying which ones of all the potential parameters have to be estimated, which ones will be calculated as the complement to 1 (there is one such parameter per multinomial) and which ones correspond to impossible events or transitions and are fixed to zero; 2. The Gemaco step: specifying the effects (time, age) acting on the active parameters; 3. The IVFV step: specifying initial values for the optimization procedure and/or fixed values for the active parameters; 4. The RUN step: launching the optimization procedure.

C.2.1 Specifying the pattern matrices using the GEPAT interface

There are three types of parameters used in the definition of a multi-event model [START_REF] Pradel | Multievent: An extension of multistate capture-recapture models to uncertain states[END_REF]):

1. the initial state probabilities; 2. the transition probabilities; 3. the event probabilities.

Each type of parameters is gathered into a row-stochastic matrix, i.e. each row corresponds to a multinomial. (Each matrix can be further decomposed into a product of several stochastic matrices allowing for example to estimate separately survival and transition parameters. However, for the current model, only 1 step is required for each type of matrix.)

To enter the GEPAT interface, click the GEPAT button at the bottom of the main window. The GEPAT interface screen for the initial state pattern matrix appears.

By default, E-SURGE let all live states available as initial steps. The state dead, last in the list, is impossible and does not even appear. The last live step is taken by default to be the one whose probability will be calculated indirectly, as the complement of those of the other live states. This is specified in the above window using the following general conventions for Gepat:

a minus sign (-) indicates that a potential parameter is unavailable in the current model (impossible transition for instance). This is equivalent to fixing it to zero but more explicit; any greek letter (strike a latin letter and E-SURGE will show its greek counterpart) indicates a free parameter, one to be estimated directly. Note that the particular greek letter entered is totally irrelevant to E-SURGE. In particular, the same greek letter is used repeatedly by default within a pattern matrix ( by default for initial states); this does NOT mean that the parameters are being forced to be equal; the asterisk (*) indicates the parameter that is calculated indirectly, as the complement to 1 of all the other parameters on the same row. There MUST always be one and only one asterisk per row because each row corresponds to a multinomial. Note that the order of the states is chosen by the user except for the dead state that is always positioned last.

Here, the default pattern is not correct for transition and event. We need to change them.

For intial states, we implement the pattern given in Figure 2. For transitions, we need two steps, one for survival and one for capture. So we set the number of steps to 2. For survival, we implement the pattern given in Figure 3. For capture, we implement the pattern given in Figure 4. For event, we implement the pattern given in Figure 5. Press the "EXIT" button to return to the main window.

C.2.2 Specifying the model using the GEMACO interface

The GEMACO interface uses keywords to create a modelling sentence that indicates how parameters vary by time, over groups, over age classes, etc. At the end of the GEMACO procedure, a design matrix is created for each type of parameters. Each row of the design matrix will correspond to a parameter of the full model (all potential variability: time, age, group) and each column corresponds to a parameter of the actual model.

The GEMACO syntax is fairly intuitive but the "sentences" you enter in the GEMACO interface must respect some priority rules that we will not develop here. We encourage the user to read the E-SURGE user manual and the paper [START_REF] Choquet | Automatic generation of multistate capture-recapture models[END_REF] in which the GEMACO syntax is fully explained.

In this example, we only want to show how to use E-SURGE to fit our model.

For the trap-dependent model, the set of initial state probabilities and the set of event probabilites are empty. Click on the top "Initial State" button to go to the "Transition" screen.

The survival and the capture modelling

To specify the model on survival, we use the phrase T(1 4 5 6,2 3) for the flood effect, the new keyword IND for the individual random effect. As we combine all, the GEMACO sentence becomed T(1 4 5 6,2 3)+IND. Select the next step for transitions which corresponds to the modelling of trap dependence. We use the phrase T(1:4)+T(5 6).F+IND. At this stage, to create the design matrices, we click on the Gemaco item in the top menu and select the "call GEMACO (all phrases)" submenu. At this stage, all the model structures are specified and the design matrices appear in the left window of each screen of the GEMACO interface; press "EXIT" button to return to the main window.

C.2.3 Specifying the initial and fixed values using the IVFV interface

In E-SURGE, the user can choose the way to generate the initial values of the optimization procedure. They can be either "constant", "randomly generated" or "equal to the estimates of a previously fitted model". Once the type of initial values is chosen, the user can also fix the values for some parameters using the IVFV interface. Press the IVFV button to enter the interface (Fig. 6). As there is no need to specify neither fixed values nor initial values. Click on the "EXIT" button to return to the main window.

C.2.4 Running the model

Before running the model, we have to specified the method of integration, here we choose the classical Gauss-Hermite method (set by default) described in the paper with 29 quadrature nodes (r=15, set by default). We also tick the "compute C-I(Hessian)" option to get confidence intervals. The model is now ready to be fitted to the data. Press the RUN button. We observe in Figure 7 that at the end of the fit, the estimate error made by the GH scheme to get the likelihood is lower than 10 -6 as -15.64 ≤ -6.

In Figure 8, we get the estimates for the model φ ( f lood + ind)p(partial(m) + ind) in the mathematical scale.

D Application 2 with E-SURGE

We will fit here the model φ ( f lood)p(partial(m) + ind) described in part 6.2 with E-SURGE.

D.1 Starting E-SURGE

From E-SURGE, start a new session and load the data set. For this application, we only need to change the number of age classes to one as for the present, we won't consider any age effect. Fig. 7 E-SURGE: The Output during the RUN step shows that 20 iterations are for convergence and that the error due to numerical integration is small (lower than 10 -6).

D.2 Fitting the model

Press the Modify button, and change the settings to specify that there is a single age class, nine groups and 2 states (see Fig. 9).

The main menu should show the changes. Fitting the models in E-SURGE involves four steps:

D.2.1 Specifying the pattern matrices using the GEPAT interface

Note that the order of the states is chosen by the user except for the dead state always positioned last. Here, the default pattern is one of the CJS model. We don't need to change them.

D.2.2 Specifying model using the GEMACO interface

In this example, we only want to show how to use E-SURGE to fit our model.

The initial states modelling

For the CJS model, the set of initial state probabilities is empty. Click on the top "Initial State" button to go to the "Transition" screen.

The survival and the capture modelling

To specify the model on survival, we use the phrase I (see Fig. 10). Click on the top "Transition" button to go to the "Event" screen.

The event modelling

Because the model conditions on the first capture occasion of each individual, the only event to model at the time of the first encounter is the site of capture. It is only later on that the event 'not encountered' becomes possible. Thus, the event probabilities at the time of the first encounter must be treated separately. This is achieved through the use of the keywords "firste", which stands for 'first encounter', and "nexte" , which stands for 'next encounters', respectively.

The output of Gepat is shown in the 'transitions pattern' subwindow at the bottom left. For each state (in row), there is only one active event which is the capture on the relevant (where the stands), and the other possible event taken as the complement is 'not encountered' (first column). For instance, for state 11, the first row, the individual can be encountered on site 1 (second column, probability to be estimated) or 'not encountered' (first column, probability calculated as 1 -the other probability). The active parameters are thus just the capture probabilities. At the time of the first encounter, the capture is certain and the capture probabilities will all be 1. At this stage, we cannot specify a fixed value, but we can specify that we need just one parameter common to all states by keeping "firste" by itself. Later, capture probability will be constant. Thus, the complete sentence is "firste+nexte+random(nexte.g)".

At this stage, to create the design matrices, we click on the Gemaco item in the top menu and select the "call GEMACO (all phrases)" submenu (see Fig. 11). All the model structures are now specified and the design matrices for fixed effect appear in the left window of each screen of the GEMACO interface; press the "EXIT" button to return to the main window.

D.2.3 Specifying the initial and fixed values using the IVFV interface

In E-SURGE, the user can choose the way to generate the initial values of the optimization procedure. They can be either "constant", "randomly generated" or "equal to the estimates of a previously fitted model". Once the type of initial values is chosen, the user can also fix the values for some parameters using the IVFV interface. Press the IVFV button to enter the interface.

The initial states probabilities

In this case, there is no need to specify neither fixed values nor initial values. Click on the top "Initial State" button to arrive at the "Transition" screen.

The survival-transitions probabilities

There is no need to fix values for the transition probabilities so this screen can be left in its default state. Click on the top "Transition" button to arrive at the "Event" screen.

The event probabilities

We can see here the different capture rate appearing in the definition of the model. The series of number indicate successively:

the line in the event matrix (corresponding to the state), the column the event matrix (corresponding to the event), the capture occasion, the age class, the group, the step in the matrix decomposition of the event matrix(here 1). Thus, the first parameter corresponds to the capture rate at the first capture occasion for the first age class, i.e. time of first encounter (A=1). This is the only parameter with A=1 because we have gathered all the capture rates relative to the first encounter into a single parameter. This parameter needs to be fixed to 1. We do this by entering the value 1 as "Initial Value" and ticking the box nearby (see Fig. 12). The other Fig. 12 E-SURGE: The first mathematical parameter corresponding to 'firste' is fixed to 1 by ticking the box nearby. The third mathematical parameter is the starting value for the square root of the standard error, it must be strictly positive. parameter correspond to the following capture rates (A=2); there is no need to fix these parameters. After all the fixed values have been specified, press the EXIT button.

D.2.4 Running the model

Before running the model, we have to specify the method of integration. Click on the button 'Advanced Numerical Options ¿ Modify' in the main window of E-SURGE. Here we choose the Adaptative Gauss-Hermite method described in the paper (the fourth value is set to 1) with 29 quadrature nodes (the fifth value is set to 15) (see Fig. 13). We also tick the submenu 'Random Effects for Independent Group Only' in the menu 'Models' (Fig. 14). We also tick the 'compute C-I(Hessian)' option to get confidence intervals and an estimated of the model rank. The model is now ready to be fitted to the data. Press the RUN button. We observe in Figure 15 that at the end of the fit, the estimate error made by the GH scheme to get the likelihood is lower than 10 -6 as -6.45 ≤ -6.

In Figure 16, we get the estimates for the model φ (.)p(i + r(g)) in the mathematical scale. The Output during the RUN step shows that only 10 iterations are needed for convergence and that the error due to numerical integration is small (lower than 10 -6).

Fig. 2

 2 Fig. 2 E-SURGE: Pattern for the intial states vector.

Fig. 3

 3 Fig. 3 E-SURGE: Pattern for survival matrix.

Fig. 4

 4 Fig. 4 E-SURGE: Pattern for the capture matrix for the trap-dependent model.

Fig. 5

 5 Fig. 5 Pattern for the event matrix.

Fig. 6

 6 Fig. 6 E-SURGE: The initial values fixed values interface. The two last parameters are the parameter associated to the two random effects.

Fig. 8

 8 Fig. 8 E-SURGE: Estimates and standard errors for the model φ ( f lood + ind)p(partial(m) + ind).

Fig. 9

 9 Fig. 9 E-SURGE: The number of age classes is set to one.

Fig. 10

 10 Fig. 10 E-SURGE: The sentence 'i' builds the model φ (.).

Fig. 11 E

 11 Fig. 11 E-SURGE: The sentence 'firste+nexte+r(nexte.g)' builds the model φ (i + r(g)).

Fig. 13 E

 13 Fig. 13 E-SURGE: We choose the Adaptative method (the fourth value is set to 1) with 15 quadrature nodes (the fifth value is set to 15).

Fig. 14

 14 Fig.14 E-SURGE: To speed-up calculations and improve the precision of the integration, we set that the random group effect is i.i.d.

Fig. 15 E

 15 Fig.15 E-SURGE:The Output during the RUN step shows that only 10 iterations are needed for convergence and that the error due to numerical integration is small (lower than 10 -6).

Fig. 16 E

 16 Fig. 16 E-SURGE: Estimates and standard errors for the φ (.)p(i + r(g)).

Table 1

 1 European Dipper example revisited. N p is the number of mathematical parameters of the model.

	Model			Np	Deviance
	φ ( f lood + ind), p(partial(m) + ind) 7	650.910
	φ ( f lood), p(partial(m) + ind)	6	650.910
	φ ( f lood + ind), p(partial(m))	6	650.910
	φ ( f lood), p(partial(m))		5	650.910
	φ ( f lood + ind), p(m + ind)	6	656.225
	φ ( f lood), p(m)			4	656.225
	φ ( f lood), p			3	660.102
	Model	Np	Deviance
	φ (i), p(i + r(g)) 3	944.013
	φ (i), p(i)	2	947.381

Table 2

 2 Bank Vole example with a random effect on group. i is for the intercept, r(g) is for a group random effect. N p is the number of mathematical parameters of the model.the Bank Vole, Myodes glareolus. To better understand the epidemics in humans, it is crucial to understand the dynamics of Puumala virus within the reservoir host populations. The data were collected to investigate the impact of Puumala virus on the demography of Bank Voles. Nine sites, typical for the optimal bank vole's habitat, were monitored in the endemic zone where most human cases have occurred in the last years in Belgium. In each site, six trapping sessions were carried out from 2004 to 2006 with 100 live traps in a 1 ha grid Tersago

(shed)

. Sites are coded as nine groups which are in fact geographic sites without movement among sites. Because the purpose here is not a full epidemiological study, we consider only two states {'Alive','Dead'} and we construct a simple survival model with coupled random group effects on capture. The sentence for capture (decomposed in E-SURGE as first capture=FIRSTE and recapture=NEXTE) FIRSTE+NEXTE+RANDOM(NEXTE.GROUP) builds a model where logit(p age=2,g ) = β 2 + b g .

Acknowledgements

The authors warmly thank Laurent Crespin for suggesting Application 2 and Tersago Katrien for providing the bank vole (Myodes glareolus) data and two anonymous referees for very useful comments that improved the presentation of the paper. This research was supported by a grant from the 'Jeunes Chercheuses et Jeunes Chercheurs' program of the French ANR (ANR-08-JCJC-0028-01).

Supplementary Materials

The


Seemandhar Jain

Aditya A Shastri
email: aditya.shastri@nmims.edu

Kapil Ahuja
email: kahuja@iiti.ac.in

Yann Busnel
email: yann.busnel@imt-atlantique.fr

Navneet Pratap Singh
email: navneet.diat@gmail.com

Cube Sampled K-Prototype Clustering for Featured Data

Keywords: Sampling, Cube Sampling, Clustering, K-Prototype Clustering, Principal Component Analysis, Clustering Accuracy

Clustering large amount of data is becoming increasingly important in the current times. Due to the large sizes of data, clustering algorithm often take too much time. Sampling this data before clustering is commonly used to reduce this time. In this work, we propose a probabilistic sampling technique called cube sampling along with K-Prototype clustering. Cube sampling is used because of its accurate sample selection. K-Prototype is most frequently used clustering algorithm when the data is numerical as well as categorical (very common in today's time). The novelty of this work is in obtaining the crucial inclusion probabilities for cube sampling using Principal Component Analysis (PCA).

Experiments on multiple datasets from the UCI repository demonstrate that cube sampled K-Prototype algorithm gives the best clustering accuracy among similarly sampled other popular clustering algorithms (K-Means, Hierarchical Clustering (HC), Spectral Clustering (SC)). When compared with unsampled K-Prototype, K-Means, HC and SC, it still has the best accuracy with the added advantage of reduced computational complexity (due to reduced data size).

I. INTRODUCTION

Due to the pervasiveness of social media and the subsequent transformation of the world digitally, an extensive amount of data is being generated. According to a survey [START_REF] Milenkovic | 30 Eye-opening big data statistics for 2020: patterns are everywhere[END_REF], more than 2.5 quintillion bytes of data is produced every day, which is increasing exponentially. Since not all of this data is useful, it is important to obtain some meaning from it to help in decision making.

Clustering is the technique to group together the data items having similar behavior. It is often used to get the initial perception of data and is widely used. Some of the popular clustering algorithms include K-Means, K-Modes, K-Prototype, Hierarchical Clustering (HC), Spectral Clustering (SC), Partitioning Around Medoids (PAM), Expectation Maximization (EM) Clustering, Affinity Propagation Clustering etc. [START_REF] Lloyd | Least squares quantization in PCM[END_REF]- [7]. Each of these has its limitations [6]. K-Prototype is a commonly used clustering algorithm because of its ability to handle numerical and categorical data together [8]. To achieve this, it uses a combination of the K-Means and the K-Modes algorithms and is easy to implement as well.

Clustering becomes difficult when the size of the data is very large. One of the ways to develop an efficient clustering algorithm is to sample this data. There are several sampling techniques that are commonly used. This include random sampling, stratified sampling, Vector Quantization (VQ), pivotal sampling, cube sampling etc. [START_REF] Olken | Simple random sampling from relational databases[END_REF]- [START_REF] Tillé | Sampling Algorithms[END_REF]. We focus on cube sampling because firstly it belongs to the class of probabilistic sampling that are accurate. Secondly, it is independent of the order of the data selected to compute the crucial inclusion probabilities [START_REF] Tillé | Sampling Algorithms[END_REF]. To sum up, we propose cube sampled K-Prototype clustering algorithm with specific focus on data that contains numerical and categorical features.

There have been a few attempts to perform clustering while using a sampling technique to reduce the data [START_REF] Shastri | Vector quantized spectral clustering applied to whole genome sequences of plants[END_REF], [START_REF] Shastri | Probabilistically sampled and spectrally clustered plant species using phenotypic characteristics[END_REF]. In [START_REF] Shastri | Vector quantized spectral clustering applied to whole genome sequences of plants[END_REF], authors have proposed use of VQ sampling technique along with SC. VQ is a very basic sampling technique with comparatively low accuracy. Moreover, the algorithm in [START_REF] Shastri | Vector quantized spectral clustering applied to whole genome sequences of plants[END_REF] is specifically designed for genome data. Another work is [START_REF] Shastri | Probabilistically sampled and spectrally clustered plant species using phenotypic characteristics[END_REF], where a different probabilistic sampling technique called pivotal sampling is used with SC. Although, pivotal sampling is good in terms of accuracy, this algorithm is tested on only plant phenotypic data. This paper also uses pivotal sampling with HC. However, it also has same limitations as pivotal sampling with SC. In another work [START_REF] Yang | Exploring the diversity in cluster ensemble generation: random sampling and random projection[END_REF], the authors have used random sampling with the ensemble clustering algorithm. However, they do not show the efficacy of their approach on large datasets.

As evident from above discussion, cube sampling has many advantages and it has not been coupled with any of the commonly used clustering algorithms, which is the focus of this work. As mentioned above, besides combining cube sampling with K-Prototype, we integrate it with K-Means, HC, and SC as well. The most novel contribution of this work is in computing the crucial inclusion probabilities for the cube sampling technique using the concept of Principal Component Analysis (PCA).

The proposed algorithm is evaluated with extensive experiments that are performed on six labeled datasets from UCI repository with varying sizes. Clustering Accuracy is used as a metric for comparison. We perform four sets of experiments. First, we compare four commonly used clustering algorithms (i.e. unsampled K-Prototype, K-Means [START_REF] Lloyd | Least squares quantization in PCM[END_REF], HC [START_REF] Painkra | Clustering analysis of soybean germplasm (Glycine max L. merrill)[END_REF], and SC [START_REF] Luxburg | A tutorial on spectral clustering[END_REF]), where the K-Prototype gives the best accuracy. Second, when using standard random sampling [START_REF] Olken | Simple random sampling from relational databases[END_REF], we again compare these algorithms. Here, again, K-Prototype with random sampling gives the best accuracy. Third, we perform same comparison with our cube sampling. Fourth and finally, we compare our cube sampled K-Prototype with unsampled other clustering algorithms mentioned above. Results for both the third and the fourth sets of comparisons show that our proposed algorithm performs the best. Also, out cube sampled K-Prototype has the added advantage of reduced computational complexity (since the input data is reduced by cube sampling before clustering).

The rest of this paper is organized as follows. Section II provides a detailed discussion of our proposed algorithm. Experimental setup is presented in Section III. Section IV gives the experimental results. Finally, conclusions and future work are provided in Section V.

II. PROPOSED ALGORITHM

In this section, we first explain the concept of cube sampling including our novel way of obtaining the inclusion probabilities. Subsequently, we discuss the K-Prototype clustering algorithm, followed by our cube sampled K-Prototype clustering.

A. Application of Cube Sampling

Cube sampling is a type of probabilistic sampling technique that selects a balanced sample [START_REF] Tillé | Sampling Algorithms[END_REF]. Consider a population U of size n such that U = {x 1 , x 2 , .., x n }, where x i ∈ R m for i = 1, 2, ..., n and m is the total number of features. A sample S of size N is said to be balanced if the following equality is satisfied [START_REF] Chauvet | A fast algorithm for balanced sampling[END_REF].

j∈S x j π j = i∈U x i , (1)
where π j is the inclusion probability of the j th unit and j = 1, 2, ..., N . Obtaining the inclusion probabilities of all the units in the population, denoted by π i , forms an important aspect of this probability based sampling technique. These probabilities are not sample invariant, however, they do not affect the working of cube sampling.

Cube sampling obtains the samples in two phases; the flight phase and the landing phase. The flight phase starts by obtaining the initial values of π i for i = 1, 2, ..., n (our unique contribution; discussed in the subsequent paragraphs). These initial probabilities are updated in each iteration of flight phase using a random walk. At each iteration, at least one unit is either included or excluded from the sample. At the end of this phase, inclusion probabilities of all the units are either equal to 0 (eliminated from the sample) or equal to 1 (included in the sample).

After the flight phase, if the balancing equation is not exactly satisfied or, inclusion probability of even one unit is neither 0 nor 1, then there is a need for the landing phase. The goal of the landing phase is to find a sample S such that it is almost balanced (i.e. as close to balanced as possible).

Let p be the number of units having inclusion probabilities neither 0 nor 1. Here, we obtain all the possible samples by taking the probability values of these p units both 0 and 1. As there are two possible values for each of these p units, we get a total of 2 p possible samples. Finally, we select the sample that closely satisfies the balancing condition.

Next, the working of cube sampling is explained using an example of a random population of size six. Let the initial probabilities values are π 0 = (0.5, 0.5, 0.5, 0.5, 0.5, 0.5), and we need to select three samples. The cube sampling transforms these initial inclusion probabilities as shown in [START_REF] Lloyd | Least squares quantization in PCM[END_REF]. Upon convergence, at the ϑ th step, units with the final probability as 0 are discarded and units with the final probability as 1 are selected in the sample.

π 0 =         0.5 0.5 0.5 0.5 0.5 0.5         →         0.6 0.6 0 0.6 0.6 0.6         →         0.6 0 0 1 0.8 0.6         →         1 0 0 1 0.5 0.5         →         1 0 0 1 1 0         = π ϑ (2)
Mathematically, the initial inclusion probabilities (π 0 ) are updated using the below equation [START_REF] Chauvet | A fast algorithm for balanced sampling[END_REF].

π ι+1 = π ι + λ * 1 (ι)u(ι) with probability q 1 (ι) π ι -λ * 2 (ι)u(ι) with probability q 2 (ι) (3)
Where

• q 1 (ι) = λ * 2 (ι)/ {λ * 1 (ι) + λ * 2 (ι)}, • q 2 (ι) = 1 -q 1 (ι), • u(ι) is random kernel of dataset U , • λ * 1 (ι) is the largest value satisfying 0 ≤ π ι +λ 1 (ι)u(ι) ≤ 1, and • λ * 2 (ι) is the largest value satisfying 0 ≤ π ι -λ 2 (ι)u(ι) ≤ 1.
Inclusion probability of an unit x i is defined as the probability with which the unit is selected in the sample S and, as earlier, is denoted by π i . Here, we propose a novel and an efficient way to calculate the initial inclusion probabilities using Principle Component Analysis (PCA) [START_REF] Karamizadeh | An overview of principal component analysis[END_REF]. Although PCA has been used earlier to sample data, it has never been combined with any probabilistic sampling techniques. PCA transforms the data by projecting it onto a single axis where elements are maximally deviated. The process of calculating the inclusion probabilities is described in the following equations [START_REF] Tillé | Sampling Algorithms[END_REF]:

π i = min 1, h -1 (n) x i X , i ∈ U (4)
where, as earlier, U is the total population and h is defined as

h(z) = i∈U min z x i X , 1 , (5)
with i∈U π i = N and X = i∈U x i .

B. Our Algorithm

The K-Prototype algorithm is an integration of the K-Means algorithm and the K-Modes algorithm to handle mixed data types. Here, the basic idea is similar to K-Means with a different distance criteria to handle numerical and categorical features. The aim is to minimize the cost function by partitioning the dataset into k clusters. This cost function is given as

C = k l=1 N j=1 p jl DIST (x j , q l ) , (6)
where p jl ∈ {0, 1} denotes the membership of data point x in cluster l, q l denotes the cluster center for cluster l and DIST(x j , q l ) is the distance criteria (dissimilarity measure), which is given as

DIST (x j , q l ) = mr r=1 (x jr -q lr ) 2 + γ l mt t=1 δ (x jt , q lt ) , (7)
where,

• r is index for numerical features and m r is the total number of numerical features for each data point x j , • x ir denotes the value of the r th numerical feature of x j , • q lr is the mean of the r th numerical feature in the cluster l, • γ l is the weight factor for categorical features for cluster l, • t is index for categorical features and m t is the total number of categorical features of x j , • x jt denotes the t th categorical feature of x j , • q lt is the mode of the t th categorical feature in the cluster l, and The proposed algorithm (Cube sampled K-Prototype algorithm) is given in Algorithm 1 and the broad framework is graphically represented in Fig. 1.

Algorithm 1: Proposed Algorithm

Input: A-Data Matrix (n × m rows), N -Sampling size, k-Number of clusters Output: Clusters 1 π denotes initial probabilities that is computed as discussed in Section II-A 2 π n denotes final probabilities that is obtained after the flight and the landing phase of cube sampling (also discussed in Section II-A) 3 Generate sample S from π n by using units with probability 1 4 Choose k cluster centers from S randomly to form q old 5 Find the distance between elements and cluster centers as described in (7) to form q. 6 while q old = q do 7 q old = q 8 Update q using (7) again Next, we perform the complexity analysis of our proposed approach.

1) Calculating inclusion probabilities: O(m 2 n + m 3 ) 2) Cube sampling: O(n)

3) K-Prototype algorithm: O(kN m) 4) Reverse mapping: O((n -N )k) Here, n is the number of data points in the original dataset, N is the sampling size, m denotes the number of features, and k is the number of clusters. Thus, the total time complexity of our algorithm is O(n(k + m 2 ) + kN m + m 3 ).

III. EXPERIMENTAL SETUP

In this section, we first briefly discuss the datasets used for our experiments and then describe the criteria used to check the goodness of our proposed algorithm. For experiments, we use six datasets with varying sizes. These include labeled datasets downloaded from the UCI repository [START_REF] Dua | UCI machine learning repository[END_REF]. Table I provides a brief summary of these datasets, where all the columns are self explanatory.

Since we use labeled datasets for our experiments, we check the quality of our clustering by using Clustering Accuracy (CA), which is calculated as given below [START_REF] Fahad | A survey of clustering algorithms for big data: taxonomy and empirical analysis[END_REF].

CA = k l=1 max (C l | L l ) |U | , (8)
where C l is the set of elements in the l th cluster, L l is the set of all the class labels in the l th cluster, max(C l | L l ) is the number of elements with the majority label in the l th cluster, and U is the population size.

The experiments are performed on a machine with Intel® Core™ i7-7500U CPU 2.90 GHz and 16 GB RAM. Python 3.7 is used to implement all the algorithms.

IV. RESULTS

In this section, we show the effectiveness of our proposed algorithm (i.e. Cube sampled K-Prototype). For this, we do four various comparisons here.

• Without performing any sampling, we first compare K-Prototype clustering with other commonly used clustering techniques (K-Means, Hierarchical Clustering (HC), Spectral Clustering (SC)). • When using frequently used random sampling, we again compare K-Prototype with same set of clustering algorithms. • Next, when using our proposed cube sampling, we further compare all the above-mentioned clustering algorithms.

• Finally, we compare our cube sampled K-Prototype algorithm with unsampled other clustering algorithms mentioned above. Table II compares four different clustering algorithms (K-Prototype, K-Means, Hierarchical Clustering (HC), and Spectral Clustering (SC)) without applying any prior sampling of data. From this table, it is clear that CA values for K-Prototype are greater than or equal to the other algorithms for all the datasets. Thus, K-Prototype performs best among all.

Next, in Table III, we again compare the four clustering algorithms mentioned above but on sampled data (using random sampling). Without loss of generality, we take sample size as 100 for the two relatively smaller datasets (German and Heart), while sample size is taken as 1000 for the other datasets. From this table, we again observe that K-Prototype with random sampling performs best when compared with all other algorithms with random sampling. Use of random sampling deteriorate the CA values substantially, which is not the case with our sampling technique (discussed below).

Further, in Table IV, we perform the same set of comparisons as in Table III except using data obtained after applying cube sampling. The sample sizes here are same as in the previous set of comparisons. We can observe from this table that the CA values for K-Prototype with cube sampling (our proposed algorithm) is greater than all the other algorithms with cube sampling.

Finally, we compare our proposed algorithm (cube sampled K-Prototype) with the unsampled four clusterings as discussed before. For this, we pick the data from the relevant columns of Table II and Table IV. This set of comparisons is given in Table V below. We observe from this table that our algorithm performs better than all others.

To demonstrate that our technique gives good result irrespective of the size of the sample, in Fig. 2, we plot the CA values of our proposed algorithm for different sample sizes and different datasets. On the x-axis, sample size is given, and on the y-axis, CA values are given. From these graphs, it is clear that sample size does not affect the CA values. Thus, use of cube sampling does not deteriorate the quality of clustering as done by K-Prototype. It also has an added advantage of reduced computational complexity (as it works on reduced dataset making it faster).

V. CONCLUSION

In this paper, we propose a cube sampled K-Prototype algorithm for grouping large data. The most innovative aspect of our algorithm is computation of inclusion probabilities in cube sampling using PCA. Experiments on labeled UCI datasets with varying sizes show that our proposed algorithm gives best clustering accuracy when compared with other commonly used unsampled clustering algorithms (K-Means, HC, and SC). Also, use of cube sampling does not deteriorate the accuracy of the K-Prototype algorithm. Furthermore, our algorithm also has the added advantage of reduced computational complexity (since the size of the data is reduced by cube sampling). In future, we plan to adapt our algorithm for datasets with millions of units. We also aim to combine cube sampling with other accurate clustering algorithms. We intend to adapt our algorithm for Hadoop framework to make it parallel as well [START_REF] Nemade | Scaled and projected spectral clustering with vector quantization for handling big data[END_REF].

•

  δ(a, b) = 1 for a = b and δ(a, b) = 0 for a = b.

9Fig. 1 :

 1 Fig. 1: Framework for the proposed algorithm.

TABLE I :

 I Dataset information from the UCI repository

	Dataset	# of records # of features # of classes Dataset Descriptions
	German	1000	20	2	Categorizes individuals based on a set of factors as either good
					or bad credit risks.
	Heart	303	13	2	Classifies if a person has heart disease or not.
	Satellite	6435	36	2	Predicts the classification of neighborhood using multi-spectral
					values of pixels in a satellite picture.
	Adult	45222	14	2	Classifies whether an individual earns more than 50k per year
					or not based on census data.
	Shuttle	43500	9	7	Indicates whether or not a NASA Space Shuttle component
					will fail.
	Kddcup99	24701	41	18	Simulates a range of intrusions in a military network environ-
					ment, distinguishing between types of intrusions and standard
					connections (using a set of auditable data).

TABLE II :

 II CA values for four different clustering algorithms on the six labeled datasets without prior application of sampling.

	Dataset	K-Prototype K-Means HC	SC
	German	0.68	0.53	0.62 0.66
	Heart	0.59	0.52	0.53 0.52
	Satellite	0.79	0.71	0.72 0.75
	Adult	0.69	0.61	0.67 0.68
	Shuttle	0.82	0.82	0.82 0.82
	Kddcup99	0.69	0.59	0.69 0.68

TABLE III :

 III CA values for four different clustering algorithms on the six labeled datasets with prior application of random sampling.