Estimation for finite mixture of simplex models: applications to biomedical data

Abstract

Simplex distribution has been proved useful for modelling double-bounded variables in data directly. Yet, it is not sufficient for multimodal distributions. This article addresses the problem of estimating a density when data is restricted to the (0,1) interval and contains several modes. Particularly, we propose a simplex mixture model approach to model this kind of data. In order to estimate the parameters of the model, an Expectation Maximization (EM) algorithm is developed. The parameter estimation performance is evaluated through simulation studies. Models are explored using two real datasets: i) gene expressions data of patients’ survival times and the relation to adenocarcinoma and ii) magnetic resonant images (MRI) with a view in segmentation. In the latter case, given that data contains zeros, the main model is modified to consider the zero-inflated setting.

Keywords

simplex distribution Finite mixture zero-inflated models simulation EM algorithm MRI

1 Introduction

Some statistical processes can be assumed to result from a mixing of several known or unknown distributions. The natural representation of these phenomena is called mixture of distributions, and it attempts to recover the information of each component in the mixture model. This recovering step is generally translated into the estimation of the parameters of those distributions. Often such estimation is the main aim of the investigation, but sometimes, the purpose of the research may be to scrutinize the original distribution of each particular observation. This step is the main goal of the model-based clustering or classification.

In model-based clustering, the empirical distribution is modelled by a finite mixture of distributions. Usually, Gaussian ones are the most ubiquitously considered (Carreira-Perpiñán, 2000; McLachlan and Peel, 2000; Contreras-Reyes and Cortés, 2016), and consequently, this kind of mixtures lie in the real space. In practice, nevertheless, most variables are bounded because they are subject to the impossility of infinite precision in measuring devices. We are concerned about this limited aspect of mixture modelling. Though the literature is limited about variables bounded in an interval, much effort has been made recently to manage variables limited to the $(0, 1)$ interval. Perhaps the most famous examples concern the beta and the uniform distributions. These distributions have even explicitly been mixtured in different contexts, as in Markitsis and Lai (2010) and Bayes et al. (2012). Further instances exist in which bounding is done in just one direction. Commonly known distributions with these features are the non-negative continuous distributions. For instance, Bagnato and Punzo (2013) applied mixtures of unimodal beta and gamma components for recovery rates to Italian bank loans and mother's age at first birth data in the United States, respectively. In the same vein, several branches of bounded distributions have been explored. Furthermore, dispersion models (DMs) are a family which is worth considering because it contains some examples of bounded variables such as simplex and von Mises distributions (Jørgensen, 1997; Artes and Jørgensen, 2000). In recent studies, DMs have been considered in their multivariate and discrete forms (Jørgensen and Martínez, 2013; Jørgensen and Kokonendji, 2016).

We present a distribution mixture model restricted to the $(0, 1)$ interval, using the simplex distribution as base of components. We exemplified this model with two real data: i) A gene expressions data of survival times of patients and their relation with adenocarcinoma and ii) magnetic resonant images (MRI) data of brain tumours, commonly used in segmentation. Traditionally, for medical image analysis, most of the researchers have typically considered mixture of Gaussian distributions (Ma and Leijon, 2009; McLachlan and Peel, 2000).

This article proceeds as follows: Section 2 briefly describes the main features of simplex distribution and the way we are going to compare models. Section 3 resumes the classical expectation maximization (EM) algorithm for mixtures model. Section 4 shows some simulations to evaluate the performance of the estimation. Later in the same section, two real data applications are presented and some conclusions are given in Section 5.

2 Background and related work

2.1 Dispersion models

The development of the simplex distribution model is highly related with DMs. The (reproductive) dispersion model is represented by densities taking the form

p (x; μ, λ) = a^{*} (λ, x) e^{λ d (x; μ)},

(2.1)

where

(λ, μ) \in Λ \times Ω \subseteq ℝ_{+} \times ℝ

. Usually, DMs are denoted with

X \sim DM (μ, σ^{2})

and

σ^{2} = 1 / λ

. The function

d (\cdot; \cdot)

is called the unit deviance and satisfies

d (x; x) = 0

\forall x \in Ω

and

d (x; μ) > 0

\forall x \neq μ

. The DM (2.1) is called regular if it meets

\frac{\partial^{2} d (μ; μ)}{\partial μ^{2}} > 0, \forall μ \in Ω .

The advantages of the parametrization of this family distribution are similar to those that natural exponential models (McCullagh and Nelder, 1989) offer in the way that parametrization allows to recover a significant amount of information from the general expression. From function $d (x; μ)$ , we can derive important quantities such as the unit variance function defined as

V (μ) = \frac{2}{\frac{\partial^{2} d (μ; μ)}{\partial μ^{2}}} .

(2.2)

In addition, DMs can be divided into several categories, for example, the (regular) proper dispersion models when the function

a^{*} (λ, x)

is factored by

a (λ) b (x)

. Furthermore, the proper DMs can be reparametrized to obtain the expression

p (x; μ, σ^{2}) = a (σ^{2}) V^{- \frac{1}{2}} (x) \exp \{- \frac{1}{2 σ^{2}} d (x; μ)\} .

(2.3)

2.2 Simplex distribution

The simplex distribution can be obtained by conditioning $m$ independent generalized inverse Gaussian random variables on their sum (Barndorff-Nielsen and Jørgensen, 1991). This has the representation

p (x; μ, σ^{2}) = [2 π σ^{2} {x (1 - x)}^{3}]^{1 / 2} \exp \{- \frac{1}{2 σ^{2}} \frac{(x - μ)^{2}}{x (1 - x) μ^{2} (1 - μ)^{2}}\},

(2.4)

which is denoted by

S^{-} (μ, σ^{2})

, and its parameters are

μ

μ \in (0, 1)

and

σ^{2}

