An estimate θˆ=(pˆ,μˆ1,μˆ2,σˆ12,σˆ22) of θ based on a sample x1,…,xn is called relevant if 0<pˆ<1,σˆ12,σˆ22>0,min(x1,…,xn)≤μˆi≤max(x1,…,xn),i=1,2. 0,\\{} & \displaystyle \min (x_{1},\dots ,x_{n})\le \hat{\mu }_{i}\le \max (x_{1},\dots ,x_{n}),\hspace{1em}i=1,2.\end{array}\]]]> Let Ω denote the set of all relevant estimates of θ. The method of moments is applied to reduce Ω to a subset Ω′ of lower dimension. The first three moments of X can be expressed as (1) ν1:=E(X)=pμ1+(1−p)μ2, (2) ν2:=E(X2)=p(σ12+μ12)+(1−p)(σ22+μ22), (3) ν3:=E(X3)=p(3μ1σ12+μ13)+(1−p)(3μ2σ22+μ23), where the last equality relies on the symmetry of the component densities f1(·) and f2(·) , see Appendix A.1 for details. Following the method of moments, we replace the parameters in (1)–(3) by their estimators while equating the theoretical moments {νi} with their sample counterparts νˆ1=1n ∑i=1nxi,νˆ2=1n ∑i=1nxi2,νˆ3=1n ∑i=1nxi3. These sample moments can be highly variable so we suggest below to replace them with the more robust trimmed means, denoted {νˆk∗},k=1,2,3 , see Section 3.1 for further details. We get the following undetermined system of equations: (4) νˆ1∗=pˆμˆ1+(1−pˆ)μˆ2, (5) νˆ2∗=pˆ(σˆ12+μˆ12)+(1−pˆ)(σˆ22+μˆ22), (6) νˆ3∗=pˆ(3μˆ1σˆ12+μˆ13)+(1−pˆ)(3μˆ2σˆ22+μˆ23). The set Ω′⊆Ω consists of all relevant estimates which solve the system (4)–(6). We define the HM-estimator as an element of Ω′ with a minimum distance criteria as (7) θˆHM=arg minθˆ∈Ω′d(F(·|θˆ),Fn(·)), where the function d(F(·|θˆ),Fn(·)) measures the distance between the fitted model distribution F(·|θˆ) and the empirical distribution Fn(·) of the sample. In this paper we use an L2 -type measure given by d(F(·|θ),Fn(·))=1n ∑i=1n(F(xi|θ)−Fn(xi))2=1n ∑i=1n(F(x(i)|θ)−i/n)2, where x(1),…,x(n) is the ordered sample. It should be noted that this choice of distance is mainly due to its computational simplicity and that a number of different measures can be considered (see [26], chap. 4)

Next we describe how the HM-estimator is obtained in practice, via a reformulation of definition (7) that is more useful for computation.

