Mixture models for gene expression experiments with two species

Abstract

Cross-species research in drug development is novel and challenging. A bivariate mixture model utilizing information across two species was proposed to solve the fundamental problem of identifying differentially expressed genes in microarray experiments in order to potentially improve the understanding of translation between preclinical and clinical studies for drug development. The proposed approach models the joint distribution of treatment effects estimated from independent linear models. The mixture model posits up to nine components, four of which include groups in which genes are differentially expressed in both species. A comprehensive simulation to evaluate the model performance and one application on a real world data set, a mouse and human type II diabetes experiment, suggest that the proposed model, though highly structured, can handle various configurations of differential gene expression and is practically useful on identifying differentially expressed genes, especially when the magnitude of differential expression due to different treatment intervention is weak. In the mouse and human application, the proposed mixture model was able to eliminate unimportant genes and identify a list of genes that were differentially expressed in both species and could be potential gene targets for drug development.

Introduction

Background

Pharmaceutical medicine is an industry with huge up-front investment for rewards that may or may not come years later. A complete drug development process, including drug discovery, preclinical research (on animals) and clinical trials (on humans), is lengthy, expensive, and risky. Determined by the US Food and Drug Administration (FDA) [1], the average total cost per drug development is about $1.9 billion. The typical development time is 10 to 15 years. The overall attrition rate of a drug compound from first-in-man to registration is approximately 80%–90% [2, 3].

FDA [1] calls the preclinical and clinical research together as the ‘critical path’ development phase, where most investment required for a successful drug launch occurs. Currently, this development phase is inherently inefficient. The goal of preclinical research is to assess how a drug is absorbed, distributed, metabolized, and excreted in animals, and to use the findings to determine potential human outcomes before starting clinical trials. Yet the rate of success after a drug candidate enters Phase I is undesirably low. As mentioned in FDA [1] and Kola and Landis [3], animal models with poor clinical relevance may be accountable for this perplexity. Hence, improving translation between two species to increase the predictive power of animal models to human studies is of tremendous value to drug discovery and development.

Homology and multiple species gene expression analysis in drug development

Microarrays are tools for gene expression analysis and can be potentially useful for investigating the mechanism of drug activities that translates across species. The utility of microarray information in the drug development process is reviewed by Braxton and Bedilion [4] who embraced the idea that gene expression analysis can be a surrogate marker for the interaction between compounds and cells and should yield information about efficacy. Debouck and Goodfellow [5] believed that microarrays can be used to generate clues to patterns of gene function that can help improve the efficiency of drug development.

As stated before, one key challenge of drug development is to successfully translate the results of preclinical findings in animal models to human beings in the clinic. Pre-clinical experiments assume that the effect of the drug tested on animals is comparable to that on humans, which can only be true if a functional equivalent of the human drug target exists in the experimental species. Orthology [6–8] is a strong indication of functional conservation and therefore provides the best functional annotation of experimentally undetermined genes across species. Holbrook and Sanseau [9] remarked that the use of orthologs has the potential to improve the understanding of biological differences between species (animals and humans).

Many of the successful applications of cross-species microarray gene expression analysis involve orthology [10–13]. Additionally, over the past decade, researchers have tried to use orthology and gene expression data to do cross-species comparison in order to understand how genes interact to perform particular biological processes [14, 15]. These studies support the idea that orthologs could be a useful tool for researchers to link experiments between species in drug development. Note that orthologous relationships can be one-to-one, one-to-many, or many-to-many [8].

The rest of the paper is structured as follows: Section ‘Joint modeling across species’ describes the proposed bivariate mixture model across species. Section ‘Simulation’ describes a simulation study undertaken to investigate the effects of different experimental designs on the power to detect important genes and on misclassification rates. Section ‘Application: the mouse and human type II diabetes experiment’ illustrates the methodology using an application to data collected in a mouse/human experiment. Section ‘Concluding remarks’ concludes.

Joint modeling across species

Let $X_{a_{\mathit{\text{ij}}}}$ and $X_{h_{\mathit{\text{il}}}}$ denote gene expression measurements from the ith orthologous gene pair for the jth animal and the lth human. The following independent linear models describe the association between gene expression and treatment:

$$\begin{array}{c}{X}_{a_{\mathit{\text{ij}}}}={\beta }_{0a_i}+{\beta }_{1a_i}{T}_{a_j}+e_{a_{\mathit{\text{ij}}}},\end{array}$$

(1)

$$\begin{array}{c}{X}_{h_{\mathit{\text{il}}}}={\beta }_{0h_i}+{\beta }_{1h_i}{T}_{h_l}+e_{h_{\mathit{\text{il}}}},\end{array}$$

(2)

where $T_{a_j}$ and $T_{h_l}$ are {0,1} treatment indicators, and $e_{a_{\mathit{\text{ij}}}}$ and $e_{h_{\mathit{\text{il}}}}$ are independent $N(0,{\sigma }_a^2)$ and $N(0,{\sigma }_h^2)$ random variables. ${\sigma }_a^2$ and ${\sigma }_h^2$ are variances for $e_{a_{\mathit{\text{ij}}}}$ and $e_{h_{\mathit{\text{il}}}}$, respectively. In drug development, the animal research and human experiments are conducted independently - one’s results do not affect the other’s. However, the treatment effects are expected to have some kind of association between the two species. This results in our choice of using two independent models for the two species to capture the effects of treatment on gene expression.

A nine-component bivariate mixture model for two species experiments