σ^{2} > 0

x \in (0, 1)

. These parameters are interpreted as location and dispersion parameters, respectively (Jørgensen, 1987). Simplex distribution serves well to handle bounded variables because it can adopt different shapes and asymmetries. Specifically, it can be used to directly model bimodality. This is the case where the dispersion parameter is relatively large. Moreover, if the location parameter is close to 0.5, each of the two modes will be around 0 and 1. Yet, under these circumstances, the estimation of the model could be problematic and questionable (Song, 2007, Chapter 2). Therefore, as we are not interested in mixturing bimodal distributions or with large dispersions, we have discarded them as individual components in the finite mixture modelling.

2.3 Model comparison through Vuong's test

In this work, we lean to a directional and symmetric test developed by Vuong (1989) for choosing between models in a non-nested context (Greene, 2012). Vuong's method is helpful to test if two underlying densities have the same Kullback–Liebler information criterion (Kullback and Leibler, 1951). The test takes classical likelihood ratio statistic theory as basis and computes the difference of maximum log-likelihood values for two competing normalized models. Specifically, let $L_{i, 0}$ and $L_{i, 1}$ be the $i$ th contribution to the likelihood function under two competitive models, say M $_{0}$ and M $_{1}$ , respectively, and let $m_{i} = \ln L_{i, 0} - \ln L_{i, 1}$ be its difference. Then, the statistic

V = \frac{\sqrt{n} (\frac{1}{n} \sum_{i = 1}^{n} m_{i})}{\sqrt{\frac{1}{n} \sum_{i = 1}^{n} (m_{i} - \bar{m})^{2}}}

(2.5)

follows a standard normal distribution when the models do not show differences, that is, when the expected quotient between the likelihood contributions tends to be 1. Consequently, if the value of

V

is greater than a constant

c

, the test recommends choosing M

_{0}

. Equivalently, if

V < - c

, M

_{1}

should be selected. Values lying in

| c |

indicate no preference for either model. Traditionally,

c \approx 1.96

for a 5% of significance level under normal assumption.

3 Proposed method

3.1 Mixture distribution

According to McLachlan and Peel (2000), $f (y; θ)$ follows a mixture of distributions if its density function can be written as

f (y; θ) = \sum_{j = 1}^{K} π_{j} f_{j} (y; θ_{j}),

(3.1)

where

f_{j} (y; θ_{j})

and

j = 1, \dots, K

are called components of the mixture and depend on parameters in

θ_{j}

. Parameter

π = (π_{1}, \dots, π_{K})

is the mixing proportion, frequently called weight and is defined under constraints

\sum_{j = 1}^{K} π_{j} = 1

and

π_{j} \in [0, 1]

j = 1, \dots, K

. Note that the mixture distribution depends on a set of parameters in

θ = (θ_{1}, \dots, θ_{K}, π)

. We are interested in determining the value of the parameters for each component from a data sample,

y = (y_{1}, \dots, y_{n})^{⊤}

, when the components

f_{j} (y; θ_{j})

obey simplex distributions (see equation 2.4), that is,

θ = ((μ_{1}, σ_{1}), \dots, (μ_{K}, σ_{K}))

, assuming the number of components

K

is known. To model the distribution, classifying observations into several groups may be imperative and is called mixture models clustering. Figure 1 shows some densities of mixture of simplex distributions.

Figure 1:

Some simplex mixture densities.

3.2 Estimation through the EM algorithm

In this section, we briefly describe the procedure to estimate the mixture distribution parameters using maximum likelihood method. Traditional log-likelihood function for model (3.1) has the form

\log L (θ) = \sum_{i = 1}^{n} \log \{\sum_{j = 1}^{K} π_{j} f_{j} (y_{i}; θ_{j})\} .

(3.2)

However, maximization of this kind of function is unstable and prone to failure. More robust results can be achieved by using the EM algorithm (Dempster et al., 1977; McLachlan and Peel, 2000). We consider

y_{c}

as an augmented dataset, that is,

y_{c}

is a pair of variables

y_{c} = (y, z)

, where

z_{ik} = 1

when

i

th subject comes from the latent

k

th component and 0 otherwise. The augmentation setting is also known as completing data (McLachlan and Peel, 2000; Pawitan, 2001). As a result, the complete log-likelihood function is turned into

\log L_{c} (θ) = \sum_{i = 1}^{n} \sum_{j = 1}^{K} z_{ij} (\log π_{j} + \log f_{j} (y_{i}; θ_{j})),

where

z

is treated as a missing value vector. Due to this structure, it is possible to use the EM algorithm. Essentially, the algorithm's nature consists of the alternation of expectation and maximization steps:

· E-step: Given the observed data, $y$ , compute the conditional expectation of equation (3.2). For an unknown set $θ$ , we considered the initial value $θ^{(0)}$ . Then, the first iteration of EM algorithm is computed as ${\hat{p}}_{ij}^{(0)} = E_{θ^{(0)}} [\log L_{c} (θ) | y, θ^{(0)}]$ . For $k = 1, 2, \dots$ , we computed the following expected value

\begin{matrix} {\hat{p}}_{ij}^{(k)} & = & E_{θ^{(k)}} [\log L_{c} (θ) | y, θ^{(k)}] = \frac{π_{j}^{(k)} f_{j} (y_{i}; θ_{j}^{(k)})}{\sum_{i = 1}^{K} π_{i}^{(k)} f_{i} (y_{i}; θ_{i}^{(k)})} . \end{matrix}

· M-step: Once the approximated values ${\hat{p}}_{ij}^{(k)}$ have been obtained, we proceed with the weighted log likelihood

Q (θ; θ^{(k + 1)}) = \sum_{i = 1}^{n} \sum_{j = 1}^{K} {\hat{p}}_{ij}^{(k)} \{\log {\hat{π}}_{j}^{(k + 1)} + \log f_{j} (y_{i}; θ_{j}^{(k)})\},