In this subsection we describe how the HM-estimator (7) can be obtained in practice. To get a representation of the solutions of the system (4)–(6), we reparametrize the problem by introducing the proportion rˆ , defined by (8) σˆ22=rˆσˆ12. Equations (4), (5), and (8) can be used to eliminate μˆ2 , σˆ12 , σˆ22 in equation (6), and as a result we obtain (9) β3μˆ13+β2μˆ12+β1μˆ1+β0=0, where the coefficients are functions of pˆ and rˆ such that β0=−νˆ3∗+3νˆ1∗νˆ2∗rˆpˆ+rˆ−pˆrˆ+(νˆ1∗)33pˆ−2(pˆ+rˆ−pˆrˆ)(1−pˆ)2(pˆ+rˆ−pˆrˆ),β1=3νˆ2∗rˆ−pˆ+rˆ−pˆrˆ(1−pˆ)(pˆ+rˆ−pˆrˆ)−(νˆ1∗)23pˆ(2pˆ−2(pˆ+rˆ−pˆrˆ)+1)(1−pˆ)2(pˆ+rˆ−pˆrˆ),β2=3νˆ1∗pˆ(2pˆ−pˆ2+pˆ2rˆ−rˆ)(1−pˆ)2(pˆ+rˆ−pˆrˆ),β3=−pˆ(2pˆ−pˆ2+pˆ2rˆ−rˆ)(1−pˆ)2(pˆ+rˆ−pˆrˆ). Furthermore, by combining (4), (5), and (8), we get that the estimated parameters μˆ2 and σˆ12 are obtained as (10) μˆ2=νˆ1∗−pˆμˆ11−pˆ, (11) σˆ12=νˆ2∗−pˆμˆ12−(1−pˆ)μˆ22pˆ+(1−pˆ)rˆ. If pˆ and rˆ are given, we see that (9) is a cubic equation for μˆ1 and so the reparameterized system has at most three solutions that correspond to relevant estimates. Define M to be the set of all pairs (pˆ,rˆ) for which at least one relevant estimate exists, and let T(pˆ,rˆ) contain all relevant estimates corresponding to the pair (pˆ,rˆ)∈M . From the definitions of Ω′ , M, and T(pˆ,rˆ) it follows that (12) minθˆ∈Ω′d(F(·|θˆ),Fn(·))=min(pˆ,rˆ)∈Mg(pˆ,rˆ), where the function g(pˆ,rˆ)=minθˆ∈T(pˆ,rˆ)d(F(·|θˆ),Fn(·)),(pˆ,rˆ)∈M, is straightforward to compute since T(pˆ,rˆ) contains at most three elements and these are solutions of the polynomial equations (8)–(11). Equation (12) reformulates the problem of deriving the HM-estimator (7) to a minimization problem for the bivariate function g(pˆ,rˆ),(pˆ,rˆ)∈M . Let (pˆHM,rˆHM) denote the point that minimizes the function g(pˆ,rˆ),g(pˆ,rˆ)∈M . Then the estimator is given as θˆHM=arg minθˆ∈T(pˆHM,rˆHM)d(F(·|θˆ),Fn(·)).

The minimization of g(pˆ,rˆ) can be obtained using a numerical optimization algorithm. Here we use the simplex algorithm in [21], which is implemented in the optim routine in R[23]. The starting value is found as the minimizer of g(pˆ,rˆ) over a finite grid of values. A schematic description of how the HM-estimator is computed is given in Figure 1.

The mixture density f(·) is typically used to model a data set x1,…,xn where n1 of the values are observations from component f1(·) and the remaining n2=n−n1 are observations from f2(·) . For such a sample, we can introduce a 0–1 vector z=(z1,…,zn) that correctly assigns each observation to either f1(·) (1’s) or f2(·) (0’s). This vector defines a true partition of the observations with respect to the density components. Usually the components represent distinct subpopulations.

The true partition defined by z=(z1,…,zn) is unobservable but can be estimated with the posterior membership probabilites, also known as the responsibilites, denoted by zˆ=(zˆ1,…,zˆn) , where (13) zˆi=pˆf1(xi|θˆ1)f(xi|θˆ)=pˆf1(xi|θˆ1)pˆfˆ1(xi|θˆ1)+(1−pˆ)fˆ2(xi|θˆ2),i=1,…,n. The value zˆi is in the interval [0,1] and estimates the probability that xi is an observation from component f1(·) . The vector zˆ=(zˆ1,…,zˆn) defines a so-called soft partition of the data x1,…,xn which serves as an approximation of the true partition given by z=(z1,…,zn) . The responsibilities in (13) can be obtained for any estimator θˆ of θ, e.g. the proposed HM-estimator or the maximum likelihood (ML) estimator used in the numerical studies in Sections 3 and 4.