${\beta }_{1a_i}$ and ${\beta }_{1h_i}$ quantify the differential expression of the ith orthologous animal and human genes due to a treatment intervention. A given gene can be classified as non-differentially expressed (NDE) - showing no signs of treatment effects, positively differentially expressed (pDE) - showing positive treatment effects, or negatively differentially expressed (nDE) - showing negative treatment effects. Therefore, for a human and animal gene pair, there are nine possibilities for categorizing this pair of genes. Further, dependency is assumed between differentially expressed orthologs, i.e., existence of association posited only for ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ in categories (1, 2, 3, 4) and zero correlation presumed for ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ in categories (0, 5, 6, 7, 8). Table 1 illustrates the nine possible categories of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$. ${({\mu }_{{\beta }_{1a_i}},{\mu }_{{\beta }_{1h_i}})}^T$ is the vector of population means of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ under each category.

Table 1 Possible categories of ${\mathbf{(}{\mathit{\beta }}_{\mathbf{1}{\mathit{a}}_{\mathit{i}}}\mathbf{,}{\mathit{\beta }}_{\mathbf{1}{\mathit{h}}_{\mathit{i}}}\mathbf{)}}^{\mathit{T}}$

In consequence of these possible patterns of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$, mixture models [16, 17] are adopted to deal with the correlation and distribution of each subgroup of genes across species. An additional advantage of mixture models is that, after prior weights for the components are specified, estimates of the posterior probabilities of population membership can be formed for each observation to give a probabilistic clustering. As a result, the pooling of information for genes across species can be exploited to better understand the underlying relationship between the treatment intervention for both species.

Tailoring the mixture model to two-species experiments with restrictions on the parameters made according to Table 1 and assuming that the treatment effects for non-differentially expressed genes are deterministically zero, i.e., ${({\beta }_{1a_i},{\beta }_{1h_i})}^T={(0,0)}^T$, the following bivariate normal mixture model is adopted as the prior distribution of the vector ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$:

$$\begin{array}{ll}\left(\begin{array}{c}{\beta }_{1a_i}\\ {\beta }_{1h_i}\end{array}\right)& \sim {\pi }_{0}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a0}\\ {\mu }_{h0}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a0}^2& {\rho }_0{\eta }_{a0}{\eta }_{h0}\\ {\rho }_0{\eta }_{a0}{\eta }_{h0}& {\eta }_{h0}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{1}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a1}\\ {\mu }_{h1}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a1}^2& {\rho }_1{\eta }_{a1}{\eta }_{h1}\\ {\rho }_1{\eta }_{a1}{\eta }_{h1}& {\eta }_{h1}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{2}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a2}\\ {\mu }_{h2}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a2}^2& {\rho }_2{\eta }_{a2}{\eta }_{h2}\\ {\rho }_2{\eta }_{a2}{\eta }_{h2}& {\eta }_{h2}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{3}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a3}\\ {\mu }_{h3}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a3}^2& {\rho }_3{\eta }_{a3}{\eta }_{h3}\\ {\rho }_3{\eta }_{a3}{\eta }_{h3}& {\eta }_{h3}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{4}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a4}\\ {\mu }_{h4}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a4}^2& {\rho }_4{\eta }_{a4}{\eta }_{h4}\\ {\rho }_4{\eta }_{a4}{\eta }_{h4}& {\eta }_{h4}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{5}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a5}\\ {\mu }_{h5}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a5}^2& {\rho }_5{\eta }_{a5}{\eta }_{h5}\\ {\rho }_5{\eta }_{a5}{\eta }_{h5}& {\eta }_{h5}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{6}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a6}\\ {\mu }_{h6}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a6}^2& {\rho }_6{\eta }_{a6}{\eta }_{h6}\\ {\rho }_6{\eta }_{a6}{\eta }_{h6}& {\eta }_{h6}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{7}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a7}\\ {\mu }_{h7}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a7}^2& {\rho }_7{\eta }_{a7}{\eta }_{h7}\\ {\rho }_7{\eta }_{a7}{\eta }_{h7}& {\eta }_{h7}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{8}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a8}\\ {\mu }_{h8}\end{array}\right),& \left(\begin{array}{cc}{\eta }_{a8}^2& {\rho }_8{\eta }_{a8}{\eta }_{h8}\\ {\rho }_8{\eta }_{a8}{\eta }_{h8}& {\eta }_{h8}^2\end{array}\right)\end{array}\right),\end{array}$$

(3)

where π _k is the probability that an observation belongs to the kth component, with ${\sum }_{k=0}^8{\pi }_k=1\text{and}{\pi }_k\ge 0.$ The following restriction of the parameter space is imposed: μ_a0 = 0, μ_h0 = 0, μ_a1 ≥ 0, μ_h1 ≥ 0, μ_a2 ≤ 0, μ_h2 ≤ 0, μ_a3 ≥ 0, μ_h3 ≤ 0, μ_a4 ≤ 0, μ_h4 ≥ 0, μ_a5 = 0, μ_h5 ≥ 0, μ_a6 = 0, μ_h6 ≤ 0, μ_a7 ≥ 0, μ_h7 = 0, μ_a8 ≤ 0, μ_h8 = 0, η_a0 = 0, η_h0 = 0, η_a5 = 0, η_a6 = 0, η_h7 = 0, η_h8 = 0, ρ₀ = 0, ρ₅ = 0, ρ₆ = 0, ρ₇ = 0, and ρ₈ = 0.

According to the theory of least squares, the marginal distribution of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$, the parameter estimates, has means equal to the prior means of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ and variances involving contributions from the prior distribution of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ and the conditional distribution of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ given ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$. The marginal distribution of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ is as follows:

