A model-based approach for clustering binned data

Data analysis often considers observations as continuous, but there are many situations where they are recorded with a very low precision during the data gathering process. This might happen since the data were measured with non high-precision instruments, due to their cost, or the lack of availability of an adequate measuring device. It could also be that the original data are only accesible in a summarized form, for simplicity or confidentiality reasons. In these and similar scenarios, assuming that the original observed data fall on an interval $(a,b)\subset\mathbb{R}$, they are only available after being processed as follows. First, the interval is divided into a series of non-overlapping subintervals, or bins, $T_{1},\ldots,T_{m}$ such that their union $T_{1}\cup\cdots\cup T_{m}$ covers $(a,b)$, and then, it is counted how many datapoints fall in each bin. Each bin $T_{l}$ is represented by its center $\tau_{l}$, whose value can be its midpoint or some statistics such as the mean or median value. This data binning process produces a sequence of pairs $(\tau_{1},f_{1}),\dots,(\tau_{m},f_{m})$ representing the center and frequency, or number of observations, assigned to the corresponding bin $T_{l}$, $l=1,\dots,m$. We call binned data to this kind of datasets.

Binned data play important roles in fields such as biology, ecology, social sciences, economy, engineering, computing, and others. For statisticians, perhaps one of the most known example of binned data is the Fisher’s Iris dataset [21], consisting of $50$ samples from each of three species of Iris; each sample records five values: sepal length, sepal width, petal length, petal width and species. The variables for lengths and widths are clearly continuous but were recorded using only one decimal place. Another classic example appears in [55] as an example of finite mixture distributions, originally published by [6], where the lengths of $256$ snappers are reported as a frequency table. In this way, many biological or ecological studies collect data using a binned structure; see, for example, [3] for line transect sampling, a method used for estimating wildlife populations, and the references in [48] and [15] for cases of data binning in fishery. In socio-economic studies, data related to income or poverty are mainly available in binned form for developing countries, or due to the high cost of sharing the raw data (see Hajargasht et al. 32 and Eckernkemper and Gribisch 14 among others); similarly for schooling and child welfare (see, Guo 28 and Friedman et al. 19). It also happen that organizations wish to protect data confidentiality, security, and integrity, as explained by [1], so only binned data are made public.

On a different direction, there are cases where raw data are binned in order to avoid computational burdens, or to simplify further processing, like in [29] and [51]. Also, the process of binning has been used to filter some issues associated to the data-generating phenomena [31]. Finally, data might be originally collected in a binned form, like in survival analysis, where it is common to have interval-censored data, or in pattern recognition and machine learning problems, where color and feature histograms are the input data (see Cadez et al. 9 and the references therein for more details).

Although descriptive analysis or some elemental estimations might be performed using the binned data as is, these can be inaccurate by some kind of bias inherent to the data, and aiming to perform more elaborated or specific statistical analyses become more cumbersome. Among the literature, a wide range of methodologies have been proposed for the statistical treatment of binned data. For example, in parametric estimation, there are methods based on maximum likelihood [54, 15, 14], in the method of moments [32], or following Bayesian approaches [14]. When the data exhibit a more complex distribution, kernel density estimation proposals are widely used; they can be found, for example, in [52], [53], [5], [9], [36], [58], [49], [59], and [31]. A different approach is given by [39]. Another popular class of methods along this line is based on a modification of the EM algorithm [13] to handle missing data. Two extensions are given by [44] and [7], and have been used in mixture modeling by [8], [29], [51] and [1]. With regard to fully Bayesian methods designed for binned data, [40] propose a density estimation model, and [2] and [24] estimate the parameters for finite mixture models.

Motivated by the study of marine species in Quintana Roo coasts, we develop a Bayesian nonparametric methodology for estimating the underlying clustering structure based on collected data, available only in a binned form. Providing a clustering structure may allow to identify different cohorts in the population, and estimate some characteristics of interest. This information can subsequently be incorporated in growth models, a common tools for fisheries management (see, for example, Lorenzen 43).