In the simulations, we restrict ourselves to the case where the component densities f1(·) and f2(·) belong to the same family of distributions. The estimators are calculated under the assumption that f1(·) and f2(·) are normal densities, which is a common assumption in practice. The data are generated from normal mixtures, for which the assumption is true, and also from mixtures of Laplace, logistic and contaminated Gaussian distributions (for details, see Appendix A.2). For the contaminated Gaussian distribution, we set the larger prior probability to α=0.9 and the variance proportion parameter to η=16 . The experiment thus includes both modest and large departures from the normal mixture assumption, allowing us to analyze the robustness of the methods with respect to imprecise model specifications.

Besides varying the family of the component densities, we consider six configurations of the parameter vector θ which correspond to a variety in shape of the mixture distributions. In addition to these configurations a negative control with a non-mixture distribution was added. The values are given in Table 1 and displayed graphically in Figure 2. Three sample sizes are considered; n=50,100 , and 500.

The mixture data were generated as follows: first we simulated the (true) partition vector z=(z1,…,zn) from the 0–1 variable Z with P(Z=1)=p and P(Z=0)=1−p . Then x1(1),…,xn(1) and x1(2),…,xn(2) were simulated from the component densities f1(·) and f2(·) , respectively. Finally a sample x1,…,xn from the mixture density f(·) was obtained as xi=zixi(1)+(1−zi)xi(2),i=1,…,n. We generated N=500 data sets for each considered scenario (mixture family, parameter configuration, and sample size), and for each scenario we obtained 500 independent realizations θˆ(1),…,θˆ(500) of an estimator θˆ , which were used for statistical evaluation of the performance of the considered methods in clustering and in estimation of the mixing proportion.

The hybrid estimator θˆHM was obtained as described in the previous section. The starting value for the simplex method was found as the minimizer of g(pˆ,rˆ) over a two-dimensional grid constructed from 10 values of pˆ and 200 values of rˆ . The values of pˆ and rˆ in the grid were evenly distributed in the intervals (0,1) and [0,20] , respectively. We used trimmed versions of the sample moments νˆk,k=1,2,3 . The 2.5% smallest and 2.5% largest of values in x1k,…,xnk were removed and the mean νˆk∗ of the resulting trimmed sample was used as the estimator of the kth moment νk .

The maximum likelihood estimator θˆML was calculated with the EM-algorithm using the mixtools package in R [1]. The EM-algorithm converges at a local maximum of the likelihood function that depends on its starting value. We chose ten starting values randomly as described in [14], and the estimator θˆML was taken as the maximizer of the likelihood among the corresponding points of convergence.

For a simulated dataset x1,…,xn the true partition z=(z1,…,zn) was known. The true partition was approximated by the soft partitions zˆHM and zˆML , calculated from the corresponding estimates θˆHM and θˆML , respectively, using Equation (13). To quantify the accuracy of an approximate partition, we used the Fuzzy Adjusted Rand Index (FARI) proposed in [3]. The FARI for z and its approximation zˆ – written as FARI( zˆ,z ) – is a number in the interval [−1,1] measuring their closeness; the higher the index the better is the approximation zˆ . A brief description of this index is given in Appendix A.3.

Let ΔFARI=FARI(z,zˆHM)−FARI(z,zˆML) denote the difference between the indices. Note that a positive difference ΔFARI>0 0$]]> implies that the partition obtained via the HM-estimator was more accurate than the partition obtained via the ML-estimator.

To determine if there was a significant difference between the methods’ clustering performance, we applied the t-test and the sign-test to the pairwise differences ΔFARI(1),…,ΔFARI(500) for the 500 simulated samples. We also made a comparison of the methods given that a difference in FARI under a certain threshold was considered as negligible, which was achieved by applying the sign-test to the differences that satisfied |ΔFARI(i)|>0.1 0.1$]]>.

The considered scenarios corresponded to problems that were more or less difficult with respect to clustering and as part of our evaluations we quantified these difficulties. Here zˆopt denotes the optimal partition obtained when the true component densities and parameter values in (13) were used. The index FARIopt=FARI(z,zˆopt) corresponded to the clustering performance obtained under correct model assumptions and a perfect estimator of θ. For each scenario, we used the mean of FARIopt(1),…,FARIopt(500) for the 500 simulated samples to measure the difficulty of the problem and as a reference value for the corresponding FARI obtained by the HM- and ML-estimators.