$$\begin{array}{ll}\left(\begin{array}{c}{\hat{\beta }}_{1a_i}\\ {\hat{\beta }}_{1h_i}\end{array}\right)& \sim {\pi }_{0}N\left(\begin{array}{cc}\left(\begin{array}{c}0\\ 0\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a0}^2& 0\\ 0& {\sigma }_{h0}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{1}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a1}\\ {\mu }_{h1}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a1}^2& {\rho }_1{\sigma }_{a1}{\sigma }_{h1}\\ {\rho }_1{\sigma }_{a1}{\sigma }_{h1}& {\sigma }_{h1}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{2}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a2}\\ {\mu }_{h2}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a2}^2& {\rho }_2{\sigma }_{a2}{\sigma }_{h2}\\ {\rho }_2{\sigma }_{a2}{\sigma }_{h2}& {\sigma }_{h2}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{3}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a3}\\ {\mu }_{h3}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a3}^2& {\rho }_3{\sigma }_{a3}{\sigma }_{h3}\\ {\rho }_3{\sigma }_{a3}{\sigma }_{h3}& {\sigma }_{h3}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{4}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a4}\\ {\mu }_{h4}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a4}^2& {\rho }_4{\sigma }_{a4}{\sigma }_{h4}\\ {\rho }_4{\sigma }_{a4}{\sigma }_{h4}& {\sigma }_{h4}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{5}N\left(\begin{array}{cc}\left(\begin{array}{c}0\\ {\mu }_{h5}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a5}^2& 0\\ 0& {\sigma }_{h5}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{6}N\left(\begin{array}{cc}\left(\begin{array}{c}0\\ {\mu }_{h6}\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a6}^2& 0\\ 0& {\sigma }_{h6}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{7}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a7}\\ 0\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a7}^2& 0\\ 0& {\sigma }_{h7}^2\end{array}\right)\end{array}\right)\\ +{\pi }_{8}N\left(\begin{array}{cc}\left(\begin{array}{c}{\mu }_{a8}\\ 0\end{array}\right),& \left(\begin{array}{cc}{\sigma }_{a8}^2& 0\\ 0& {\sigma }_{h8}^2\end{array}\right)\end{array}\right),\end{array}$$

(4)

An EM algorithm is developed to accomplish the nontrivial likelihood maximization, along with methodology for handling singular covariance matrices that arise during the implementation of the algorithm. (See the Appendix for details). Gene membership is determined according to the maximum posterior probability that an observation ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ comes from the kth component of the mixture.

Simulation

The following Monte Carlo simulation studies investigated the performance of the proposed mixture model using information across two species in comparison to the traditional microarray method using just one-species information when identifying genes associated with treatment stimulus under several different scenarios.

Several factors influence the sampling properties of the estimated treatment effects ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ using the mixture model were considered in the simulation:

•  Replicates (number of arrays) per treatment for each species: n _a and n _h for animals and humans, respectively.

•  Number of orthologous genes in each category: n _k , k = 0,…,8, the kth category.

•  Array noise: $e_{a_{\mathit{\text{ij}}}}$ and $e_{h_{\mathit{\text{ij}}}}$ in (1) and (2). Also recall that, by assumption, $e_{a_{\mathit{\text{ij}}}}$ and $e_{h_{\mathit{\text{il}}}}$ are independent $N(0,{\sigma }_a^2)$ and $N(0,{\sigma }_h^2)$ random variables.

•  Parameters in (3) by which the sampling distribution of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ is determined.

With so many variables, it is impractical to study the sampling properties of the fitted model without fixing some variables. That is, an experimental design for the simulation in which these factors are completely considered is not feasible. The simulation study instead focused on three aspects. First, although high-density microarrays provide useful genome-wide data, they are often associated with a substantial amount of experimental noise that could affect the performance of the analysis. Hence, it is of interest to investigate how the array noise would affect the model efficiency on gene identification across species.

Second, the sample size of cross-species experiment is likely to be different, and may be one of the deciding factors of the power associated with the modeling approach. In particular, the efficiency of gene identification, whether the proposed model gains power over one-species experiments through pooling information across species, should be examined carefully, especially when the sample size of the experiments is small.

Third, over-fitting may be of concern. The proposed mixture model is, by its nature, highly structured and data driven. If the data are not driven by all nine categories as the model suggests, will the mixture model fail? Is the proposed model flexible enough to handle different types of data structure? To examine the model performance systematically and to test if the proposed model will fail when too many components are used to fit the data where there are actually fewer clusters, two types of data were generated: all nine categories non-empty (case I) and some of the nine components empty (case II).

In simulation studies case I and case II, two methods of gene identification were implemented: the proposed mixture model, utilizing information across two different species, and the traditional t statistics adjusting for multiple comparisons based on single-species data. Five hundred data sets were generated for each different scenario under each simulation study.

Parameter determination and data generation

Theoretically, the number of genes in category 0 (non-differentially expressed in both species) should dominate others, and every other category may comprise some genes. Orthology information from HomoloGene of the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/homologene) and Mouse Genome Informatics (MGI) of Jackson Laboratory [18] and the practical experience gained from analysis of two-species gene expression experiments at GlasoSmithKline (GSK) was used as reference to determine a reasonable number of genes in each category.