Our proposed methodology, on the one hand, models heterogeneity in the data following a clustering perspective by means of random partitions [17, 37, 38], in particular, a modification of the distribution over partitions for the Dirichlet process [18] is used. On the other hand, data within groups are modeled under a probabilistic framework, also known as a model-based approach [23]. However, since the data is provided in a binned form, it is necessary to extend this framework; which we tackle by treating original data as missing values. Once the underlying clustering structure is estimated, we can easily estimate group’s parameters, as well as some other characteristics, like the distribution of the data. We test the performance of our methodology using simulated and real data. Finally, we analyze samples of the length of Lobatus (Strombus) gigas, collected at different times in Punta Gavilán, Quintana Roo.

For any dataset $y=\{y_{1},\dots,y_{n}\}$ of $n$ elements, cluster analysis, also known as clustering, aims to provide an estimate of their grouping structure, say $\pi$, which partitions $y$ into nonempty subsets $\pi_{j}$, $j=1,\dots,k$. Elements belonging to the same subset, or group in the context of clustering, have more similar characteristics among them than those in any other group. Given the importance of the task of clustering, many methodologies have been developed. Some of them cluster the elements according to their distance, in some sense, to the other elements in the same group or to some representative of the group; hierarchical clustering and the $k$ means algorithm are examples of these methods. Other approaches make use of a probabilistic component to associate the elements to a single group. There, it is assumed that each element is a realization of some random variable, so all elements belonging to the same group $\pi_{j}$ were generated from the same probability distribution $g_{j}$; this approach is known as model-based [23]. A common example of model-based methods are mixture models, i.e. convex linear combinations of probability distributions.

A different probabilistic approach for cluster analysis is based on random partitions. Early works on this topic can be traced back to [17] and [37, 38] for the study of population genetics, and, nowadays, they are widely used in Bayesian nonparametric models; see, e.g. [33] and [11] for a more extensive overview. A random partition is a random object taking values over the space generated by all the possible groupings of the $n$ elements, i.e. every possible value of $\pi$. When they are combined with the model-based approach, we can determine the grouping structure underlying to a given dataset. By providing a probabilistic framework for all the possible groupings, we can have a priori and a posteriori probabilities for each possible value of $\pi$, as well as some other estimates of interest. An example of this class of models is known as product partition models [30]. Throughout this work, we refer this approach for clustering as random partition models.

A natural way to describe random partition models is through Bayesian hierarchical models. Let us denote by $\mathcal{P}$ the space of all the possible values of $\pi$, in terms of the observations’ indices, and assume that the probability distribution for each group, $g_{j}$, also called kernel, is the same and only its parameter, $\theta_{j}$, varies, so $g_{j}(\cdot):=g(\cdot\mid\theta_{j})$. Then, the hierarchical model can be written as

where $\nu_{0}$ and $\rho_{0}$ are prior distributions for the kernel parameters and the random partition, respectively.

The prior $\nu_{0}$ in Model (1) can be chosen such that, together with $g$, both form a conjugate model. On the other hand, the prior for the random partition $\pi$, $\rho_{0}$ seems less straightforward. For example, [30] defines partition distributions as proportional to the product of non-negative cohesion functions, $\kappa$, applied to each group in $\pi$, i.e.

for all partition $\pi_{1}/\cdots/\pi_{k}:=\{\pi_{1},\dots,\pi_{k}\}\in\mathcal{P}$.

A different class of partition distributions arises when working with discrete almost surely (a.s.) random probability measures (RPMs), daily used tools in Bayesian nonparametrics; see, for instance, [50], [35]. Any sample of size $n$ from an a.s. discrete RPM $\tilde{p}$, taking values in $\mathbb{X}$, is bounded to have ties with positive probability, and, as consequence, to form $K=k$ groups with distinct representative values, each one having frequency $N_{1},\dots,N_{k}$, with $1\leq k\leq n$. Therefore, a partition distribution is obtained through