We compared the methods in terms of their accuracy for estimating the mixing proportion p. Details on how we defined the point estimators of p are given in Appendix A.4.

The following standard characteristics for evaluating an estimator pˆ of p based on N simulations were used: (14) mean=1N ∑i=1Npˆ(i),biasˆ=1N ∑i=1N(pˆ(i)−p),MSEˆ=1N ∑i=1N(pˆ(i)−p)2. These characteristics were calculated for the simulated estimates {pˆHM(i)}i=1N and {pˆML(i)}i=1N obtained from the HM- and ML-estimators, respectively.

To determine if there was a significant difference in efficiency between the methods, we applied the t-test to the difference in estimated mean squared error ΔMSE=MSEˆML−MSEˆHM , where the subscripts HM and ML refer to the methods used in (14). Note that a positive value of ΔMSE suggested that the hybrid method was more efficient as an estimator of p.

This section includes a detailed treatment of the results for sample size n=50 . The results for sample sizes n=100 and n=500 led to similar conclusions, and are given in Appendix A.5.

The HM- and ML-estimators were evaluated mainly by their ability to cluster the samples in agreement with the true partition vector z and their ability to estimate the proportion parameter p. The corresponding estimators for the difference in mean μ2−μ1 were compared in a similar way as for the mixing proportion p. This was done for the case when the data were generated from a normal mixture distribution (i.e. under correct model assumptions), a logistic mixture distribution (modest violation of model assumptions), a Laplace mixture distribution and contaminated Gaussian distribution (serious violation of model assumptions). For each case six values of the parameter vector were evaluated, Table 1 and Figure 2.

In [13] a microarray experiment on human mRNA samples for measuring gene expression levels in two subtypes of acute leukemia is described, namely acute lymphoblastic leukemia (ALL) and acute myeloid leukemia (AML). The experiment contained n=72 samples, of which 47 were of type ALL and 25 were of type AML, and the expression levels of 6,817 genes were observed. In this study we used the preprocessed and filtered version of this data considered in [8], which contains the expression values of 3,571 genes.

We applied a supervised procedure to get a subset of the 3,571 genes which were expressed differently with respect to the ALL/AML grouping. For each gene i we calculated the normal mixture clustering zˆ(i) in (13) with parameter estimates obtained under known group memberships (i.e. the sample means and variances). Then we used (the measure) FARI(z,zˆ(i)) , where z is the true ALL/AML grouping expressed as a 0–1 vector, to quantify the extent to which the mean expression of gene i differs between the groups. The 342 genes that met the criterion FARI(z,zˆ(i))≥0.1 were regarded as truly differently expressed genes and included in our test set. Alternatively, we could have used the two-sample t-test to make this selection.

We applied the HM and ML clustering methods to each of the 342 test variables to compare their ability to divide the 72 cancer samples into the ALL and AML groups. The analysis was carried out as described in Section 3.2.

Independent of method the overall performance was rather poor; most of the clusterings had a FARI between 0 and 0.4, Figure 7. This implies that it is hard to cluster the samples accurately based on single genes.

The observed differences between the HM and ML methods were all significant and in favor of the hybrid method. For 213 of the 342 test genes the HM method clustered the samples more accurately than the ML method (i.e.  propHM=0.623 , p-value < 10−5 ). The mean difference in FARI (i.e.  Δ‾FARI ) was 0.020 (p-value <10−5 ). Moreover, of the 38 genes for which there was a considerable difference between the methods (i.e. an absolute difference larger than 0.1), the HM clustering was superior over the ML clustering in 32 of the cases (i.e.  propCHM=0.842 , p-value <10−4 ). For notation, see Section 3.4.1.