The vectors of the number of genes in each category (n₀,n₁,…,n₈)^T determined for simulation studies case I and case II are categories (6,000,30, 30,30,30,100,100,100,100)^T and categories (6,000,30, 30,0,0,100,0,100,0)^T, respectively. For experiments across species, sample sizes may differ. Considering that this proposed bivariate method could benefit from pooling information across species, especially when the sample size is small, and the practical situation, two scenarios were implemented: the number of replicates per treatment for each species is equal and small, and the number of replicates per treatment for animals is greater than humans. In addition, to evaluate how robust the proposed method is against array noise, two situations were considered: the two experiments are equally noisy and the human data are noisier than the animal’s. Furthermore, values of parameters in (3) for the sampling distribution of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ were predetermined. Variances of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ in each component were assumed to be the same: ${\eta }_{a1}^2={\eta }_{a2}^2={\eta }_{a3}^2={\eta }_{a4}^2={\eta }_{a7}^2={\eta }_{a8}^2={\eta }_{h1}^2={\eta }_{h2}^2={\eta }_{h3}^2={\eta }_{h4}^2={\eta }_{h5}^2={\eta }_{h6}^2=0.25$. The correlation between ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ in categories (1,2,3,4) was assumed to be 0.9, i.e., ρ₁ = ρ₂ = ρ₃ = ρ₄ = 0.9. Nonzero component means (μ _ak ,μ _hk )^T, k = 0,…,8, were determined so that |μ/η| = 0.5 or 1.5. The combination of these parameters resulted in eight different scenarios for each case as presented in Table 2. Note that $|{\mu }_{{\beta }_{1a_i}}|$ and $|{\mu }_{{\beta }_{1h_i}}|$ represent the absolute value of the mean vector of ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$.

Table 2 Combination of parameters for simulation studies case I and case II

After generating ${({\beta }_{1a_i},{\beta }_{1h_i})}^T$ accordingly, the next step is to simulate the two species gene expression data X _a and X _h based on linear models (1) and (2). Note that ${\beta }_{0a_i}$ and ${\beta }_{0h_i}$ are independent N(8,1) random variables for differentially expressed genes and deterministically 0 for non-differentially expressed genes.

Simulation results

It is of interest to compare how effective the mixture model is on gene identification using information across two species with the conventional one-species approach. The conventional two-sample t test for gene selection was performed using just single species data (animals or humans) and a multiplicity adjustment was made according to the procedure proposed by Benjamini and Hochberg [19].

The results are presented in Table 3. The first section of Table 3, categories (1, 2, 3, 4, 5, 6), manifests the number of genes classified into categories (1, 2, 3, 4, 5, 6) using the the mixture model (Mixture) and the number of genes selected using human data only (Human only), with the corresponding nominal FDR controlled at FDR _I . For each simulated data set, FDR _I was calculated as (number of genes that are erroneously classified into categories (1, 2, 3, 4, 5, 6))/(total number of genes classified into categories (1, 2, 3, 4, 5, 6)). This was to ensure a fair comparison between the mixture model and the conventional one-species method. FDR _II was calculated in the same fashion. Avg FDR _I and Avg FDR _II are simply the averaged values of FDR _I and FDR _II across the 500 simulated data sets. When the estimated nominal FDR = 0, i.e., the mixture model did not falsely categorize any genes, the nominal FDR for the single-species method was controlled at 0.0001. The second section of Table 3, categories (1,2,3,4,7,8), represents the number of genes classified into categories (1,2,3,4,7,8) using the mixture model (Mixture) and number of genes selected using animal data alone (Animal only). Beneath each set of eight simulation cases is Tukey’s Honestly Significant Difference (HSD) for a familywise error rate of 0.05 for the results obtained using the proposed mixture model. Tukey’s HSD was calculated as q_0.05(8,3,992) × (MS(Error) /500)^1/2, since there were eight simulation cases (500 simulated data sets in each situation) in each simulation study (case I and case II) and the error degree of freedom was 3,992. MS(Error) denotes the error mean square (= Error sum of squares/Error degree of freedom, see Table 4). q_0.05(8,3,992) = 4.29.

Table 3 The number of genes selected based on (a) bivariate mixture model, (b) conventional one-species approach

Table 4 ANOVA table to quantify variability

The effect of array noise on the mixture model can be easily seen from the column of Avg FDR _I and by comparing the results of sim1 vs. sim3, sim5 vs. sim7, sim9 vs. sim11, and sim13 vs. sim15. These are the cases with smaller means and variances for humans which change from 0.1 to 0.3 for each pair of comparisons. The differences between the observed FDRs for these four groups were at least 0.2 and 0.3 for case I and case II, respectively. In contrast, when means were larger, changes of variances did not seem to affect the results in the sense that the corresponding observed FDRs had barely changed while the variances of human increased. The observed FDRs for animals were not as sensitive to the array noise on humans as the observed FDRs for humans.

Under both simulation studies case I and case II, increasing the number of replicates in the animal experiment helped the gene identification for animals: more animal genes were identified for sim5 to sim8 than for sim1 to sim4, and for sim13 to sim16 than for sim9 to sim12. The corresponding FDRs were also lower for sim5 to sim8 and sim13 to sim16. In contrast, increasing the number of replicates in the animal experiment did not significantly improve the results of gene identification for humans.

Table 4 is the summary of a three-way analysis of variance (ANOVA). Three factors, each with two levels, were used in the analysis: replicates ((n _a ,n _h )^T = (10,10)^T or (100,10)^T), mean magnitude (${(|{\mu }_{{\beta }_{1a_i}}|,|{\mu }_{{\beta }_{1h_i}}|)}^T={(0.25,0.25)}^T$ or (0.75,0.75)^T), and array noise (${({\sigma }_a^2,{\sigma }_h^2)}^T={(0.1,0.1)}^T$ or (0.1,0.3)^T). This analysis, performed independently, quantifies the variability among the results (gene counts and observed FDRs) obtained using the proposed mixture model in Table 3 for the 16 different simulated situations under simulation studies case I and case II.