(3.3)

and the vector

π

is estimated independently by

{\hat{π}}_{j}^{(k + 1)} = \sum_{i = 1}^{n} {\hat{p}}_{ij}^{(k)} / n

j = 1, \dots, K

. Therefore, by maximizing equation (3.3) with respect to

θ_{j}

, the update

θ_{j}

is obtained by

{\hat{θ}}_{j}^{(k + 1)} = {argmax}_{θ_{j}} Q (θ; θ^{(k)}) .

However, since the solution of equation (3.2) is not available in a closed form, we use a numerical approximation to get the updated estimates. Among several available methods, we settled for a finite-difference scheme to obtain the derivatives, such as the general purpose Nelder–Mead's optimization algorithm. Another possible option presents the BFGS quasi-Newton method (Nocedal and Wright, 2006). Finally, iterations must repeatedly alternate between the E- and M-step until the difference $| L ({\hat{θ}}^{(k + 1)}) - L ({\hat{θ}}^{(k)}) |$ be less than a prefixed small value $ϵ$ . In particular, we have obtained satisfactory results by using $ϵ = 0.001$ , and those were practically the same when experimenting with the BFGS and Nelder–Mead algorithms. For further details on EM algorithm's application to finite mixture models, see McLachlan and Peel (2000).

As in the generalized linear model, a transformation on parameters is often defined using strictly monotonic and twice differentiable link functions that map parameter space into ℝ. Then, to connect the structural part of the model with the predicted parameter, we have chosen the logit and logarithm functions for location and dispersion parameters, respectively. We have not attempted to measure the efficiency of this selection, although some studies have been conducted in this regard (Andrade, 2007). This selection also provides comparative benefits because some current alternative modelling softwares use those functions as default links (Grün et al., 2012).

Furthermore, we used as starting values for the location parameters, the means produced by an exploratory $k$ -means clustering analysis, $μ_{j}$ , $j = 1, \dots, K$ , (McLachlan and Peel, 2000 §2.12). In earlier stages of the research, we had tried from valid random initial values, but results were sensible to particular choices. This is, however, a very known issue for mixture models when maximizing the likelihood function can even produce consistent estimates, when function has a global maximum but otherwise just offer the local maximizing value (Dempster et al., 1977).

On the other hand, we have used the moment estimator for dispersion parameters (Song, 2007)

{\hat{σ}}_{j}^{2} = \frac{1}{n_{j} - 2} \sum_{i = 1}^{n_{j}} d (y_{ij}; {\hat{μ}}_{j}),

(3.4)

where

n_{j}

represents the size of cluster

j

{\hat{μ}}_{j}

the estimation of the mean and

y_{ij}

denotes the

i

th observation for the

j

th group obtained by the

k

-means analysis. In addition, we have included a computational flag over the dispersion parameters such that when some of them will be greater than 10, all estimation process will be restarted afresh until all components have a smaller dispersion. This also has, as a consequence, the exclusion of bimodal distributions. In the following experiences, restartings have occurred in less than 1% of the estimations.

4 Experimental results

4.1 Simulation

To evaluate the performance of simplex mixture distribution estimation, we use a Monte Carlo simulation as outlined further. The mixtures are simulated with 2, 3, 4 and 5 components using the same parameters of Figure 1, and then they are estimated by carrying out the EM algorithm described in Section 3.2, using particularly the BFGS algorithm for the maximization step. The simulations were repeated for $r = 1 000$ times and for three different sample sizes, $n = 50, 100, 500$ . To measure the quality of the estimations, the mean squared error (MSE) and the relative bias (RB) are used. In particular, we have used the expressions

MSE = \frac{1}{r} \sum_{j = 1}^{r} ({\hat{θ}}_{j} - θ)^{2} and RB = \frac{| \bar{θ} - θ |}{θ} \times 100 %, \bar{θ} = \frac{1}{r} \sum_{j = 1}^{r} {\hat{θ}}_{j} .

Note that RB and MSE are expressed in percentage and squared units, respectively. In addition, classical quantities for the classification problem, such as percentage of overall accuracy (% OA), sensitivity (Sens.) and specificity (Spec.), are also considered. Furthermore, in order to measure the tightness of data and estimated density model, we use the Hellinger distance (HD). Particularly, for a pair of distributions with densities $f$ and $g$ , HD is defined as

HD (f, g) = \frac{1}{2} \int_{- \infty}^{\infty} {(\sqrt{f (x)} - \sqrt{g (x)})}^{2} dx

(Nikulin, 2001). HD has an easy interpretation, it ranges inside

[0, 1]

HD = 0

when both densities are equal and, oppositely,

HD = 1

when both densities are absolutely different. In this article, we use as

f

the empirical distribution of data (estimated using a kernel density method) and as

g

the simplex mixture distribution whose parameters have been estimated from data.

The simulation results are summarized in Table 1. The second column shows true parameter values: the first half are locations and the second half dispersion parameters. As expected, the estimation tends to be more precise, in RB terms, when sample size increases. Nonetheless, the dispersion parameters show the most inaccurate estimation. For instance, the parameter whose real value is 0.1 in the fifth mixture model, was estimated more than two times beside its correct value. This is represented under the respective RB values, 164.57, 154.60 and 142.18%.

On the other hand, classification results (% OA, Sens. and Spec.) perform pretty well when number of clusters are of moderate size. For instance, the 93% of individual membership was correctly recovered, on average, for three clusters of 50 size each. Generally, the values are over 70% of correct classification. The same behaviour is present in the HD measure where, moreover, in all cases, the estimation are as closer or less than 30%. However, mixture model rate classification decreases largely when the number of components jumps from fourth to fifth group. In our experience, this distortion is due to the numerical method's inability to manage different and increasing sources of dispersion (Yao, 2010).