by letting $\Pr(\pi=\pi_{1}/\cdots/\pi_{k})=p(\#\pi_{1},\dots,\#\pi_{k})$, with $\#A$ the cardinality of $A$. In the computation of the expression above, the expectation is with respect to the distribution of $\tilde{p}$. The resulting class of distributions is known as exchangeable partition probability functions (EPPFs), where exchangeable stands for symmetric since it can be noticed that the probability is the same under permutations of the groups. Further, a priori, it does not matter what the sampled values of $\tilde{p}$ are in order to compute the probabilities, and more importantly, it can be seen that EPPF-based models also allow us to infer about the number of groups, denoted by the random variable $K$. The Dirichlet process [18], the daily workhorse in Bayesian nonparametrics, is the canonical example for the distribution of $\tilde{p}$, and its EPPF is given by

for some $\alpha>0$, where $n=n_{1}+\cdots+n_{k}$, and $(x)_{n}=\Gamma(x+n)/\Gamma(x)$ is the Pochhammer symbol. More general constructions can be found in, for example, [41, 42], [22] and [25].

We illustrate the performance of our model by using some datasets. In these cases, a normal density is set for modeling individual groups, i.e. $g$, with unknown mean and precision, so $\theta_{j}=(\mu_{j},\lambda_{j})$, and a conjugate normal-gamma prior $\nu_{0}$ is chosen. Also, the restricted EPPF for the Dirichlet process is set for the prior over partitions, $\rho_{0}$. Therefore, the model is

The full conditional distribution for the latent variables $y_{i}$ are, then, a $T_{e_{i}}$-truncated normal distribution of parameters $(\mu_{j},\lambda_{j})$. For the kernel parameters $(\mu_{j},\lambda_{j})$, we have a normal-gamma distribution of parameters

with $\bar{y}_{j}=\frac{1}{n_{j}}\sum_{i\in\pi_{j}}y_{i}$. Additionally, a gamma prior is assigned to the total mass parameter $\alpha$ for the prior distribution of $\pi$; its posterior is a mixture of gamma distributions as explained in [16].

It is important to say that in many applications, the bins $T_{l}$, for $l=1,\dots,m$, are not provided, but only a representative value $\tau_{l}$ as explained in the Introduction. Without any further details, we assume the bins take the form $T_{l}=(t_{l-1},t_{l}]$, with $t_{l}=\tau_{l}+(\tau_{l+1}-\tau_{l})/2$, for $l=1,\dots,m-1$, $t_{0}=\tau_{1}-(\tau_{2}-\tau_{1})/2$, and $t_{m}=\tau_{m}+(\tau_{m}-\tau_{m-1})/2$.

As a first example, a synthetic dataset is used. A sample of size $n=500$ is drawn from the following mixture of normal distributions

where the second parameter of the densities is expressed in terms of the variance. Given the sample, data was binned taking the interval $(5,35)$, and $m=30$ bins where obtained, each one of length one, together with their frequencies. The histogram is depicted in Figure 1. The MCMC was run $30,000$ iterations long, discarding the firsts $20,000$. Base measure parameters were set to $(\omega,a,b)=(0,1.1,1)$ and different values for $c$ were used, in particular $1$, $0.1$, and $10$. Hyperparameters for the total mass parameter $\alpha$ are set to $(1,1.1)$. For each configuration, the posterior modal partition $\hat{\pi}$ was computed, and conditioned on it, the density and groups’ parameters where also obtained. Figure 1 also displays the estimated density for each case.

When parameter $c=1$, the estimated partition is $\hat{\pi}=(150,92,117,141)$, from which the mixing weights can be estimated by dividing each group size by the sample size, being them $(0.300,0.184,0.234,0.282)$. Moreover, conditioned on $\hat{\pi}$, the estimated means are $(8.00,15.51,23.79,29.91)$, and the standard deviations, $(1.14,2.63,2.63,3.11)$; the estimated density is the continuous line in Figure 1.

The other values for parameter $c$ allow us to see how it acts as a smoothing parameter. When $c=0.1$, the posterior modal partition contains a few groups, three, and the estimated density is coarser or softer (dotted line in Figure 1) than the case $c=1$. On the other hand, more groups are detected when taking $c=10$, the posterior mode contains nine groups, and the estimated density follows closer all bin heights (dashed line in Figure 1).

The second dataset used corresponds to the measurements of the length, in inches, of $256$ snappers [6], referred in the Introduction. Figure 2 presents the histogram. The same three sampling schemes as before were used, and a similar behaviour is observed when the value of parameter $c$ varies. For the case $c=1$, four groups are estimated, whose means are $(3.24,5.17,7.24,9.72)$. It is worth mentioning that these estimates are similar to those presented in [55].

Studies of population demographics of marine resources are important for the management and to guarantee the sustainability of exploited marine populations. Fisheries management requires high quality estimates of length- and weight-at-age to incorporate such information in harvest and optimal management models (see, for example, Petitgas et al. 48 and Ogle and Isermann 47). Growth models commonly use length-at-age relationships obtained from modal progression analysis; see, e.g. [12]. Most bony fishes produce growth rings in their scales, which can be used to obtain length-at-age keys. Other marine resources, such as molluscs and shellfishes, do not produce growth rings and the usual way to obtain the length-at-age keys is through the analysis of modes from histograms constructed from samples from commercial catches or from field surveys. Estimating the number and the location of those modes is not always trivial due to the existence of size overlapping from different cohorts. In this section, we apply our proposed model to analyze the modal progression of Lobatus (Strombus) gigas using data collected at 10 points in time in Punta Gavilán, Quintana Roo, México. The data were binned using 10-millimeter intervals (Figure 3).

The method was run for each sample, one by date, taking $10\,000$ iterations after discarding a first lot of $20\,000$. Base measure parameters were fixed $(\omega,a,b)=(0,1.1,1)$; regarding parameter $c$, it was found that it has some relationship with the sample size, so a different value was chosen for each one. Figure 3 displays, for each sample, the modal estimated density (dashed line) together with its respective mixture components (continuous lines); Table 1 contains the estimated modal partition, and estimated posterior means and standard deviations given this partition.

The histograms and the estimated densities for each sample, as presented in Figure 3, show the presence of at least two groups (modes). It is believed that each mode corresponds to what is called a “year-class” in fisheries. The continuous growth process of L. gigas can be inferred from the displacement of the modal peaks to the right of the $x$-axis. The posterior means in Table 1 show that for example, the mode at $69.88$ mm in March 2010 moved to $82.05$ mm in April 2010 and reached a value of $144.02$ mm in February 2011. A similar pattern can be tracked for other modes. Moreover, in March and September 2010, and in February 2011, a third mode appears at the left side of the histogram, which corresponds to the incorporation of a new cohort to the study area.

The modal progression obtained from the model we propose can also be incorporated in growth models, such as those by [27] and [56] (cf. Hernandez-Llamas and Ratkowsky 34 and Lorenzen 43). The resulting growth rates and parameters can be incorporated to stock and management assessment process to help resource managers to maintain a sustainable fishery for L. gigas. The method proposed here has the advantage of detecting the number of modal peaks and their location in an automatic way, unlike methods like [4], where the user has to detect the modes visually and guess their location. Other methods like kernel density estimation, the user relies on additional procedures for detecting the modal peaks.

The statistical analysis of binned data is of interest in different applied fields. We have presented a Bayesian nonparametric approach aiming, for univariate binned data exhibiting a complex structure, to perform cluster analysis together with parameter estimation for detected groups. In particular, we applied it to analyze samples of the length of L. gigas, detecting cohorts and their numerical characteristics.

Under a Bayesian perspective, the absence of the original observations was handled by a simple data augmentation mechanism for missing data, that is by incorporating an extra level in the hierarchical model. This approach places our method in the model-based paradigm, which allows to provide an estimated model for the individual groups. Furthermore, it is also possible to infer about the complete model, meaning using the complete mixture model, for example, through its estimated density. Even more, the usage of random partitions in the modeling let us to infer about the number of groups, a feature hardly available in other methodologies. Computationally, restricting their support to only consider clusters with no gaps, eases the exploration of the whole space of possible groupings during the posterior sampling.

The results obtained for the real dataset, L. gigas’ length, allows to understand the composition of the species at different times of the year. As explained, besides of automatically estimating the number of cohorts, our method do not require any extra procedure for detecting groups’ numeric characteristics, like the mean value. Also, all this information can be further explored by incorporating it in growth models, or similar tools; their results might be useful for resource management.

Finally, our model can be extended in several ways. It is straightforward handling non-homogeneous bins, i.e. bins having different lengths, or using another univariate data model for groups. Moreover, the extension to analyze data in higher dimensions mainly requires an appropriate sampling scheme of the truncated data distribution.