Throughout the 16 simulation cases, with nominal FDR controlled at FDR _I , the bivariate mixture model outperformed the single-species method for human gene identification by always recognizing more genes with lower observed FDRs. For the animal part, the mixture model performed at least as well as the single-species method on gene selection by identifying at least as many genes. Notice that selecting genes related to humans (categories (1, 2, 3, 4, 5, 6)) seemed to be associated with higher false discovery rate than selecting genes related to animals (categories (1, 2, 3, 4, 7, 8)). Furthermore, the observed FDRs were lower for case I than for case II, regardless of the type of genes interested (differentially expressed for humans or animals).

The comprehensive simulation study suggested that the proposed model, though highly structured, offered advantages over single-species analyses, especially when the magnitude of differential expression due to different treatment intervention was weak.

Application: the mouse and human type II diabetes experiment

Background introduction

A systems biology study was completed by GlaxoSmithKline (GSK) to study the efficacy of type II diabetes drugs in both preclinical (mice) and clinical (humans) experiments. The mouse and human data were collected and preprocessed using Affymetrix MAS 5.0 at the probe set level on Affymetrix MOE430A array and Human Genome U133 Plus 2.0 array, respectively. For mice, the total number of probe sets was 22,690. Mice were fed a diet enriched in fat (58% kcal from fat) for 8 weeks prior to treatment. Most of the mice of this susceptible strain developed obesity and mild hyperglycemia and hyperinsulinemia. Control mice on an 11% low fat diet remained normal. The mouse treatment arm consisted of a diabetes drug at multiple dose levels with vehicle controls over a 2-week period. The study was a full factorial design, where 40 animals in high fat diet and 40 animals in low fat diet were randomized to receive either placebo, or different dosages of the type II diabetes drug (low, medium or high). Note that the results for mice were measured at one time point, the end of the study. For simplicity, only integration results on animals treated with placebo (10 mice) or high dose of the drug (9 mice) were demonstrated.

On the other hand, there were 54,676 probe sets in the human data set and all 59 subjects were type II diabetes patients. The gene expression measurements for human subjects were collected twice during the experiment, one at baseline before the treatment started (week 0) and the other one at the completion of the study (week 8). Among the 59 subjects, 14 only had data for only one time point and hence were not included in further data analysis. Subjects in the clinical trial were treated by either placebo, or three other type II diabetes drugs, including the one given to the mice. Titrated dosing was implemented to ensure that each person received his/her dose based on his/her body profile. Human subjects treated with placebo (number of subjects = 11) and the same type II diabetes drug given to mice (number of subjects = 13) were used in this analysis. Besides being measured at two time points, other information contained in the human data and used in the data analysis included: prior therapy (four categorical levels) and concomitant medication status (two categorical levels).

Data for both species were logarithmic transformed (base equal to 2). Methods for combining information from multiple probe sets that were not identical have been discussed in many publications [20, 21]. A gene-level transcript value in order to pair the mouse and human genes through orthology was obtained by averaging probe sets across a gene. This resulted in 13,483 and 20,252 genes for mouse and human, respectively. Missing values were replaced with array means in both data sets.

The mouse and human orthology information, MGI release 4.32, was used to map the mouse and human orthologous genes. Human and mouse orthologs (17,834) were included in release 4.32. Combining the information among MGI release 4.32 and mouse and human genes from the GSK data gave a total of 11,922 orthologs.

In addition to the measurements of gene expression for both mice and humans, different efficacy endpoints, such as blood glucose, insulin, hemoglobin A1c (HbA1c), and others were also measured during the preclinical and clinical experiments by GSK in order to evaluate the effect of treatment intervention on both species.

The purpose of this data analysis was to evaluate the capability of the proposed mixture model, an approach utilizing the information across two species, to identify genes that may help scientists reveal the biological similarity between two species (ex: mouse and human), essentially genes in categories (1,2,3,4), and so to improve the efficiency of drug development by decreasing the compound attrition rate from preclinical trials to clinical trials.

Data analysis

The estimated treatment effects for mice and humans are obtained by the following independent simple linear models:

$$\begin{array}{l}{X}_{a_{\mathit{\text{ij}}}}={\beta }_{0a_i}+{\beta }_{1a_i}{T}_{a_j}+e_{a_{\mathit{\text{ij}}}},\end{array}$$

(5)

$$\begin{array}{ll}{X}_{h_{\mathit{\text{il}}}}^{\mathit{\text{week}}8}& ={\beta }_{0h_i}+X_{h_{\mathit{\text{il}}}}^{\mathit{\text{week}}0}+{\beta }_{1h_i}{T}_{h_l}+{\mathit{\beta }}_{2h_i}\text{PriorTherapy}\\ +{\beta }_{3h_i}\text{ConMed}+e_{h_{\mathit{\text{il}}}},\end{array}$$

(6)

where i = 1,…,11,922 (number of orthologs); j = 1,…,19 (number of mice); l = 1,…,24 (number of humans). $X_{a_{\mathit{\text{ij}}}}$, $X_{h_{\mathit{\text{il}}}}^{\mathit{\text{week}}0}$ and $X_{h_{\mathit{\text{il}}}}^{\mathit{\text{week}}8}$ are gene expressions from the ith orthologous gene for the jth mouse, gene expression from the ith orthologous gene for the lth human before treatment intervention, and gene expression from the ith orthologous gene for the lth human at the completion of the clinical experiment, respectively. $T_{a_j}$ and $T_{h_l}$ are {0,1} treatment indicators, and $e_{a_{\mathit{\text{ij}}}}$ and $e_{h_{\mathit{\text{il}}}}$ are independent $N(0,{\sigma }_a^2)$ and $N(0,{\sigma }_h^2)$ random variables. Additionally, there are two more covariates in the human model: PriorTherapy (a four-level categorical variable indicating patients’ therapy prior to the clinical trial) and ConMed (a two-level categorical variable for concomitant medication status).