Table 1:

Estimation (Estim.), relative bias (RB) and MSE of simulations

		$n = 50$			$n = 100$			$n = 500$
$K$	True	Estim.	RB	MSE	Estim.	Rel.	MSE	Estim.	RB	MSE
	0.9	0.90	0.06	0.00	0.90	0.02	0.00	0.90	0.02	0.00
	0.6	0.60	0.04	0.00	0.60	0.02	0.00	0.60	0.07	0.00
2	1.0	0.96	3.94	0.01	0.99	0.52	0.00	1.00	0.13	0.00
	0.5	0.49	2.66	0.00	0.49	2.15	0.00	0.50	0.33	0.00
	(% OA, Sens., Spec., HD)	(100, 1.00, 1.00, 0.23)			(100, 1.00, 1.00, 0.20)			(100, 1.00, 1.00, 0.15)
	0.9	0.90	0.38	0.00	0.90	0.05	0.00	0.90	0.05	0.00
	0.8	0.79	1.40	0.00	0.80	0.32	0.00	0.80	0.02	0.00
	0.6	0.60	0.56	0.00	0.60	0.41	0.00	0.60	0.04	0.00
3	1.0	1.01	0.92	0.05	0.98	2.33	0.02	0.98	1.96	0.00
	0.5	0.56	12.34	0.04	0.57	13.56	0.03	0.50	0.70	0.00
	0.5	0.47	6.15	0.00	0.48	3.09	0.00	0.50	0.69	0.00
	(% OA, Sens., Spec., HD)	(93, 0.97, 0.89, 0.13)			(95, 0.98, 0.92, 0.11)			(98, 0.99, 0.97, 0.09)
	0.9	0.90	0.55	0.00	0.90	0.20	0.00	0.90	0.03	0.00
	0.8	0.78	2.25	0.00	0.79	0.98	0.00	0.80	0.03	0.00
	0.6	0.59	1.43	0.00	0.60	0.54	0.00	0.60	0.06	0.00
	0.4	0.40	0.04	0.00	0.40	0.03	0.00	0.40	0.01	0.00
4	1.0	1.03	2.68	0.05	1.01	1.07	0.03	0.99	0.79	0.00
	0.5	0.55	10.76	0.03	0.55	9.32	0.02	0.51	1.65	0.00
	0.5	0.47	5.86	0.01	0.48	4.09	0.01	0.50	0.62	0.00
	0.1	0.10	1.25	0.00	0.10	0.53	0.00	0.10	0.01	0.00
	(% OA, Sens., Spec., HD)	(92, 0.97, 0.89, 0.26)			(96, 0.98, 0.94, 0.25)			(98, 0.99, 0.97, 0.22)
	0.9	0.89	1.01	0.00	0.89	1.18	0.00	0.89	0.77	0.00
	0.8	0.77	3.50	0.01	0.77	4.00	0.01	0.78	2.87	0.00
	0.6	0.61	1.41	0.01	0.60	0.08	0.01	0.59	2.40	0.00
	0.4	0.46	13.80	0.01	0.45	13.33	0.01	0.43	8.32	0.01
	0.3	0.31	4.17	0.00	0.31	3.44	0.00	0.31	4.85	0.00
5	1.0	1.07	6.97	0.07	1.07	6.66	0.06	1.06	5.69	0.04
	0.5	0.55	10.67	0.07	0.54	8.30	0.06	0.50	0.96	0.00
	0.5	0.47	6.77	0.08	0.44	11.48	0.07	0.42	16.06	0.02
	0.1	0.26	164.57	0.11	0.25	154.60	0.09	0.24	142.18	0.12
	1.0	0.93	6.56	0.05	0.94	6.00	0.02	0.99	1.25	0.00
	(% OA, Sens., Spec., HD)	(65, 0.94, 0.61, 0.21)			(65, 0.89, 0.44, 0.18)			(67, 0.92, 0.57, 0.16)

Table 2:

Estimation of the gene expression models

Parameter	Simplex	Beta	Normal
Loc. or mean.	0.458 7 (0.001 5) 0.675 0 (0.002 6)	0.447 5 (0.064 8) 0.657 5 (0.019 4)	0.465 1 (0.003 9) 0.726 9 (0.090 7)
Disp., Prec. or St. dev.	0.875 6 (0.008 9) 1.373 1 (0.020 1)	26.144 6 (2.378 4) 12.874 7 (1.101 1)	0.102 0 (0.004 9) 0.088 2 (0.005 9)
π	(0.68, 0.32)	(0.45, 0.66)	(0.25, 0.75)
log-likelihood	$1 300.52$	$- 1 170.77$	$- 3 805.01$
AIC	$- 2 593.03$	2 282.03	7 618.02
BIC	$- 2 565.54$	2 309.51	7 645.51
Vuong's test	$-$	$21.96$	$54.80$
HD	0.024	0.020	0.029

4.2 Applications

4.2.1 Gene expression data

Ji et al. (2005) used a set of genes to investigate the relation between the survival times of patients and adenocarcinoma. The authors arbitrarily selected one gene, tyrosinehydroxylase (TH), from the original list and computed 7 128 correlation coefficients. After applying the transformation $y_{i} = (x_{i} + 1) / 2$ , $i = 1, \dots, 7128$ , to generate correlations between 0 and 1, they modelled it with a mixture of beta distributions.