Equation (4) was used to model the distribution of the estimated treatment effects. Genes in categories (1, 2, 3, 4) are believed to be potential biomarkers that can greatly improve the design of the process of drug development as these orthologous genes interact with drugs in a way that shows some relationship between two different species (in this case, human and mouse) and so studying the behavior of these genes in preclinical trials might help scientists better understand the mechanism of drug activities in clinical trials.

Results

Parameter estimation

The maximum likelihood estimates of the parameters in the bivariate normal mixture model using the EM algorithm are given in Table 5. The estimated mixture weight for category 0 was ${\hat{\pi }}_0=0.889$ indicating approximately 10,599 (11,922 × 0.889) pairs of uninteresting mouse and human orthologs. (μ _a ,μ _h )^T denotes the mean vector of each mixture component. The estimated treatment effect means of mice were in general larger than that of humans, indicating that the magnitude of the difference of expression in genes due to treatment intervention in mice tended to be larger than in genes where differential expression is exhibited in humans. The estimated variance for mice ${\hat{\sigma }}_{\mathit{\text{ak}}}^2$ tended to be larger than that for humans ${\hat{\sigma }}_{\mathit{\text{hk}}}^2$, suggesting that overall, the variability in the mouse data set was larger than the human data set. The estimated correlation coefficients between the treatment effects for both species were (-0.531,0.380,-1,0.102) for categories (1,2,3,4), respectively. The estimated correlation for category 3 was based on only two observations. Standard errors of the parameter estimates from the two-species experiment were obtained by bootstrapping with 1,000 bootstrap replicates. Based on the bootstrap standard errors, there was little evidence of bias for the parameter estimates.

Table 5 Parameter estimates of the bivariate mixture model

Gene identification

The vector of the number of genes identified for each category was ${({\hat{n}}_0,{\hat{n}}_1,\dots ,{\hat{n}}_8)}^T$ = (10,814, 41, 12, 2, 20, 168, 578, 12, 275)^T. Figure 1 displays the scatter plots of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ before and after gene membership identification, including the scatter plot of ${({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ for all orthologs, orthologs after eliminating the uninteresting ones (orthologs not in category 0), and orthologs reacting to the treatment stimulus for both species (orthologs in categories (1, 2, 3, 4)).

Figure 1

Scatter plots of the estimated treatment effects${\mathbf{(}{\stackrel{\mathbf{ˆ}}{\mathit{\beta }}}_{\mathbf{1}{\mathit{a}}_{\mathit{i}}}\mathbf{,}{\stackrel{\mathbf{ˆ}}{\mathit{\beta }}}_{\mathbf{1}{\mathit{h}}_{\mathit{i}}}\mathbf{)}}^{\mathit{T}}$ before and after gene membership identification. From left to right: all orthologs, orthologs differentially expressed in either species (categories (1,2,3,4,5,6,7,8)), and orthologs differentially expressed in both species (categories (1,2,3,4)).

Based on the bivariate mixture model, among the 11,922 pairs of orthologs, 10,814 pairs did not react to the drug treatment for either species, 75 pairs of orthologs (sum of the gene counts in categories 1 through 4) showed evidence of differentiation between treatments for both species. Genes in categories (5,6,7,8) are also potential candidates for further investigation to improve the process of drug development since studying these genes might uncover the myth of the overall high attrition rate of a drug compound from preclinical trials to clinical trials.

In comparison, an attempt was made to identify differentially expressed human genes using solely the human data, i.e., traditional t statistics were used to test whether or not ${\beta }_{1h_i}=0$ and the approach of Benjamini and Hochberg [19] was used to adjust for multiple comparisons. With the nominal FDR controlled at 0.01, 0.1, 0.5, and 0.9, this single-species method failed to identify any differentially expressed human genes at any levels of FDR. With the p values histograms in Figure 2 showing an obvious difference between the observed significance of differential expression for mice and humans, specifically, the p values were nearly uniformly distributed for humans, the result was not surprising.

Figure 2

Histograms of p values from tests of no treatment effects. From left to right: p value histogram for mice and p value histogram for humans.

Concluding remarks

This research was motivated by a fundamental yet still not well-understood problem in the drug development process. The results obtained in preclinical animal trials do not seem to translate well enough to make inferences for human clinical trials, resulting in an undesirably high attrition rate in human experiments. A bivariate mixture model which utilizes information across two species was proposed to identify genes that exhibit similar patterns of expression across species, with the hope that studying genes could help understand biological differences across species at the molecular level and ultimately help reduce attrition in drug development. It is also of great interest to identify genes active in animal, but not in human since studying this group of genes might lead to answers that explain why some drug trials fail in translation to humans. The comprehensive simulation study suggested that the proposed model, though highly structured, can accommodate various configurations of differential gene expression, especially when the magnitude of differential expression due to different treatment intervention was weak.

In the application of the bivariate mixture model on the GSK type II diabetes experiment, the mixture model was able to separate differentially expressed genes from non-differentially expressed genes. A potential multi-gene predictor may be developed according to the genes identified by the bivariate mixture model to benefit patients in therapeutic decision making.

The mixture model is highly structured, with strong but somewhat simplifying assumptions. The grouping of all genes for which the expression difference is positive in both species into a single category parameterized by a single bivariate mean may be somewhat of a simplistic approach. In practice, the data may not be normally distributed or be comprised by exactly nine groups and lead to bias and inefficiency. Forcing the normality assumption and the grouping may be inefficient. However, by modeling the least squares estimates of the expression differences as bivariate realizations from a distribution with a single mean vector, some flexibility in the model is retained. It is perhaps less important to precisely quantify the magnitude of the expression difference than to determine whether it is positive or negative and whether or not the direction is preserved across species.

Using the proposed mixture model for gene identification is completely data-driven. The simulation study in Section ‘Simulation’ indicated that the mixture model at times report a poor observed false discovery rate due to large variability in measurement of expression. Currently, gene membership classification is determined by maximizing the posterior probability that observation y _i belongs to the kth cluster. It is possible to choose costs to attach weight to different types of misclassification and then choose a classification rule to minimize expected cost. This rule may lead to a more desirable list that can better accommodate goals for future research.

The mixture model currently only handles data with two treatments/cancer types/drugs, i.e., data containing only two-level variables. In practice, it is often the case that an experiment is run under multiple conditions. An extension of the model to handle experiments with factors with three-or-more levels is an important future work.

This approach focuses on identifying orthologs that show same/opposite mechanism between two species. However, genes identified by the bivariate mixture model (genes in categories (1, 2, 3, 4)) may not lead to the most powerful model for prediction of cancer types/response status for either one of the species. The most powerful model may be based on those human genes that lie in categories 5 and 6. For those genes, the corresponding animal genes show no signs of differential expression. Hence, incorporating prediction ability into the model or developing a prediction model that can help utilize the genes selected by the mixture model is of great interest and is an obvious candidate for future research.

Appendix

Mixture models and the EM algorithm

Estimating the parameters in (4) is nontrivial. Redner and Walker [22] offer an excellent review of estimating the parameters which determine a mixture density. In particular, the paper is devoted to a particular iterative procedure for numerically approximating maximum likelihood estimates (MLE) of the parameters in mixture densities. This method was formalized by Dempster et al. [23] and termed the EM (Expectation-Maximization) algorithm, and is used for numerically approximating the maximum likelihood estimates for (4).

The EM algorithm with no constraints

For a finite mixture model with C components, given data y with independent multivariate observations y₁,…,y _g , each y _i is taken to be a realization of the mixture probability density function,

$$f({\mathbf{\text{y}}}_i|\mathbf{\Psi })=\sum _{k=1}^C{\pi }_k{f}_k({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k),$$

where Ψ = (θ₁,…,θ _C ,π₁,…,π _C )^T, a vector of unknown parameters. f _k and θ _k are the density and parameters of the kth component in the mixture, respectively. π _k is the probability that an observation arises from the kth component. Note that π _k  ≥ 0 and ${\sum }_{k=1}^C{\pi }_k=1$. For classification purpose and to achieve minimum misclassification rates, y _i is assigned to the population (category) for which the posterior probability that y _i belongs to the kth cluster (the kth component of the mixture) is maximized. The posterior probability is given by

$${\tau }_k({\mathbf{\text{y}}}_i;\mathbf{\Psi })=\frac{{\pi }_k{f}_k({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k)}{{\sum }_{h=1}^C{\pi }_h{f}_h({\mathbf{\text{y}}}_i|{\mathit{\theta }}_h)}.$$

In summary, the EM algorithm may be implemented to maximize the likelihood of a multivariate normal mixture model by following these two steps:

E-step: The E-step on the (j+1) iteration takes the conditional expectation of the complete-data log likelihood, given the observed data (Q( Ψ|Ψ^(j))).

At the (j + 1) iteration, the E-step results in

$$Q(\mathbf{\Psi }|{\mathbf{\Psi }}^{(j)})=\sum _{i=1}^g\sum _{k=1}^C{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})(log{\pi }_k+log{f}_k({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k)),$$

(7)

where $f_k({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k)={(2\pi )}^{-\frac{p}{2}}{|{\mathbf{\Sigma }}_k|}^{-\frac{1}{2}}{\mathrm{e}}^{-\frac{1}{2}{({\mathbf{\text{y}}}_i-{\mathit{\mu }}_k)}^T{\mathbf{\Sigma }}_k^{-1}({\mathbf{\text{y}}}_i-{\mathit{\mu }}_k)}$ and ${\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})=\frac{{\pi }_k^{(j)}{f}_k({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k^{(j)})}{{\sum }_{h=1}^C{\pi }_h^{(j)}{f}_h({\mathbf{\text{y}}}_i|{\mathit{\theta }}_h^{(j)})}$.

M-step: The M-step on the (j+1) iteration requires the global maximization of Q(Ψ|Ψ^(j)) with respect to Ψ over the parameter space to give the updated estimate Ψ^(j+1).

$${\hat{\pi }}_k^{(j+1)}=\frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}{g};$$

(8)

$${\stackrel{\mathbf{ˆ}}{\mathit{\mu }}}_k^{(j+1)}=\frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)}){\mathbf{\text{y}}}_i}{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})};$$

(9)

$${\stackrel{\mathbf{ˆ}}{\mathbf{\Sigma }}}_k^{(j+1)}=\frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)}){({\mathbf{\text{y}}}_i-{\mathit{\mu }}_k)}^T({\mathbf{\text{y}}}_i-{\mathit{\mu }}_k)}{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}.$$

(10)

The iterations of the EM algorithm continue until some stopping criterion is met, such as the difference of the conditional expectation of the complete-data log likelihood at the (j + 1) step and the conditional expectation of the complete-data log likelihood at the j step is sufficiently small, i.e.,

$$Q({\mathbf{\Phi }}^{(j+1)}|{\mathbf{\Phi }}^{(j)})-Q({\mathbf{\Phi }}^{(j)}|{\mathbf{\Phi }}^{(j)})\le \epsilon .$$