We take the same gene to compare several models: mixture of simplex, mixture of beta and mixture of normal using two components, as in the original work. Results of the fitted models appear in Table 2. Elapsed times for estimations were 33.37, 323.58 and 1.29 seconds, respectively. The process was obtained on an Intel Core i5 processor and 7.7 GB of random-access memory running GNU/Linux (kernel 4.10). The table has both pairs of parameters for three models, that is, for simplex distribution location and dispersion parameters, for beta distribution mean and precision, and the mean and variance parameters for normal distribution. Location parameters are similar across the models. Estimated precision parameters of mixture beta distribution are notably high and, according to their standard deviations, they are not significant. Additionally, Vuong's test recommends the simplex over the beta selection model. Values of likelihoods for each model are added to Table 2. However, Ji et al. (2005) have found in practice that most traditional model selection criteria, mainly based in log-likelihood values, do not discriminate correctly the size of $K$ under beta-based models. Popular criteria, such as AIC and BIC (also appended to Table 2), are defined in terms of log-likelihood values. Specifically,

AIC (\hat{θ^}) = - 2 L (\hat{θ^}) + 2 q and BIC (\hat{θ^}) = - 2 L (\hat{θ^}) + q \log (n),

where

\hat{θ^}

is the maximum likelihood estimator of

θ^

L (θ^)

is the log-likelihood function,

n

is the sample size and q is the number of estimated parameters (Burnham and Anderson, 2002).

Vuong's test relates to the first column model, that is, simplex mixture is the M $_{0}$ (see equation 2.5). Accordingly, the simplex model best explains the gene correlations (Figure 2). Simplex mixture captures reliably the shape of the distribution around the second component on the right. Towards the centre of data, beta mixture density is closer to the histogram than the simplex model. Normal mixture also reports satisfactory results. Note that, regardless of test results, the three models produce very similar approximations of the densities, as is reflected in the HD (0.024, 0.020 and 0.029).

4.2.2 Glioblastoma multiforme data

We considered images of 30 different patients diagnosed with glioblastoma multiforme (GBM) taken from the Cancer Imaging Archive (Clark et al., 2013; Smith et al., 2016). Each patient has at least one image from an MRI procedure, making for a total of 122 images. Images have been converted into greyscale and scaled to live in the unit interval of

100 \times 100

pixels size. However, as images data include a few zero values, it was considered a fairer model turning the first component of equation (3.1) into a degenerate distribution. This converts the model into the following zero-inflated mixture simplex model (ZIMS)