Throughout this research, ε = 0.0001.

The EM Algorithm with constraints

With constraints on the parameter space, the EM algorithm derived in Section ‘The EM algorithm with no constraints’ cannot be used directly. To accommodate the restricted parameter space introduced for (4), first note that (4) is a bivariate normal mixture with nine components and the density for each component can be written as

$$\begin{array}{l}{f}_k\left({\mathbf{\text{y}}}_i|{\mathit{\theta }}_k\right)\\ =\frac{e^{-\frac{1}{2}\frac{1}{{\sigma }_{a_k}^2{\sigma }_{h_k}^2\left(1-{\rho }_k^2\right)}\left({\sigma }_{h_k}^2{\left(y_{a_i}-{\mu }_{a_k}\right)}^2-2{\rho }_k{\sigma }_{a_k}{\sigma }_{h_k}\left(y_{a_i}-{\mu }_{a_k}\right)\left(y_{h_i}-{\mu }_{h_k}\right)+{\sigma }_{a_k}^2{\left(y_{h_i}-{\mu }_{h_k}\right)}^2\right)}}{2\pi \sqrt{{\sigma }_{a_k}^2{\sigma }_{h_k}^2\left(1-{\rho }_k^2\right)}},\end{array}$$

where ${\mathbf{\text{y}}}_i={(y_{a_i},y_{h_i})}^T={({\hat{\beta }}_{1a_i},{\hat{\beta }}_{1h_i})}^T$ and ${\mathit{\theta }}_k={({\mu }_{a_k},{\mu }_{h_k},{\sigma }_{a_k}^2,{\sigma }_{h_k}^2,{\rho }_k)}^T$.

The solution for the E-step remains the same. The membership probability τ _k (y _i ;Ψ^(j)) can be obtained by taking the conditional expectation of the complete-data log likelihood, given the observed data.

As for the M-step, the MLE for π _k remains the same, ${\hat{\pi }}_k^{(j+1)}=\frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}{g}$. However, the MLEs for μ _k and Σ _k need to be modified according to the constrained parameter space for each mixture component.

Essentially, the MLEs on the (j + 1) iteration for (μ _ak ,μ _hk )^T, k = 1,…,4, are

$$\begin{array}{l}{\hat{\mu }}_{a1}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_1({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_1({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_1({\mathbf{y}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{h1}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_1({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_1({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_1({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{a2}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{h2}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_2({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{a3}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{h3}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_3({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{a4}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{h4}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_4({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\end{array}$$

The MLEs of Σ _k , k = 1,…,4, at the (j + 1) M-step remain unchanged as in (10).

The MLEs on the (j + 1) iteration for (μ _ak ,μ _hk )^T, k = 5,…,8, are

$$\begin{array}{l}{\hat{\mu }}_{h5}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_5({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_5({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_5({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{h6}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_6({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}}{{\sum }_{i=1}^g{\tau }_6({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_6({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{a7}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_7({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_7({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_7({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\ge 0;\\ 0& \text{otherwise.}\end{array}\right.\\ {\hat{\mu }}_{a8}^{(j+1)}=\left\{\begin{array}{ll}\frac{{\sum }_{i=1}^g{\tau }_8({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}}{{\sum }_{i=1}^g{\tau }_8({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& \text{if}{\sum }_{i=1}^g{\tau }_8({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}\le 0;\\ 0& \text{otherwise.}\end{array}\right.\end{array}$$

The MLE for Σ _k , k = 0,5,6,7,8, on the (j + 1) iteration is

$${\hat{\mathbf{\Sigma }}}_k=\left(\begin{array}{cc}\frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{a_i}^2}{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}& 0\\ 0& \frac{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})y_{h_i}^2}{{\sum }_{i=1}^g{\tau }_k({\mathbf{\text{y}}}_i;{\mathbf{\Psi }}^{(j)})}\end{array}\right).$$

Regularized covariance matrices in the EM algorithm

The component covariance matrices of a mixture model may be singular or near singular in the EM iterative process. When the covariance matrices corresponding to one or more components are ill-conditioned (singular or near singular), the EM algorithm breaks down. Particularly, applying the EM algorithm for a mixture model with large numbers of components when there are actually fewer groups often results in the failure of the EM algorithm due to ill-conditioning [24, 25]. Indeed, this break-down of the EM algorithm may imply that clusters contain insufficient observations and too many components are used to fit the data set where there are actually fewer clusters, or clusters contain points that are of very little variation compared to other clusters. Hence, an intuitive solution to this is to decrease the number of the mixture components. However, this immediately leads to another question: how many clusters are needed. Though an active area of research, it is beyond the scope of this study. Nonetheless, various approaches have been proposed to generate numerically non-singular covariance matrices [26–33]. Among such, the regularization method proposed by Sato and Ishii [32] has been adopted to obtain numerically non-singular covariance matrices throughout this research. In Sato and Ishii [32], the regularized covariance matrix for the kth mixture component in the (j + 1) M-step is

$${\mathbf{\Sigma }}_{\mathit{\text{Rk}}}^{(j+1)}={\mathbf{\Sigma }}_k^{(j+1)}+\alpha \frac{\text{tr}({\mathbf{\Sigma }}_k^{(j+1)})}{p}{I}_p,$$

(11)

where 0 ≤ α ≤ 1 is a small constant and I _p is a p-dimensional identity matrix. If ${\mathbf{\Sigma }}_k^{(j+1)}$ equals 0, then $\text{tr}({\mathbf{\Sigma }}_k^{(j+1)})$ is set to be a small threshold value (0.0001 in this research).