g (y_{i}; θ) = \{\begin{matrix} π_{1} & if y_{i} = 0 \\ π_{2} f_{2} (y_{i} | θ_{2}) + \dots + π_{K} f_{K} (y_{i} | θ_{K}) & if y_{i} > 0 . \end{matrix}

Figure 2:

Fitted results in gene expression example: _: simplex, : beta, : normal.

The same modification was also made with the beta distribution mixture to get a zero-inflated mixture beta distribution model (ZIMB; see details of zero-inflated beta distribution in Ospina and Ferrari, 2010). Models were fitted for all the images, scanned from

K = 2

to 6 clusters, and results were compared using the log-likelihood and Vuong's test. Note that the goal of these fittings is to recover the underlying density that is assumed to explain the data. Other approaches, not considered here, try to remove the noise in a corrupted image conserving most of the important details (Nair et al., 2008). Table 3 shows Vuong's test where the ZIMS and ZIMB models

Table 3:

Totalized Voung's test for GBM images. Greatest log-likelihood results in parenthesis. Median Hellinger Distance (HD) values are appended

	K
	2	3	4	5	6
# $V \leq - 1.96$	52(49)	51(26)	60(21)	61(11)	68(15)
# $V \geq 1.96$	69(75)	66(20)	59(14)	59(6)	51(7)
# $V \in (- 1.96, 1.96)$	1	5	3	2	3
HD beta	0.34	0.34	0.37	0.37	0.40
HD simplex	0.30	0.30	0.32	0.31	0.33

Figure 3:

Three GBM images (patients) at 100 $\times$ 100 pixels size: (a), (b) and (c) refer to an original image, ZIMS and ZIMB model fit, respectively, for $K = 2$ . Images represented in batches (d)–(f) and (g)–(i) follow the same guidelines for $K = 3$ and $K = 4$ , correspondingly.

are compared. The value 52(49) in the first row mean that, for models with two clusters, 52 times the simplex-based distribution model produced better adjustment (as earlier, simplex mixture-based model is the M

_{0}

) and under the log-likelihood criterion, 49 times the models with two clusters stood out among all models with a different number of clusters under examination. For

K = 2

, the beta-based model produced almost 50% more than the simplex model (75). This apparently induces a simplification of data structure into only two clusters. Dubious cases are relatively few for all numbers of clusters examined (i.e., 1, 5, 3, 2 and 3). Figures 3 and 4 show five MRI images for different patients, where the simplex-based model is supported by Vuong's test. Note that for

K = 3

, fitted values are almost indistinguishable (Figures 3e–3f). For the fifth patient, however, beta results look better than the simplex-based model (Figures 4e and 4f, respectively). Normal mixture (omitted) performed poorer than ZIMS and ZIMB models, with 1.3% in Vuong's test (where normal mixture was superior). Similarity between observed and estimated densities are displayed in Table 3 through the HD measure. For both kind of models, although at different rates, it seems that HD increases as number of clusters increases.

Figure 4:

Continuation of Figure 3. (a)–(c) represent a GBM image (patient) for K = 5 and its estimation under ZIMS and ZIMB. This is the case when K = 6 is shown in (d)–(f).

In the next section, an individual image taken from Burzynski et al. (2014) in the treatment of an inoperable, recurrent, persistent or progressed GBM is explored in detail.

4.2.3 Burzynski GBM image

An MRI was undertaken by patient number 4 of Burzynski et al.(2014) before the treatment of ‘pazopanib and other agents, in combination with administered sodium phenylbutyrate’, and is shown in Figure 5a. There were estimated five zero-inflated models mixturing beta and simplex distributions, respectively (Table 4). Estimation times have ranged from 35 to 7 194.02 seconds for simplex models and from 98 to 2 039.33 seconds for beta models. According to the greatest log-likelihood value, we chose the model with two components.

Table 4:

Results for different ZIMS and ZIMB models of the Burzynski image example

		K
		2	3	4	5	6
	log-likelihood	4 641.65	2 183.39	3 531.34	1 559.03	2 407.73
ZIMS	AIC	$- 9 275.30$	$- 4 354.78$	$- 7 046.67$	$- 3 098.05$	$- 4 791.47$
	BIC	$- 9 246.45$	$- 4 311.52$	$- 6 988.99$	$- 3 025.95$	$- 4 704.94$
	log-likelihood	4 725.19	4 638.98	3 794.34	4 275.90	1 570.03
ZIMB	AIC	$- 9 442.39$	$- 9 265.95$	$- 7 572.68$	$- 8 531.80$	$- 3 116.07$
	BIC	$- 9 413.55$	$- 9 222.69$	$- 7 515.00$	$- 8 459.70$	$- 3 029.54$
	Vuong's test	$- 0.75$	$- 52.44$	$- 3.69$	$- 66.87$	26.36

In addition, a ZIMB mixture distribution with the same structure is estimated, as well as a normal mixture model with two components (see Table 5). Based on this, it can be noted that one of the two groups has very similar values in mean or location value (around $0.31$ ), while the other component values are not so similar across models (that is, $0.18, 0.02, 0.21$ ). In relation to dispersion and precision estimates of inflated models, the first component has a very high value for ZIMB, and a moderate value for ZIMS model. Components of weight parameter, $π$ , are rather dissimilar.

Table 5:

Results for Burzynski image example for K = 2

Parameter	ZIMS	ZIMB	Normal Mixture
Loc. or mean.	0.188 4 (0.018 0)	0.021 1 (0.007 6)	0.213 6 (0.002 5)
	0.316 0 (0.002 9)	0.312 5 (0.007 7)	0.314 8 ( $8 \times 10^{- 4}$ )
Disp., Prec. or	6.307 1 (0.010 0)	424.109 0 (0.032 2)	0.199 5 (0.001 7)
St. dev.	0.440 5 (0.010 2)	8.913 5 (0.015 3)	0.033 8 ( $9 \times 10^{- 4}$ )
$π$	$(0.02, 0.48, 0.50)$	$(0.02, 0.35, 0.63)$	$(0.69, 0.31)$
Vuong's test	$-$	$- 0.75$	$59.91$
(% OA, Sens., Spec., HD)	(62, 1.00, 0.18, 0.36)	(28, 1.00, 0.10, 0.31)	(49, 1.00, 0.14, 0.15)

Figures 5b–5d suggest the use the ZIMS model because it recovers the tumour area. It is interesting to note that the Vuong's test criterion gives a non-conclusive result between the first and second model ( $- 0.75$ ). The log-likelihood value for the ZIMB model was $4 725.19$ , also suggesting the selection of the beta mixture model.

Figure 5:

(a): Original Burzynski image at 100 × 100 pixels size. Models with two components, (b) ZIMS estimation, (c) ZIMB estimation and (d) mixture normal estimation. Models with four components: (e) ZIMS estimation and (f) ZIMB estimation.

Some other explorations are made. In particular, the second model with higher log-likelihood (3 531.34) is that with four components (Figures 5e–5f). Vuong's test between inflated models yields $- 3.69$ , indicating the beta model should be preferred over the simplex mixture model. Fitted models also points out that both models capture new groups (Table 6). However, none of those models jointly retain all those clusters, that is, the cluster between $0.1$ and $0.3$ with mean $0.25$ is found by the ZIMS model, whereas the beta mixture model catches a cluster with mean $0.51$ .

Table 6:

Results for Burzynski image example for K = 4

Parameter	ZIMS	ZIMB	Normal Mixture
	0.020 2 (0.005 9)	0.062 4 (0.016 3)	0.018 4 ( $2 \times 10^{- 4}$ )
	0.245 9 (0.051 0)	0.314 9 (0.003)	0.116 1 (0.001 9)
Loc. or mean.	0.274 (0.015 3)	0.437 2 (0.030 0)	0.312 6 ( $7 \times 10^{- 4}$ )
	0.316 1 (0.002 8)	0.507 1 (0.004 9)	0.474 1 (0.012 9)
	1.722 1 (0.017 3)	18.820 6 (0.026 5)	0.008 2 ( $2 \times 10^{- 4}$ )
	9.548 9 (0.028 3)	98.281 8 (0.019 5)	0.054 5 (0.001 8)
Disp., Prec. or St. dev.	2.566 9 (0.012 6)	2.966 7 (0.037 7)	0.043 9 ( $7 \times 10^{- 4}$ )
	0.401 4 (0.010 7)	582.070 9 (0.082 9)	0.197 3 (0.005 4)
$π$	$(0.02, 0.06, 0.17, 0.31, 0.44)$	$(0.02, 0.03, 0.10, 0.33, 0.52)$	$(0.69, 0.31)$ $(0.69, 0.31)$
Voung test	$-$	$- 3.69$	$74.82$
(% OA, Sens., Spec., HD)	(78, 0.98, 0.24, 0.46)	(35, 0.810, 0.00, 0.08)	(98, 1.00, 0.84, 0.15)

Figure 6:

Burzynski image example: Models with two components: _ zero-inflated mixture simplex estimation. : zero-inflated mixture beta estimation. ____: mixture normal estimation

Assuming that the tumour area in this particular example lie between the values 0.45 and 0.60, we are able to estimate the performance measures as in the simulation section. They are shown at the last line of Tables 5 and 6. Figures 6 and 7 offer the approximation to the histogram of this image under $K = 2$ and 4.

Figure 7:

Burzynski image example: Models with four components: : zero-inflated mixture simplex estimation. : zero-inflated mixture beta estimation. __: mixture normal estimation.

5 Conclusion

In this article, we proposed a new model for data that is restricted to the $(0, 1)$ interval and with several modes. Generally, this problem has been addressed using normal mixture distribution, the most widely used method, albeit not conceptually incorrect. To address this issue, we modified the mixture structure to incorporate a more convenient family distribution, such as the simplex distribution. To estimate it, an EM algorithm for the simplex mixture model was developed, whose overall performance was evaluated by executing several simulations. We conducted our simulation under the advice of not using high dispersion parameter values because these will produce bimodal distributions. Regarding simulation setting, the results were good. More precisely, classification rates for up to four components were over 90% of successful classification. However, when five components are considered, the overall accuracy drops by up to 70%. Additionally, precision parameters were the most affected in the estimation procedure. On the other hand, location parameter estimates were not far from real values in a relative magnitude greater than 15%, while both situations improve augmenting the sample size. Thus, we recommend the use of simplex mixture models when dimension of data, and suspected number of clusters, are of moderate size.

The proposed methods were applied to two kinds of real datasets. The first application was concerned with an already worked dataset about correlation among genes. Our proposed model provides slightly better adjustment than the referenced models, according to log-likelihood values and Vuong's test. In the second application, we worked with medical images. Colours in images are naturally limited to an interval and, in some cases, it is possible to get values at extremes of these intervals. For instance, in a greyscale image, the image can contain absolutely white or black colour in a pixel. Therefore, we have modified the original model of equation (3.1) to allow zero values to be treated within images. Depending on the situation, the model can be modified to include other inflated values. The estimation times for presented simplex mixture models have been mostly similar to the established beta models, although both have notably increased when zero-inflated values have been incorporated. Normal mixtures, which are not correct for this kind of data, have times of computation really low, around 10 seconds for the most expensive data.

This study has not been focused on proving the superiority of one model over any competitors. Rather, the goal has been its introduction and its assessment using simulation and real data. Compared to the widely applied beta model, results for these examples pointed out that no one model clearly captures all modes and clusters in the data. Evaluated criteria do not lead to a generally superior model to be employed in most cases. This is because criteria suggest the use of the simplex-based model for some values of $K$ , and recommend the opposite in other cases. Selection criteria deserve, though, a cautionary note in light of the non-conclusive results obtained. We suggest practitioners reinforce their conclusions visually exploring the fitted values for all mixture models before taking a final decision. One possible way to take advantage of restricted models is combining their predictions into a single one, ’à $|$ a’ committee models (Bishop, 2006). Further research should address the location and dispersion submodels related to a set of regressors through predictors, including concomitant variables into the weight parameter and incorporate the distance among pixels (Wedel, 2002; Levada et al., 2008).

Footnotes

Acknowledgments

F.O. López Quintero research was partially supported by UTFSM under PIIC No. 017/2014. J.E. Contreras-Reyes research was funded by CONICYT No. 21160618 (Res. Ex. 4128/2016). We would like to thank the editor and an anonymous reviewer for their helpful comments, suggestions and valuable discussion of this work.

References

Andrade

ACGd

(2007) Efeitos da especificação incorreta da função de ligação no modelo de regressão beta [Effects of the incorrect specification of the link function in Beta regression]. Master's thesis, Universidade de São Paulo.

Artes

Jorgensen

(2000) Longitudinal data estimating equations for dispersion models. Scandinavian Journal of Statistics , 27, 321–34.

Bagnato

Punzo

(2013) Finite mixtures of unimodal beta and gamma densities and the k-bumps algorithm. Computational Statistics , 28, 1571–97.

Barndorff-Nielsen

Jorgensen

(1991) Some parametric models on the simplex. Journal of Multivariate Analysis , 39, 106–16.

Bayes

Bazán

Garcia

(2012) A new robust regression model for proportions. Bayesian Analysis , 7, 841–66.

Bishop

(2006) Pattern Recognition and Machine Learning . New York, NY: Springer.

Burnham

Anderson

(2002) Model Selection and Multimodel Inference: A Practical Information-theoretic Approach, 2nd edition. New York, NY: Springer.

Burzynski

Janicki

Burzynski

Brookman

(2014) Preliminary findings on the use of targeted therapy with pazopanib and other agents in combination with sodium phenylbutyrate in the treatment of glioblastoma multiforme. Journal of Cancer Therapy , 5, 1423.

Carreira-Perpiñán

MÁ

(2000) Mode-finding for mixtures of Gaussian distributions. IEEE Transactions on Pattern Analysis and Machine Intelligence , 22, 1318–23.

10.

Clark

Vendt

Smith

Freymann

Kirby

Koppel

Moore

Phillips

Maffitt

Pringle

M, et. al.

(2013) The Cancer Imaging Archive (TCIA): Maintaining and operating a public information repository. Journal of Digital Imaging , 26, 1045–57.

11.

Contreras-Reyes

Cortés

(2016) Bounds on Rényi and Shannon entropies for finite mixtures of multivariate skew-normal distributions: Application to swordfish (Xi-phias gladius Linnaeus). Entropy , 18, 382.

12.

Dempster

Laird

Rubin

(1977) Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B , 21, 1–38.

13.

Greene

(2012) Econometric Analysis: International Edition: Global Edition (Pearson series in economics). Upper Saddle River, NJ, US: Pearson Education Limited.

14.

Grun

Kosmidis

Zeileis

(2012) Extended beta regression in R: Shaken, stirred, mixed, and partitioned. Journal of Statistical Software , 48, 1–25.

15.

Liu

Wang

Coombes

(2005) Applications of beta-mixture models in bioinformatics. Bioinformatics , 21, 2118–22.

16.

Jørgensen

(1987) Exponential dispersion models. Journal of the Royal Statistical Society Series B , 49, 127–62.

17.

Jørgensen

(1997) The Theory of Dispersion Models . London, UK: CRC Press.

18.

Jørgensen

Kokonendji

(2016) Discrete dispersion models and their Tweedie asymptotics. AStA Advances in Statistical Analysis , 100, 43–78.

19.

Jørgensen

Martínez

(2013) Multivariate exponential dispersion models. In Multivariate Statistics: Theory and Applications. Proceedings of the IX Tartu Conference on Multivariate Statistics & XX International Workshop on Matrices and Statistics , pages 73–98. Singapore: IEEE.

20.

Kullback

Leibler

(1951) On information and sufficiency. The Annals of Mathematical Statistics , 22, 79–86.

21.

Levada

Mascarenhas

Tannús

(2008) Pseudolikelihood equations for Potts MRF model parameter estimation on higher-order neighborhood systems. IEEE Geoscience and Remote Sensing Letters , 5, 522–26.

22.

Leijon

(2009) Beta mixture models and the application to image classification. In 16th IEEE International Conference on Image Processing (ICIP) , pages 2045–48. New Jersey, US: IEEE.

23.

Markitsis

Lai

(2010) A censored beta mixture model for the estimation of the proportion of non-differentially expressed genes. Bioinformatics , 26, 640–46.

24.

McCullagh

Nelder

(1989) Generalized Linear Models. Vol. 37. London, UK: CRC press.

25.

McLachlan

Peel

(2000) Finite Mixture Models . New York, US: John Wiley & Sons.

26.

Nair

Revathy

Tatavarti

(2008) Removal of salt-and pepper noise in images: A new decision-based algorithm. In Proceedings of the International Multi Conference of Engineers and Computer Scientists. Vol. 2. Hong Kong: Citeseer.

27.

Nikulin

(2001) Hellinger distance. Encyclopedia of Mathematics . Berlin Heidelberg and New York: Springer.

28.

Nocedal

Wright

(2006) Numerical Optimization, 2nd edition. New York, NY, US: Springer Science & Business Media.

29.

Ospina

Ferrari

(2010) Inflated beta distributions. Statistical Papers , 51, 111–26.

30.

Pawitan

(2001) In All Likelihood: Statistical Modelling and Inference Using Likelihood. Oxford, UK: Oxford University Press.

31.

Smith

Clark

Bennett

Nolan

Kirby

Wolfsberger

Moulton

Vendt

Freymann

(2016) Radiology data from the Cancer genome atlas glioblastoma multiforme (TCGA-GBM) collection. URL https://dx-doi-org.web.bisu.edu.cn/10.7937/K9/TCIA.2016.RNYFUYE9(last accessed 1 March 2016)

32.

Song

PX-K

(2007) Correlated Data Analysis: Modeling, Analytics, and Applications. New York, NY, US: Springer Science & Business Media.

33.

Vuong

(1989) Likelihood ratio tests for model selection and non-nested hypotheses. Econometrica , 57, 307–33.

34.

Wedel

(2002) Concomitant variables in finite mixture models. Statistica Neerlandica , 56, 362–75.

35.

Yao

(2010) A profile likelihood method for normal mixture with unequal variance. Journal of Statistical Planning and Inference , 140, 2089–98.

Estimation for finite mixture of simplex models: applications to biomedical data

Abstract

Keywords

1 Introduction

2 Background and related work

2.1 Dispersion models

3.1 Mixture distribution

Some simplex mixture densities.

4.1 Simulation

Table 1:

Estimation (Estim.), relative bias (RB) and MSE of simulations

Estimation of the gene expression models

4.2.1 Gene expression data

4.2.2 Glioblastoma multiforme data

Figure 2:

Fitted results in gene expression example: _______: simplex, ______: beta, ____: normal.

Totalized Voung's test for GBM images. Greatest log-likelihood results in parenthesis. Median Hellinger Distance (HD) values are appended

Three GBM images (patients) at 100 × 100 pixels size: (a), (b) and (c) refer to an original image, ZIMS and ZIMB model fit, respectively, for K = 2 . Images represented in batches (d)–(f) and (g)–(i) follow the same guidelines for K = 3 and K = 4 , correspondingly.

Continuation of Figure 3. (a)–(c) represent a GBM image (patient) for K = 5 and its estimation under ZIMS and ZIMB. This is the case when K = 6 is shown in (d)–(f).

Table 4:

Results for different ZIMS and ZIMB models of the Burzynski image example

Results for Burzynski image example for K = 2

(a): Original Burzynski image at 100 × 100 pixels size. Models with two components, (b) ZIMS estimation, (c) ZIMB estimation and (d) mixture normal estimation. Models with four components: (e) ZIMS estimation and (f) ZIMB estimation.

Results for Burzynski image example for K = 4

Burzynski image example: Models with two components: _____ zero-inflated mixture simplex estimation. ______: zero-inflated mixture beta estimation. ______: mixture normal estimation

Burzynski image example: Models with four components: ______: zero-inflated mixture simplex estimation. ______: zero-inflated mixture beta estimation. ______: mixture normal estimation.

Footnotes

Acknowledgments

References

Fitted results in gene expression example: _: simplex, : beta, : normal.

Three GBM images (patients) at 100 $\times$ 100 pixels size: (a), (b) and (c) refer to an original image, ZIMS and ZIMB model fit, respectively, for $K = 2$ . Images represented in batches (d)–(f) and (g)–(i) follow the same guidelines for $K = 3$ and $K = 4$ , correspondingly.

Burzynski image example: Models with two components: _ zero-inflated mixture simplex estimation. : zero-inflated mixture beta estimation. ____: mixture normal estimation

Burzynski image example: Models with four components: : zero-inflated mixture simplex estimation. : zero-inflated mixture beta estimation. __: mixture normal estimation.