Bias-Reduced Estimation of Finite Mixtures: An Application to Latent Group Structures in Panel Data

Bias-Reduced Estimation of Finite Mixtures : An Application to Latent Group Structures in Panel Data

\thispagestyle{empty}

Introduction

\setstretch{1.3} \pagenumbering{arabic} Finite mixtures are extensively used in statistics, computer science, and machine learning for pattern recognition and unsupervised classification to account for various types of unobserved heterogeneity bishop_pattern_2006, fruhwirth-schnatter_finite_2006. Several applications of finite mixtures can also be found in labor and health economics heckman_method_1984,deb_demand_1997, keane_career_1997, jones_healthcare_2015. For instance, deb_demand_1997 use finite mixtures to distinguish two unobserved, latent types (i.e. the “healthy” and the “ill”) regarding the demand for medical care and find substantial differences in fitted distributions for each type. Methods to account for unobserved heterogeneity can reduce bias hsiao_analysis_2014, improve inference and forecast boot_optimal_2018, and also allow for the estimation of heterogeneous treatment effects ahn_dierence_2024. A finite mixture distribution is a convex combination of a small number of distinct parametric densities where the combination weights, known as mixing weights, correspond to the proportion of observations that originate from each density in the population. The resulting density, known as the mixture density, fully describes the distribution of the observed data. Finite mixtures are related to unsupervised clustering methods since each component density within the mixture fully describes its corresponding group of observations. Finite mixtures are also known to be highly flexible since they can “approximate any distribution to a sufficient level of accuracy” and easily accommodate the presence of covariates and nonlinear components nguyen_approximations_2019. Estimation of the parameters contained in the mixture density is usually performed using maximum likelihood estimation (MLE) procedures, including nonparametric MLE compiani_using_2016, kitamura_nonparametric_2018. Although it has been shown that maximizing the (constrained) likelihood of a finite mixture of parametric densities yields consistent estimates, little is known about the finite-sample behavior of the MLE of finite mixtures balakrishnan_statistical_2017,kwon_minimax_2021. In this paper, I show that maximizing the likelihood of a mixture density using standard parametric MLE methods can lead to finite-sample biases that are much larger than those that are typically expected to be obtained from standard MLE approaches under generic identifiability and correct specification of the mixture density. The biases come from the presence of outliers that are generated by the tails of component densities with unbounded or large support. When components are not sufficiently separated, the global maximum of the likelihood shifts away from the true parameters such that at least one component density becomes strongly influenced by these “abnormal observations” chen_statistical_2023. The presence of large finite-sample biases in the mixture parameters is confirmed via Monte Carlo simulations under mixtures of normal, Poisson, and exponential distributions. The size of the bias depends on the components' size and on the degree of separation between each component density.

The contributions of the paper are twofold. First, I show why standard MLE of finite mixtures can lead to severely biased estimates of all parameters in the mixture under weak regularity conditions and generic identifiability fruhwirth-schnatter_finite_2006.\footnote{By “all parameters” in the mixture, I refer to the mixing weights and the set of parameters governing each density in the mixture. Note also that “components” and “groups” are used interchangeably throughout the paper.} Given the analytical complexities associated with the general behavior of mixture densities in finite samples, I only present the general intuition behind the existence of the finite-sample bias using the two-component mixture of normal distributions as a practical example. Extension of the same intuition to mixtures of other distributions in the exponential family is straightforward. By standard MLE of finite mixtures, I refer to procedures that attempt to globally maximize the mixture likelihood function, as defined in Section (ref). The well-known expectation-maximization (EM) algorithm of dempster_maximum_1977 and all Newton-type algorithms, when applied to the mixture likelihood, both fall into this category. The second contribution of the paper is to show under which conditions the maximization of a different objective function, the classification-mixture likelihood (C-ML) mclachlan_iterative_1975,mclachlan_finite_2019, will lead to less biased, consistent, and asymptotically efficient estimates of all parameters in any generically identifiable mixture models. Similar to dzemski_convergence_2021, I show that consistency in estimation and asymptotic efficiency can be obtained by maximizing the C-ML function if it is combined with a consistent classifier, meaning that the proportion of misclassified observations goes to zero in the limit when evaluated at the true parameter values. Although MLE is asymptotically unbiased in general, it is not unbiased in finite samples. Therefore, obtaining unbiased estimates (up to a term of order $N^{-1}$) of all parameters in the mixture in finite samples is impossible without bias-reduction techniques kosmidis_bias_2014. In the case of finite mixtures, I show that it is possible to reduce the size of the bias from the standard MLE by classifying all observations in their true component density. This latter requirement is much stronger than the one used for consistency, and it cannot be guaranteed for all data-generating processes (DGPs) with a given sample size. Nonetheless, conditions under which maximizing the C-ML objective will lead to the correct classification of all observations in the sample are derived. Such conditions can take the form of restrictions on group memberships, as typically done in the literature on latent group panel structures bonhomme_grouped_2015, su_identifying_2016, or through additional assumptions on the general structure of the outcome vector and/or the available covariates, which is the general strategy that is retained in this paper due to its greater flexibility. I illustrate those two theoretical contributions with both simulated and real-world data. First, I use simulated data to characterize the finite-sample biases of standard MLE procedures when applied to mixtures of normal, Poisson, and exponential distributions. Then, I use both simulated and real-world data to compare the finite-sample performance of the proposed estimation strategy to standard MLE procedures via the EM algorithm. Both types of data correspond to latent group panel structures with a univariate outcome and $p$ strictly exogenous covariates. Simulation results confirm that estimates obtained from the maximization of the mixture likelihood are more biased and less efficient than those obtained from the proposed estimation strategy under weak regularity conditions. On the contrary, reducing the number of covariates or increasing the overlap between component densities leads to a higher misclassification rate, which results in larger finite-sample biases. Simulation results show that no algorithm outperforms the other when the misclassification rate increases away from zero. The empirical part of the paper uses the EM algorithm and the proposed estimation algorithm to model individual healthcare expenditure (HCE) from administrative data over time. Both algorithms use a latent group two-part model (LGTPM) for estimation and an identical specification for each component density. Results show that the proposed algorithm leads to a reduction in out-of-sample prediction error of 17.6% compared to the best result obtained when using the EM algorithm, and of 56.6% compared to the single-component two-part model. The proposed algorithm also allows the recovery of the true group memberships for almost every observation in the sample, which can be modeled in a second step for forecasting purposes. The remainder of the paper is as follows. Section (ref) briefly reviews the related literature. Section (ref) shows why maximizing the mixture likelihood leads to large finite-sample biases in simple mixture models under weak regularity conditions. It also details the proposed estimation strategy that is used to reduce the size of such biases. Section (ref) presents simulation scenarios and results, confirming the theoretical insights from the previous section. Section (ref) presents the empirical application and the related results, while Section (ref) concludes.

Related Literature

Finite Mixtures and the EM Algorithm. Although there exist several optimization algorithms that can maximize the mixture likelihood, the EM algorithm is almost always used by practitioners when estimating mixture models due to its simplicity and empirical performance kwon_global_2024.\footnote{Note that the EM algorithm is not a maximization algorithm per se, but a general strategy for the maximization of any incomplete-data likelihood function. Indeed, it is often desirable to use a numerical optimization method at each M-step of the EM algorithm. See Section (ref) for an example, and Section 2.4.4 of fruhwirth-schnatter_finite_2006 for more details on the subject.} The EM algorithm consists of the consecutive repetition of an expectation step (i.e., the E-step) and a maximization step (i.e., the M-step), each conditional on the results obtained from the previous step. The E-step computes assignment probabilities for each observation and each component density, whereas the M-step maximizes the likelihood function conditional on the most recent assignment probabilities. In their seminal paper, dempster_maximum_1977 showed that such an algorithm never decreases the incomplete-data likelihood function between two consecutive iterations, the mixture likelihood being a special case of the former when the mixing distribution has discrete support. It is also well known that the mixture likelihood is not a convex function of the parameter space mclachlan_finite_2000. Hence, the multimodal nature of the mixture likelihood makes it difficult for any optimization algorithm to find its global maximum, which has given rise to various procedures with the explicit goal of reducing the probability of being “trapped” in a local spurious maximum fruhwirth-schnatter_em_2019. Moreover, it is acknowledged that “the sample size [...] has to be very large, before asymptotic theory of maximum likelihood applies”, especially when the “component densities are poorly separated” redner_mixture_1984,aragam_uniform_2023. To provide statistical guarantees, balakrishnan_statistical_2017 derived sufficient conditions for the initialization of the EM algorithm such that it converges “geometrically to a fixed point that is within statistical precision of the unknown true parameter” balakrishnan_statistical_2017. However, these conditions do not imply that the true parameters will be close to the global maximum of the mixture likelihood for any fixed sample size and any true mixture density. In fact, I show in the next section that searching for the global maximum of the mixture likelihood can yield highly misleading results in finite samples due to the presence of “natural” outliers, which correspond to outliers that are generated by the tails of at least one component density. This intuition is confirmed via simulations that are presented in Section (ref). This paper also relates to a series of papers that followed the publication of balakrishnan_statistical_2017, and where the authors studied the convergence properties of the EM algorithm when applied to mixtures of Gaussian densities and linear regressions wu_randomly_2022, kwon_em_2019, kwon_em_2020, kwon_global_2024. Globally, these papers show that the EM algorithm will asymptotically converge to the true parameter values with high probability -- regardless of the initialization procedure -- for two-component Gaussian and linear regression mixtures with equal mixing weights and known variances. However, the general behavior of the EM algorithm in finite samples remains elusive when mixing weights are unequal or variances are unknown, which is often the case in practice. \\ K-means and the C-EM Algorithm. The K-means algorithm is one of the most widely used clustering algorithms in unsupervised machine learning hastie_elements_2009. It can be used to estimate Gaussian mixture models just as the EM algorithm, but is less flexible than the EM since it assumes spherical covariance matrices of the outcome(s). Moreover, it is known to yield inconsistent estimates of all parameters in the mixture unless the cluster centers are infinitely distant from each other pollard_strong_1981,bryant_large-sample_1991. The C-EM algorithm (for Classification-EM) is the likelihood generalization of K-means and is as flexible as the EM, but still leads to inconsistent estimates of the mixture parameters bryant_asymptotic_1978, celeux_classification_1992.\footnote{The acronym “C-EM” is used to avoid potential confusion with the Coarsened Exact Matching (CEM) procedure proposed by iacus_causal_2012 or with the Component-wise EM algorithm. It also mimics the C-LASSO acronym of su_identifying_2016, even though the “CEM” acronym was originally employed by celeux_classification_1992. Thanks to Giuseppe Moscelli for bringing this point to my attention.} Since the K-means algorithm corresponds to the special case of the C-EM algorithm when all within-cluster errors are independent, normally distributed, and homoskedastic, the rest of the paper focuses on the general and more flexible C-EM algorithm. Contrary to the EM algorithm, the C-EM algorithm classifies each observation into the group that maximizes its corresponding density value. This practical difference leads to an important theoretical difference, which is that the C-EM algorithm maximizes the C-ML function rather than the mixture likelihood function. The C-ML function is the estimated counterpart of the \textit{complete-data} likelihood function, which is the likelihood function one would get if all group memberships were known mclachlan_finite_2019. If this is the case, then it is said that the estimator achieves \textit{oracle efficiency} and all estimated parameters are asymptotically normal and efficient as the corresponding maximum likelihood estimates reach the semiparametric efficiency bound newey_semiparametric_1990,su_identifying_2016. Akin to the EM algorithm, it has been shown that the C-ML function never decreases between two consecutive iterations of the C-EM algorithm.\footnote{This is however only true for a certain class of classifiers. See, for instance, celeux_classification_1992, bottou_convergence_1994, and the proof of Lemma (ref) for more details.} Given that the C-ML function is also multimodal, convergence of the C-EM algorithm to the global maximum of the C-ML function is rarely guaranteed, and similar issues are encountered in practice when maximizing either the mixture likelihood or the C-ML function same_online_2007. It is also important to stress that global convergence of the estimation algorithm is not very helpful if the global maximum of the objective function does not lie in a shrinking neighborhood of the true parameter values when the sample size remains below a certain threshold. Simulation results presented in section (ref) shows that such a phenomenon occurs with finite mixtures of parametric densities that overlap significantly. \\ \textbf{Latent Group Panel Structures}. Several recent studies in econometrics have used different variants of the K-means algorithm to account for unobserved heterogeneity in panel data through the introduction of grouped fixed-effects (GFE) and/or group-specific coefficients bonhomme_grouped_2015,bonhomme_distributional_2019,bonhomme_discretizing_2022,okui_heterogeneous_2021,lumsdaine_estimation_2023, liu_identification_2020, wang_panel_2024. The GFE estimator is also closely related to factor models since it allows time-fixed effects to vary across groups, similar to models with interactive fixed effects bai_panel_2009. However, factor models typically assume homogeneity of all the other parameters in the model, whereas the GFE estimator imposes restrictions on group memberships or on the form taken by the unobserved heterogeneity (across units or over time) to achieve consistency bonhomme_grouped_2015,bonhomme_discretizing_2022, wang_panel_2024. Contrary to the recent developments on the subject, the estimation procedure proposed in this paper leaves group membership completely unrestricted across units and over time, while still achieving consistency under the “many covariates” asymptotics where the number of covariates $p$ is allowed to grow at a strictly slower rate than the sample size cattaneo_inference_2018, kitamura_nonparametric_2018. Several authors also relied on the LASSO device or on binary segmentation algorithms for combining unit-level coefficients into group-level coefficients su_identifying_2016, su_sieve_2019, qian_shrinkage_2016, wang_homogeneity_2018, wang_identifying_2021. All estimation strategies presented in this strand of the literature impose certain restrictions on group memberships over time. For instance, in the case of the Classifier-LASSO originally developed by su_identifying_2016, leaving group membership completely unrestricted implies the estimation of $NT$ parameters for each time-varying covariate in the sample, which is impractical in regular panel datasets. Contrary to this strand of literature, this paper focuses exclusively on regular finite mixtures where no restriction is imposed on group memberships along any dimension of the dataset.

Maximum Likelihood Estimation of Finite Mixtures

This section is divided into three parts. Section (ref) presents the notation that will be used throughout the paper. Section (ref) then intuitively explains why standard MLE of finite mixtures leads to sometimes heavily biased estimates of all parameters under weak regularity conditions in finite samples. For simplicity, the case that is considered in Sections (ref) and (ref) is a cross-section of $N$ independent observations without covariates. Extension of the following results to panel data or time series with strictly exogenous covariates is straightforward. Finally, Section (ref) describes the proposed estimation strategy and shows under which conditions it will lead to less biased, consistent, and asymptotically efficient estimation of all parameters in the mixture. All proofs are deferred to Appendix (ref).

General Setup and Notation

The mixture density function of any observation $y_{i} \in \mathcal{Y} \subseteq \mathbb{R}^{d_y}$ is generally defined as follows

where $y_{i}$ is the (possibly multivariate) observed outcome of individual $i$ with $i = 1,...,N$. The function $f: \mathcal{Y} \subseteq \mathbb{R}^{d_y} \rightarrow \mathbb{R}_{>0}$ is the mixture density, defined as a function of $\theta = (\theta_1,...,\theta_g,...,\theta_G)$, which represents the set of component parameters, and of $\pi = (\pi_1,...,\pi_G)^{\top}$, the vector of mixing weights, where $\xi = (\theta,\pi)$ represents the whole set of parameters, and where $G \in \mathbb{N}_{>1}$ is the total number of components. The set of component parameters $\theta$ is assumed to lie inside the compact parameter space $\Theta = \Theta_1 \times ... \times \Theta_G \subset \mathbb{R}^{G \times d_{\theta_g}}$, whereas the vector of mixing weights $\pi$ is assumed to lie within the open space $\Pi = (0,1)^{G}$, where the only additional constraint imposed on the mixing weights is that $\sum_{g=1}^G \pi_g =1$. Each component density $f_g(.)$ corresponds to a well-defined parametric density with respect to a $\sigma$-finite measure, generally denoted by $\upsilon(d y_{i})$. All component densities do not necessarily belong to the same family of distributions, although this is common in practice to facilitate the estimation of $\theta$, and can be extended to include (exogenous) covariates provided an appropriate specification of the link between $y_i$ and $x_i$ for each component density. Throughout the paper, the set of true parameter values is denoted by $ \xi^0 = (\pi^0, \theta^0)$ and lies in the interior of the true parameter space $\Pi \times \Theta\backslash E \equiv \Xi^0$, where $\pi^0 = (\pi^0_1,...,\pi^0_g,...,\pi^0_G)^{\top}$, $\theta^0 = (\theta^0_1,...,\theta^0_g,...,\theta^0_G)$, and where I exclude the set $E = \{\Theta: \cup_{j\ne g}\cup_{g=1}^G (\theta^0_j = \theta^0_g) \}$ from the true parameter space to avoid identification issues. Any dataset generated by this true density is denoted by $\mathbf{y}_N = \{y_{i}\}_{i=1}^N$. Note that the true number of components, $G$, is assumed to be discrete, finite, and known by the econometrician unless stated otherwise. This assumption allows us to rule out identification issues caused by the misspecification of $G$, a topic that has been extensively covered in both statistics and econometrics redner_note_1981,kasahara_testing_2015,budanova_penalized_2025. In the finite mixture framework, each observation is assumed to originate from only one of the $G$ densities. The unobserved, true binary assignment (or grouping) variable, denoted by $z^0_{ig}$, is defined as follows

Consequently, the true mixing weight of the $g^{th}$ component density, $\pi_g^0$, is defined accordingly

which corresponds to the unconditional probability of any observation to belong to the $g^{th}$ component/group. Note that $z^0_{ig} = 1$ for at least one value of $i \in \{1,2...,N\} = [N]$ for any $g \in \mathbb{G}$ since $\pi^0_g \in (0,1)$ for any $g \in \mathbb{G}$. This also implies that the value $z^0_{ig}$ can be seen as the realization of the random variable $Z_{i}$ drawn from a univariate multinomial distribution with $G$ categories and vector of probabilities $(\pi^0_1,...,\pi_G^0)$ mclachlan_finite_2019. Several objective functions can be written using the above equations. The three most important objective functions are described below in Definition (ref).

The complete-data log likelihood function is the objective one would obtain were the true group memberships known for every observation in the sample. Note that maximizing $l^C(\xi,\mathbf{z}^0)$ with respect to $\theta$ is similar to maximizing $G$ distinct log likelihood functions. The vector of mixing weights may then be estimated using Eq. ((ref)). Since it is assumed that $z^0_{ig}$ is unobserved for all $i \in [N]$, the objective functions shown in Definition (ref) $(ii)$ and $(iii)$ can be used to estimate both $\pi^0$ and $\theta^0$. If $z_{ig} = z^0_{ig}$ for all pairs $(i,g) \in [N] \times \mathbb{G}$, then $l^{CM}(\theta, \mathbf{z}) = l^C(\xi,\mathbf{z}^0) - \sum_{i=1}^N z^0_{ig} \log(\pi_{g})$ and oracle efficiency is achieved for any value of $N$ since maximizing $l^{CM}(\theta, \mathbf{z})$ with respect to $\theta$ will yield estimates as efficient as the ones obtained from maximizing $l^C(\xi,\mathbf{z}^0)$ with respect to $\theta$. Consequently, the estimates $\hat{\theta}^{CM} = \arg \max_{\theta \in \Theta} l^{CM}(\theta, \mathbf{z})$ will be consistent and reach the semiparametric efficiency lower bound asymptotically if $f_g(y_{i}|\theta)$ is regular for any $g \in \mathbb{G}$, meaning that $f_g(y_{i}|\theta)$ has to be a smooth function of $y_i$ for any $\theta \in \Theta$ with nonsingular information matrix (see newey_semiparametric_1990 for more details). Note also that all three objective functions do not put any restriction on group memberships across units, and that the estimated membership $\hat{z}_{i}(\theta^0):= \arg \max_{g \in \mathbb{G}} f_g(y_i|\theta^0)$ will not always lead to $\hat{\mathbf{z}}(\theta^0) = \mathbf{z}^0$ for any value of $N$ unless all component densities $f_g(\cdot|\theta^0)$ are infinitely distant from each other bryant_large-sample_1991,celeux_classification_1992.\footnote{An infinite distance between all component densities implies that the Kullback-Leibler divergence between each pair of different component densities is equal to infinity, meaning that $\int_{\mathcal{Y}} f_g(y_{i}|\theta^0_g) \times f_j(y_{i}|\theta^0_j) \upsilon(dy_{i})$ is arbitrarily close to zero for any pair $(g,j) \in \mathbb{G}\times\mathbb{G}\backslash g$ .} Although it might not be efficient in finite samples, it is commonly acknowledged that maximizing the mixture likelihood function $l^{ML}(\xi)$ with respect to $\xi$ will yield consistent and asymptotically efficient estimates of the mixture parameters provided appropriate constraints on the parameter space redner_mixture_1984, chen_consistency_2017, chen_statistical_2023. However, it remains unclear in the literature on finite mixtures if the estimated parameters $\hat{\xi}^{ML} = \arg \max_{\xi \in \Pi \times \Theta} l^{ML}(\xi)$ feature large biases in finite samples, and what would drive the size of such biases. The next section attempts to shed some light on these aspects.

Finite-Sample Bias of MLE of Finite Mixtures

This subsection shows why maximizing the mixture likelihood with respect to $\xi$ sometimes leads to large finite-sample biases in all estimated parameters under the following regularity conditions.

\setstretch{1.3} Assumption (ref)$(i)$ assumes that all observations are independently and identically distributed according to the true mixture density, which is standard in the literature. Assumption (ref)$(i)$ however rules out cases where the log likelihood function is not bounded from above. As emphasized recently by chen_statistical_2023, this assumption is crucial for the consistency of the standard MLE of finite mixtures. In practice, it rules out standard MLE of mixtures with component densities that belong to the location-scale family tanaka_strong_2009. This issue is well-studied in the particular case of a mixture of Gaussian densities, where it is possible to set the mean value of the first component equal to any observation in the sample while making the variance of this same component arbitrarily close to zero, which leads to an infinite density value li_non-finite_2009. To bound the likelihood function from above, hathaway_constrained_1985 proposes to bound from below the minimum ratio of the variances from any two component densities in the mixture. If the true parameters lie in the constrained parameter space, redner_note_1981 showed that maximizing the constrained mixture likelihood will yield consistent estimates of all parameters in the mixture. However, constraining the mixture likelihood does not guarantee small biases or estimation variances of the estimated parameters in finite samples. Furthermore, it is not clear how one should determine the exact value of the lower bound for a given sample size tanaka_strong_2009. Instead of bounding the minimum ratio of any two variances, tanaka_strong_2009 and chen_statistical_2023 suggest using penalized MLE with a sufficiently severe penalty -- with appropriate upper and lower bounds -- to ensure consistency of all parameters in the mixture. Nevertheless, little is known about the finite sample properties of such a strategy given the multimodal nature of the penalized likelihood function. Assumption (ref)$(ii)$ states that the value taken by every component density is restricted to be strictly positive on the whole support, where this support is assumed to be identical across densities. Note also that this assumption implies that the support of $f_g(\cdot)$ is unbounded, which excludes mixtures of distributions with bounded support, such as the uniform, beta, and Dirichlet distributions. Although this assumption could be relaxed by enforcing bounds on the support of all component densities, these bounds has to be large enough such that there exists outliers $\tilde{y}_i$ that lead to $f_g(\tilde{y}_i|\theta^0_g) > 0$ with $\left|\frac{\tilde{y}_i-\mathbb{E}[y_i]}{\sigma_y}\right| >> 0$ and $\sigma^2_y= Var[y_i]$. Assumption (ref)$(ii)$ also rules out cases where any two component densities have identical parameter values, which prevents the identification issues arising from overspecification of the number of components budanova_penalized_2025. Assumption (ref)$(iii)$ is standard in MLE problems and is satisfied for several families of parametric densities. This assumption guarantees that $\hat{\theta}_{g}:= \arg \max_{\theta \in \Theta} \sum_{i=1}^N \log(f_g(y_{i}|\theta))$ is point identified and always corresponds to a single set of estimated parameter values. Assumption (ref)$(iv)$ is commonly known as generic identifiability of the mixture density. Section 1.3 of fruhwirth-schnatter_finite_2006 describes the three types of identification issues that can arise when modeling finite mixtures, including generic identifiability. The two other identification issues are characterized as “weak” since they can be overcome by using appropriate constraints on the component labels and the number of components. Specifically, teicher_identifiability_1961 and yakowitz_identifiability_1968 show that mixtures of Gaussian, Gamma, exponential, Poisson, negative binomial, and binomial distributions (under certain conditions) are generically identifiable provided that the number of components $G$ is correctly specified. Note that the absence of covariates also avoids the intra-component label switching issue that arises when label switching in covariates leads to identification issues grun_finite_2008. Assumption (ref)$(v)$ is similar to Assumption (ref)$(iii)$, but applied to the mixture likelihood function. It is different from Assumption (ref)$(iii)$ since it eliminates situations where the mixture likelihood would have two identical global maxima, but at two different locations in $\Xi$ for a given dataset $\mathbf{y}_N$. Finally, Assumption (ref)$(vi)$ is made for convenience and is naturally satisfied for most mixtures of well-defined parametric densities. Although there exists some well-known cases where the Fisher information matrix $\mathcal{I}(\xi)$ is not finite for given values of $\xi \in \Xi$, such pathological cases occur when $\xi$ is located on the boundary of the parameter space and are therefore unlikely to be encountered in most practical settings li_non-finite_2009. If Assumption (ref) is satisfied, it is widely recognized that the standard maximum likelihood estimate $\hat{\xi}^{ML} := \arg \max_{\xi \in \Xi} l^{ML}(\xi)$ is $\sqrt{N}$-consistent, and asymptotically normal and efficient with chen_statistical_2023 :

However, this distributional behavior is only asymptotic, and it is unclear if the above normal distribution represents a good approximation uniformly across the parameter space $\Xi$ in finite samples. As noted by redner_mixture_1984, maximizing the mixture likelihood with respect to $\Xi$ can often be characterized as an ill-conditioned problem in the sense that “its solution is very sensitive to some perturbations in the data while being relatively insensitive to others” redner_mixture_1984, which often makes the asymptotic distribution shown above a very bad approximation of the finite-sample distribution of MLE of finite mixtures. Such a large discrepancy between the finite-sample and asymptotic behavior of the estimator arises because the global maximum of the mixture likelihood function need not be located near the true parameter values $\xi^0$ for any given sample size $N$ below a certain threshold. This intrinsic yet overlooked feature of the MLE of finite mixtures is exemplified below for a mixture of two normal distributions.

Example (ref) shows that the MLE of finite mixtures can mistake apparent contamination for $\epsilon$-contamination in finite samples huber_robust_1964, at least for mixtures of normal distributions. $\epsilon$-contamination refers to the situation where the observed data originates from the following “contaminated” density

where $\epsilon > 0$ is close to zero and where $f_{\epsilon}(\cdot|\theta_{\epsilon})$ represents the “contamination” distribution, which is governed by the set of parameters $\theta_{\epsilon}$. In practice, mixture models are well-suited to estimate the parameters of an $\epsilon$-contaminated density, but can also lead to higher likelihood values in finite sample by fallaciously treating natural outliers from well-balanced mixtures as a product of $\epsilon$-contamination. Contrary to true cases of $\epsilon$-contamination, natural outliers will asymptotically be “smoothed out” by the tails of the mixture density. Therefore, the probability of observing a large difference between the estimates and the true parameter values will depend on the sample size and on the tails of the true density : mixture densities with “fatter” tails will generate more extreme outliers in finite samples, and more extreme outliers will lead to a larger probability of observing a large difference between the estimated parameters and the true values. Bounding the component variances from below only puts a mild restriction on the number of outliers included in the contamination density.

The effect of these “natural” outliers on the mixture likelihood is illustrated in Figure (ref). For values of $N$ that are below a certain threshold $N^*$, the true parameters $\xi^0$ can be arbitrarily far away from the global maximum for any $N \le N^*$, depending on the location and the number of outliers generated by the true mixture density, as shown in Panel (a) of Figure (ref). As $N$ remains below $N^*$, the $\epsilon$-contamination model may better represent the observed data than a mixture model evaluated near the true parameters. As $N$ increases and surpasses $N^*$, natural outliers are “smoothed out” and the consistency of $\hat{\xi}^{ML}$ lead to a displacement in the location of the global maximum of $l^{ML}(\xi)$, which will then become much closer to the true parameter values, as shown in Panel (b) of Figure (ref). Therefore, the consistency property of $\hat{\xi}^{ML}$, combined with Assumption (ref)$(v)$, leads to an asymptotic global maximum of $l^{ML}(\xi)$ that is located in a shrinking neighborhood of the true parameter values. However, this questions the uniform consistency of $\hat{\xi}^{ML}$ over $\Xi^0$ since $\hat{\xi}^{ML}$ may converge at an arbitrarily slow rate as $\xi^0$ approaches the boundary of the parameter space $\Xi^0$. Validating this last intuition goes beyond the scope of this paper. In summary, the previous example and the above discussion demonstrated that it is impossible to guarantee that the bias of $\hat{\xi}^{ML}$ will shrink toward zero as $N$ increases but remains below $N^*$ for a given $\xi^0 \in \Xi^0$. The bias could even increase with the sample size, as shown in Figure (ref) of Section (ref). Note that the value of the threshold $N^*$ is data-dependent and typically increases with the total degree of overlap between the true component densities redner_mixture_1984. Given the analytical complexities related to the finite-sample distribution of MLE of finite mixtures, formally deriving the analytical expression of such a bias is left to future research.

Bias-Reduced Estimation of Finite Mixtures

Although it is asymptotically unbiased under standard regularity conditions, it is well known that the maximum likelihood estimator exhibits finite-sample bias even for well-defined, single-component densities. Formally, the bias of the MLE can be generally represented as follows firth_bias_1993

for a general set of true parameters denoted by $\beta$, and where $\frac{b_d(\beta)}{N^d}$ with $d \in \mathbb{N}_{>0}$ represents the $d^{th}$-order term of the bias. Several bias-reduction methods exist in order to get rid of the first-order term $\frac{b_1(\beta)}{N}$ : direct bias-correction, jackknifing, bootstrap, etc.; see kosmidis_bias_2014 for an introductory review. The bias of the MLE of finite mixtures is different. As argued in the previous subsection, the finite-sample bias of the MLE of finite mixtures might be represented by Eq. ((ref)) only once $N$ is above a certain threshold $N^*$ due to the multimodal nature of the mixture likelihood. Given that $N^*$ is data-dependent, it is not clear if this finite-sample bias necessarily reduces when $N$ increases but remains below $N^*$. It is also not clear how traditional bias-reduction techniques can be directly applied to finite mixtures estimated by traditional MLE techniques, such as the EM algorithm. The strategy that I am proposing to reduce the bias in the MLE of finite mixtures relies on the correct classification of almost every observation in the sample. By correct classification, I mean the estimation of the true group memberships $\mathbf{z}^0$ such that the estimated memberships $\hat{\mathbf{z}} = \mathbf{z}^0$ up to an appropriate permutation in the group labels for any given sample size $N$. Once this is achieved, the finite-sample bias can be generally represented by Eq. ((ref)) for each group separately. A second bias-reduction step may then be applied to the within-group estimates using traditional bias-reduction methods. The rest of this article focuses exclusively on the first step of this two-step procedure. The proposed bias-reduction strategy relies on the availability of a consistent unsupervised classification procedure. More precisely, I show that the use of the Mahalanobis distance leads to a consistent classifier such that, under the assumptions stated below, its misclassification rate goes down to zero as $p$ increases, where $p$ can either represent the size of the set of available (exogenous) covariates $x_{i}$ if the density of interest is the conditional mixture density $f(y_i|x_i;\xi)$, or the size of the outcome vector when $\dim(y_i) > 1$. In this section, I consider the second case where $y_i$ is a $p$-sized vector with no covariate for notational simplicity. All results can however be easily extended to the first setup with $f(y_i|x_i;\xi)$, where $p = \dim(x_i)$ and $x_i$ is used instead of $y_i$ for classification (see Definition (ref)). In either case, additional information is required to reduce the misclassification rate and the estimation bias in finite samples. If no additional information is available, then traditional MLE can still be applied -- namely via the EM algorithm -- but this can lead to severe biases if the sample size is not sufficiently large or if the true component densities feature a strong degree of overlap. Another option would be to restrict the group memberships across observations, but it is unclear which restrictions should be applied to guarantee consistency of the estimates in cross-sectional data.

\setstretch{1.3} Assumption (ref)$(i)$ is the analog of Assumption (ref)$(i)$ but adapted to the C-ML objective, where the expected log likelihood of the $g^{th}$ component density is bounded from above. It is worth noting that this assumption is less restrictive than Assumption (ref)$(i)$ since there is no value of $\theta_g \in \Theta_g$ and $y_i \in \mathcal{Y}$ such that $f_g(y_{i}|\theta_g) = \infty$ if $f_g(\cdot|\theta_g)$ is a well-defined, single-component density. Assumption (ref)$(ii)$ is analogous to Assumption (ref)$(iv)$ but expressed at the level of the components. This assumption implies that the C-ML objective will be generically identifiable if each $\theta_g$ is point identified. Contrary to the mixture likelihood, overspecifying the number of components in the C-ML objective does not lead to any identification issue since there is no mixing weight involved in the function. In practice, it is possible to test for the presence of a “useless” component by testing if $\theta_g = \theta_j$ for any pair $(g,j) \in \mathbb{G}\times \mathbb{G}\backslash g$. Contrary to standard MLE of finite mixtures, it is possible to do so with a likelihood ratio test without facing any non-standard asymptotics provided that $\theta_g$ lies in the interior of $\Theta_g$ chen_likelihood_2014. Testing the presence of a “useless” component in this context is also facilitated by the notion that estimators based on classification objectives tend not to mix observations from different groups when $G$ is overspecified liu_identification_2020. Assumption (ref)$(iii)$ is the component-level analog of Assumption (ref)$(vi)$, and is slightly weaker than Assumption (ref)$(vi)$ given the absence of any mixing weight in the density $f_g(y_i|\theta_g)$ and the score function $s_i(\theta_g)$. Assumption (ref)$(iv)$ states that the true mean values and variance-covariance matrix of $y_{i}$ are identical within each group, and that only the variance-covariance matrix has to differ across groups. This is not stronger than assuming that $y_i$ is distributed according to $f_g(\cdot|\theta^0)$ if and only if $z^0_{ig} = 1$, which has already been assumed in Section (ref). This assumption however states that each $\Sigma^0_g$ is positive-definite, which implies that no element of $y_i$ is perfectly collinear with any other set of outcomes within each group. This is crucial to guarantee that $\Sigma^0_g$ is invertible for any $g \in \mathbb{G}$. The same assumption also implies the absence of cross-sectional dependence between vectors of outcome values, which could be relaxed at the expense of tedious theoretical complications (see the proof of Corollary (ref)). It also assumes that the fourth moment of each element in $y_{i}$ is bounded from above, which is naturally satisfied for most parametric densities. Finally, Assumption (ref)$(v)$ implies that the proportion of elements that differ between each lower triangular matrix $W_g$, where $W_g$ is such that $\{\Sigma^0_g\}^{-1} = W_gW_g^\top$ for any $g\in \mathbb{G}$, does not go to zero as $p \to \infty$. This assumption puts a relatively weak restriction on each true density $f_g(\cdot|\theta^0_g)$ since it can be satisfied even when a relatively small and fixed number of elements in the covariance matrix $\Sigma^0_g$ are different across groups. However, the magnitude of the differences in $\Sigma^0_g$ across groups will influence classification error in finite samples depending on the chosen classifier (or group membership indicator), which is generally defined below.

Although they are ruled out by assumption (ref) $(ii)$, eventual ties between values of $h_j(y_i|\theta_j)$ require a tie-breaking rule in order to complete the above definition, but it is left unspecified for simplicity. Using Definition (ref), the bias of any classifier $z_{ig}(\theta)$ is defined as follows.

Remark 1. An unbiased classifier does not correspond to a classifier such that $\mathbb{E}_0[z_{ig}(\theta^0)] = \mathbb{E}_0[\arg\underset{j\in \mathbb{G}}\max \ h_j(y_{i}|\theta^0_j)] = z^0_{ig}$, which actually corresponds to a uniformly consistent classifier (see Definition (ref)). Instead, the bias of any classifier refers to the capacity of its discriminant function to correctly classify any observation after averaging $h_j(\cdot|\theta^0_j)$ over an infinite number of random draws from $f(y_i|\theta^0_j)$. A consistent classifier is then defined as follows.

Akin to su_identifying_2016, I now define a classifier that is uniformly consistent.

Remark 2. A uniformly consistent classifier is a classifier that groups every observation into its true group with probability approaching 1 as $N \to \infty$ : it is consistent uniformly over all possible realizations of $y_{i}$. Note that a uniformly consistent classifier will not necessarily yield a misclassification rate that is equal to zero in finite samples. Analogously, a misclassification rate equal to zero in finite samples does not necessarily imply that the employed classifier is uniformly consistent. If the rate of convergence of a uniformly consistent classifier is (very) fast, then it should lead to an estimated misclassification rate $\hat{E}_{Np}(\theta)$ that is (very) close to zero when $\theta$ is close to or equal to its true value $\theta^0$. Unlike su_identifying_2016, the latter definition does not distinguish between Type I and Type II errors since they can be controlled simultaneously in finite samples if no observation is left unclassified, which is implicitly assumed throughout the paper due to Definition (ref). I now define three different classifiers : the joint density classifier, the Euclidean distance classifier, and the Mahalanobis distance classifier.

Using the joint density classifier is a natural choice when an appropriate distributional assumption is made on the outcome vector. This classifier is known as the Bayes' classification rule and reduces to a naive Bayes' classifier when all covariance elements in $\Sigma = (\Sigma_1,...,\Sigma_G)$ are assumed to be null hastie_elements_2009. If the density of interest is the conditional mixture density $f(y_i|x_i;\xi)$ where $x_i$ is multivariate (and may be composed of continuous and/or discrete elements), then it is necessary to assume that Assumption (ref)$(iv)$-$(v)$ also applies to the vector of covariates $x_i$. In this case, the discriminant function of the joint density classifier is written as $h_{j}(y_{i},x_i|\theta_j,\psi_j) = f_j(y_i|x_i;\theta_j)p_j(x_i|\psi_j)$, where $p_j(\cdot|\psi_j)$ refers to the covariates' density for the $j^{th}$ component, and where $\psi = (\psi_1,...,\psi_G) \in \Psi$ represents the set of parameters that governs each $p_j(\cdot)$ density. The Euclidean distance classifier is a nonparametric classifier that relies on the first moment of $y_i$ and is the one that is normally used with the K-means algorithm. The Mahalanobis distance classifier is also nonparametric but relies on the first and second moments of $y_i$ to classify each observation. Several other classifiers exist and are used in practice, but those three classifiers are among the most prevalent in the literature on unsupervised classification and clustering techniques. The next two lemmas show that all three classifiers are unbiased under specific conditions.

Lemma (ref) shows that all three classifiers are unbiased if the vector of true mean values $\boldsymbol{\mu}^0_g$ is different across all groups, whereas Lemma (ref) shows that the Mahalanobis distance classifier can be unbiased even if the vectors of true mean values are identical across components. This makes this classifier more robust to cases where the component densities feature strong overlap between each other, compared to the Euclidean distance classifier. Even though the joint density classifier is unbiased if $\boldsymbol{\mu}^0_j = \boldsymbol{\mu}^0_g \Leftrightarrow j = g$, it is generally not possible to show that it is (uniformly) consistent without specifying the distribution of the joint density. Nonetheless, Theorem (ref) and Theorem (ref) show that the Mahalanobis distance classifier is uniformly consistent if the size of $y_i$ (or $x_i$) grows at a strictly faster rate than the number of observations for a fixed value of $G$.

The conclusion of Theorem (ref) does not necessarily imply that $p > N$ is necessary to correctly classify every observation in the sample; this will however be required as $N$ tends to infinity. In this case, each estimated matrix $\hat{\Sigma}_g$ will not be invertible, and it will not be possible to compute $z^{M}_{ig}(\hat{\boldsymbol{\mu}},\hat{\Sigma})$ for any $g \in \mathbb{G}$. However, Corollary (ref) shows that the misclassification rate given by the Mahalanobis distance classifier will go to zero as $p \to \infty$ if the number of covariates grows at a strictly faster rate than the number of groups, regardless of the rate at which the sample size increases. This means that $p$ could be much smaller than $N$ in finite samples, and the Mahalanobis distance classifier would still be able to correctly classify every observation in the sample under Assumption (ref)$(ii)$-$(iii)$ and Assumption (ref). When $N < \infty$ is fixed, this corollary also indicates that $p$ has to be larger than $G$ to ensure that the misclassification rate $\hat{E}_{Np}(\theta^0)$ remains close to or equal to zero.

Remark 3. As emphasized by dzemski_convergence_2021, uniform consistency of a classifier is not necessary to obtain consistent estimates of $\theta^0$ when estimating latent group panel structures or finite mixtures. Such a statement is confirmed below in Theorem (ref) where consistency of the classifier is sufficient to yield consistent estimates of $\theta^0$. Corollary (ref) shows that the estimates obtained from the maximization of the C-ML objective, when equipped with the Mahalanobis distance classifier, will be identical to those obtained from the maximization of the complete-data log likelihood $l^{C}(\theta,\pi,\mathbf{z}^0)$ when the size of $y_i$ (or $x_i$) tends to infinity. This corollary is a direct consequence of Theorem (ref) when $p \to \infty$ and $N$ is fixed, which makes the two objective functions $l^{CM}(\theta,\mathbf{z}^M(\boldsymbol{\mu}^0,\Sigma^0))$ and $l^{C}(\xi,\mathbf{z}^0)$ identical up to a constant factor, where $\mathbf{z}^M(\boldsymbol{\mu}^0,\Sigma^0) = (\mathbf{z}^M_1(\boldsymbol{\mu}^0,\Sigma^0),...,\mathbf{z}^M_G(\boldsymbol{\mu}^0,\Sigma^0))$, with $\mathbf{z}^M_g(\boldsymbol{\mu}^0,\Sigma^0) = (z^M_{1g}(\boldsymbol{\mu}^0,\Sigma^0),...,z^M_{Ng}(\boldsymbol{\mu}^0,\Sigma^0))^{\top}$. As a result, the finite-sample bias of the estimated parameters will decrease linearly with $N$, as long as $p/N \to \infty$ such that no observation is misclassified.

In practice, small values of $p$ will lead to larger biases in both $\hat{\theta}^{CM}(\mathbf{z}^M(\cdot,\cdot))$ and $\hat{\pi}_g^{CM}$ compared to their infeasible counterparts $(\hat{\theta}^{C},\hat{\pi}^{C})$. However, results from the simulations presented in Section (ref) tend to show that the size of $p$ does not have to be very large to yield relatively small finite-sample biases and low variances in the estimated parameters when Assumption (ref) is satisfied. The next theorem and corollary show that the Euclidean distance classifier does not offer the same guarantee as the Mahalanobis distance classifier under Assumption (ref). In fact, the Euclidean distance classifier may misclassify a nonvanishing proportion of the observations in the sample if the distance between each vector $\boldsymbol{\mu}^0_j$ does not become sufficiently large as $p \to \infty$. This is because the Euclidean distance is based only on the mean values, which is not sufficient to guarantee that any single observation in the sample is correctly classified as $p \to \infty$ under Assumption (ref). For instance, if a fixed number of true mean values differ across groups, Theorem (ref) and Corollary (ref) show that the ratio $\frac{||\boldsymbol{\mu}^0_g - \boldsymbol{\mu}^0_{j}||^2}{\text{tr}(\Sigma^0_{g})}$ will converge to 0 as $p \to \infty$, which will make the upper bound on $\mathbb{P}[\cup_{g=1}^G (z^{E}_{i(g)}(\boldsymbol{\mu}^0) \ne z^0_{ig})]$ converge to one as $p \to \infty$. This upper bound may still be sharpened under the assumption that each observation is randomly assigned to any given group when all true mean values are identical across groups.\footnote{For simplicity, I assumed that identical distances across groups necessarily imply misclassification, although this depends on the chosen tie-breaking rule. Note also that identical true mean values do not violate Assumption (ref)$(ii)$ as long as Assumption (ref)$(v)$ is satisfied.} On the contrary, the Mahalanobis distance classifier is robust to identical true mean values across groups, and is a consistent classifier as long as a sufficiently large number of elements in the covariance matrix $\Sigma^0_g$ differ across groups.

Remark 4. When all elements in $y_{i}$ are normally distributed, the joint density classifier can be expressed as a function of the Mahalanobis distance given that the squared Mahalanobis distance is embedded into the exponential part of the multivariate normal distribution. In this case, the joint density classifier is both unbiased and consistent under the same conditions as the ones stated in Lemma (ref) and Corollary (ref). Simulations confirm that the joint density classifier leads to a smaller misclassification rate than the Mahalanobis distance classifier in finite samples when the normality assumption is verified (results not shown).

If the normality assumption is not verified, the Mahalanobis distance classifier offers a very good nonparametric alternative to the joint density classifier. Extending the use of the Mahalanobis distance classifier to a fully nonparametric estimation strategy is straightforward via the use of a kernel density estimator for estimating each $f_g(y_{i}|\theta_g)$ at each maximization step of the C-EM algorithm (see Algorithm (ref)). Note that the proposed C-EM algorithm is more general than the CEM algorithm of celeux_classification_1992 since it leaves the classifier $z_{ig}(\theta)$ unspecified, as long as it is consistent. Note also that the C-EM algorithm can be adapted to any conditional component density of the form $f_g(y_i|x_i;\theta_g)$, where the chosen classifier may be based on $y_i$ and/or $x_i$. For instance, the Mahalanobis distance classifier could be applied to the vector of covariates $x_i$, where the maximization step would have to also include the estimation of $\hat{\psi}^{(k+1)} = \arg \max_{\psi \in \Psi} \sum_{i=1}^N \sum_{g=1}^G z^M_{ig}(\hat{\psi}^{(k)}) \log(p_g(x_{i}|\psi))$. I now derive the asymptotic distribution of the estimates given by Algorithm (ref) under the following assumptions.

\setstretch{1.3} Assumptions (ref)$(i)$-$(ii)$ are necessary to ensure that the algorithm leads to consistent estimates of $\theta$. Those assumptions can be adapted to the case of asymptotically biased estimators, such as MLE applied to nonlinear fixed effects models, by using a bias-corrected estimator. Assumption (ref)$(iii)$ implies that the number of observations within each group goes to infinity as the sample size increases. Assumption (ref)$(iv)$ is a special case of Assumption (ref)$(iii)$ and is standard in most MLE problems. Finally, Assumptions (ref)$(v)$ states that the true Hessian matrix $\mathcal{H}_g$ is finite and non-singular for each $g \in \mathbb{G}$ and that the Hessian matrix $\mathcal{H}(\theta_g)$ can be computed by using the derivative of the score vector with respect to $\theta_g$ for each group separately. The next lemma and two theorems follow from Assumption (ref) and Algorithm (ref), and can be adapted to panel data and time series by modifying the score function, the Information matrix, and the Hessian matrix appropriately. Lemma (ref) shows that the C-ML function never decreases between two consecutive iterations of the C-EM algorithm when the joint density classifier is employed. Theorem (ref) shows that the C-EM algorithm will yield consistent estimates of all parameters in the mixture after a finite number of iterations provided appropriate initial parameter values, and a sufficiently large size of $y_i$ (or $x_i$) and sample size. If those conditions are satisfied, Theorem (ref) shows that each estimated component parameter $\hat{\theta}^{(k)}_g$ will be normally distributed asymptotically.

Remark 5. It is worth noting that Theorem (ref) refers to the asymptotic behavior of $\hat{\theta}^{(k)}$ conditional on the previous set of estimates $(\theta^{(0)},\hat{\theta}^{(1)},...,\hat{\theta}^{(k-1)})$. This can be very different from the asymptotic behavior of the whole sequence of estimated parameters, which implies that $\hat{\theta}^{(1)}$ is such that $\mathbf{z}({\theta}^{(0)}) = \mathbf{z}(\hat{\theta}^{(1)})$ asymptotically given a consistent classifier and initial parameter values that are sufficiently close to the true values. On the other hand, it may be that $\hat{\theta}^{(k)}$ does not converge in probability to the true parameter values if we fix the $k-1$ sets of previous estimates, even if the algorithm is initialized at the true parameter value. This will happen if the entire sequence $(\theta^{(0)},\hat{\theta}^{(1)},...,\hat{\theta}^{(k)})$ converges, in finite samples, to a global maximum of the C-ML function that does not coincide with the global maximum asymptotically. This is why Theorem (ref) states that $N$ and $p$ have to be sufficiently large to guarantee the consistency of $\hat{\theta}^{(k)}$ due to the conditional nature of the estimated parameters.

Remark 6. The convergence of $\hat{\theta}^{(k)}$ is not uniform over $\pi^0$ since arbitrarily small values of $\pi^0_g$ can make the convergence of $\hat{\theta}_g^{(k)}$ arbitrarily slow as well. This is why all convergence results presented here are “group-wise” under Assumption (ref)(iii). It is possible to make $\hat{\theta}^{(k)}$ converge to $\theta^0$ uniformly over $\pi^0_g$ by modifying Assumption (ref)(iii) such that $\lim_{N \to \infty} N^0_g/N = \alpha_g$, where $\alpha_g$ is bounded away from zero for every $g \in \mathbb{G}$, which is a stronger assumption than the original one and might not be appropriate depending on the context. Note also that the estimates provided by the C-EM algorithm will asymptotically reach the semiparametric efficiency bound even if the employed classifier is not uniformly consistent since $\hat{\theta}^{(k)}$ will converge to the same distribution regardless of the employed classifier, as long as it is consistent. As in standard MLE problems, misspecification of the component densities will make the C-EM algorithm converge to a pseudo-true value $\bar{\theta} \ne \theta^0$ that minimizes the Kullback-Leibler divergence between the “working” and the true component densities gourieroux_pseudo_1984. Misspecification of the component densities might also lead to the inconsistency of the joint density classifier if the working density is too far from the truth. On the other hand, the Mahalanobis distance classifier remains consistent under Assumption (ref)(ii)-(iii) and Assumption (ref) even if the component densities are misspecified. This is one of the main benefits of using the Mahalanobis distance classifier compared to the joint density classifier. \\ Remark 7. Under misspecification of the component densities and/or lack of efficiency, estimation of $\text{Var}[\hat{\theta}^{(k)}_{g}]$ can be done using the sample counterpart of the asymptotic variance given in Theorem (ref). This estimated variance is defined as follows

where $s'(\hat{\theta}_g^{(k)}) = \sum_{i=1}^N \frac{\partial}{\partial \hat{\theta}^{(k)}_{g}} s_{i}(\hat{\theta}_g^{(k)})$, and where $\hat{\theta}^{(k)}$ is the set of parameters that maximize the C-ML function after searching for the global maximum of the objective function. Note that the expressions for both the asymptotic variance and its estimated counterpart will not account for the finite-sample uncertainty in group memberships. This uncertainty can be taken into account by the use of resampling techniques, such as the bootstrap, or by constructing confidence sets that account for group memberships dzemski_confidence_2024, higgins_bootstrap_2024. Implementing such methods goes however beyond the scope of this paper. The next section confirms the theoretical intuitions and results presented in this section via simulations. Those theoretical intuitions and results are summarized as follows : First, standard MLE of finite mixtures yields estimates that can be heavily biased in finite samples, with the size of the bias increasing with the overlap between the component densities. Second, maximizing the C-ML function equipped with a consistent classifier will yield consistent (as $N,p\to \infty$) and less biased estimates in finite samples compared to estimates obtained from standard MLE. Finally, the C-EM algorithm presented above can be used to maximize the C-ML function, which will yield consistent estimates provided that the algorithm is initialized with parameter values that are sufficiently close to their true counterparts.

Monte Carlo Simulations

Data-generating Processes and Estimation Strategies

In this section, I describe the objectives and details of two distinct simulation exercises. The main objective of the first simulation exercise is to confirm the presence of a large finite-sample bias in the MLE of finite mixtures for relatively small values of $N$. The main objective of the second simulation exercise is to compare the finite-sample performance of the EM and the C-EM algorithms in the context of a mixture of linear panel data, where each algorithm is used to maximize the mixture likelihood and the C-ML function, respectively.

Finite-Sample Bias of Standard MLE

This first simulation exercise uses two-component mixtures of various distributions to show that maximizing the mixture likelihood function $l^{ML}(\xi)$ may yield heavily biased estimates of all parameters in the mixture. All datasets are generated using the following DGP :

where $y_i$ is univariate and $f(\cdot)$ represents the normal, Poisson, or exponential distribution, and where each $z^0_{i1} = \{0,1\}$ is set arbitrarily such that the set of true mixing weights $\pi^0 = (\pi^0_1, 1-\pi_1^0)$ corresponds to $(N^{-1}\sum_{i=1} z^0_{i1},N^{-1}\sum_{i=1} (1-z^0_{i1}))$. To ensure the consistency of the obtained estimates in the case of the mixture of normals, a penalized version of $l^{ML}(\xi)$ is maximized with respect to $\xi$ via the EM algorithm of dempster_maximum_1977.\footnote{The EM algorithm is similar to the C-EM algorithm shown in Algorithm (ref) but with $z_{ig}(\theta)$ replaced by $\tau_{ig}(\theta)$ in both steps, the latter corresponding to the probability that the $i^{th}$ observation belongs to the $g^{th}$ group. There are several ways to compute this probability, but the most popular and the one that I use in this paper is the following : $\tau_{ig}(\theta) = \frac{\pi_gf_g(y_i|\theta)}{\sum_{g=1}^G\pi_gf_g(y_i|\theta)}$, where $\pi_g$ is estimated by averaging $\tau_{ig}(\theta)$ (from the previous iteration) over $i$.} The penalized mixture likelihood for the mixture of normal distributions is described as follows

where

which corresponds to the recommendation of chen_statistical_2023 for mixtures of normal distributions in large samples. The other simulations based on mixtures of Poisson and exponential distributions maximize the standard mixture likelihood $l^{ML}(\xi)$ since they do not suffer from the same consistency problem as with the mixture of normal distributions. Finally, to find the global maximum of ${l}^{ML}(\xi)$ (or $\tilde{l}^{ML}(\xi)$ for the mixture of normals), I use the following initialization and maximization strategy for each simulated dataset :

The goal of this strategy is to compare the results provided by the maximization of ${l}^{ML}(\xi)$ (or $\tilde{l}^{ML}(\xi))$ via the EM algorithm for initial parameter values that are either equal to the true values, or equal to values that are influenced by the presence of natural outliers in the dataset. Since it is well-known that the EM algorithm does not decrease the value of the mixture likelihood between two consecutive iterations, initializing the algorithm with the true parameter values $\xi^0$ should yield estimated parameters that are close to these true values unless other parameter values can yield a higher likelihood value. According to the discussion in Section (ref), the parameter values that might yield a higher likelihood value are the ones where one component density is strongly influenced by a small set of outliers while the other component density is strongly influenced by every other observation. Dividing the sample into two sets of observations such that the deemed outliers (of the same sign) are all in the same set leads to the computation of initial values that are likely to identify such parameters, while also yielding a relatively high value of $\tilde{l}^{ML}(\xi)$. To identify the bias related to the maximization of the mixture likelihood, 1,000 replications with different simulated data were performed for each type of mixture and each set of true parameters. Label switching across replications was also avoided by sorting the component densities by their estimated mean values (or estimated standard deviations when $\mu_1^0 = \mu_2^0$) for each replication. The sample size was also allowed to vary between $N=100$ and $N=10,000$ for each set of true parameter values to confirm the idea according to which there exists a value $N^*$ for each simulation scenario such that the total estimation bias and mean squared error of $\hat{\xi}^{ML}$ drop when $N > N^*$. Convergence of the EM algorithm for each scenario was assumed to be reached when the value of the log likelihood function was less than 1e-10 between two consecutive iterations, or when the number of iterations reaches 100.

Latent Group Linear Panel Structure

I consider here the typical case in latent panel structure where the conditional density of the outcome value follows a normal distribution for all values of $g \in \mathbb{G}$. The goal of this simulation exercise is to determine which one of the two objective functions, $l^{ML}(\xi)$ or $l^{CM}(\theta,\mathbf{z})$, yields the best finite-sample results when these functions are maximized by the EM and the C-EM algorithms, respectively. In this simulation exercise, $p$ refers to the number of time-varying covariates that are used for classification jointly with the conditional density of the univariate outcome. The exact DGP used to simulate the outcome variable is described as follows

where $z^0_{it}$ represents the true group membership of the $it^{th}$ observation and is generated by a categorical distribution where the vector of probabilities follows a Markovian process,\footnote{The transition matrix that is used to create such a process is generated via a Dirichlet distribution that varies across replications. This Markovian process is coincidental and is not central to further developments.} $\tilde{\beta}_{g} = (\beta_{g}^{\top}, \gamma_{g}^{\top},\delta_{1g},...,\delta_{Tg})^{\top}$, $X_{it,1} = (x_{it,1}, \bar{x}_{i,1},\mathbbm{1}[t=1],....,\mathbbm{1}[t=T])$, $x_{it,1}$ is a time-varying covariate and stands as the first element in the vector of covariates $x_{it}$, $\bar{x}_{i,1}^{\top} = T^{-1} \sum_{t=1}^T x_{it,1}$ is the time-averaged value of $x_{it,1}$, $\beta_{g}$ and $\gamma_{g}$ are both scalar parameters to be estimated for each value of $g \in \mathbb{G}$, $\alpha_{ig} \sim N(0,\sigma^2_{\alpha,g})$ where $\sigma^2_{\alpha,g} = g$, and where $\delta_{tg}$ is a time-fixed effect that is normally distributed with mean $\bar{x}_{tg,1} := \frac{\sum_{i=1}^N\mathbbm{1}[z^0_{it}=g]x_{it,1}}{\sum_{i=1}^N\mathbbm{1}[z^0_{it}=g]}$ and unit variance. For simplicity, I specify $\epsilon_{it} \sim N(0,\sigma^2_{\epsilon})$ with $\sigma^2_{\epsilon} = 1$, and also $\mathbb{E}[\epsilon_{it}\epsilon_{ls}] = 0$ for any pair $(l,s) \ne (i,t)$. The vector of covariate $x_{it}$ is generated according to

where $\mathcal{N}_p$ denotes the $p$-variate normal distribution, $\boldsymbol{\mu}^0_{g}$ is the vector of true mean values for the $g^{th}$ group where the $j^{th}$ element $\mu^0_{gj}$ is generated by a standard normal distribution for every $j \in \{1,2...,p\}$, and where $\Sigma^0_{g}$ is generated according to the following equation

with

where $\omega_{g,ij} \sim N(0,1)$ for all triples $(g,i,j) \in \mathbb{G}\times \{1,2...,p\}^2$. Both $\boldsymbol{\mu}^0_{g}$ and $\Sigma^0_{g}$ vary across replications for every $g \in \mathbb{G}$. Defining the matrix $P_g$ as such guarantees that $\Sigma^0_{g}$ will always be positive-definite for each $g \in \mathbb{G}$. Finally, the parameters $\beta_g$ and $\gamma_g$ are drawn from a standard normal distribution for all values of $g \in \mathbb{G}$. Note that only the first element of $x_{it}$ is included as a regressor in Eq. ((ref)) for simplicity. Including the other generated covariates in the specification of Eq. ((ref)) does not alter significantly the conclusions of the simulation exercise, but makes the results less interpretable due to a larger number of elements in $\beta_g$. Eq. ((ref)) corresponds to the Mundlak specification where the time-average values of the covariates plus unit-random effects are used in place of unit-specific fixed effects. It has been shown that it is equivalent to the least-square dummy variable (LSDV) estimator wooldridge_two-way_2021 but with the advantage of being more computationally tractable since it greatly reduces the number of parameters to estimate. Moreover, using the Mundlak approach allows the non-random part of the unit-fixed effects, $\bar{x}_{i,1}^{\top}\gamma_{z^0_{it}}$, to vary across groups for any given individual. The LSDV estimator also allows for such flexibility, but each unit has to be assigned to the same group for at least two periods in order to avoid identification issues when using the C-EM algorithm. The Mundlak specification also facilitates cross-validation since the non-random part of the unit-fixed effect can be predicted for any unit in the test set. This framework also avoids the biases generated by the incidental parameter problem encountered in nonlinear panels with unit-fixed effects and fixed $T$. As mentioned earlier, estimation is carried out using both the EM and the C-EM algorithms for comparison. The C-EM algorithm uses the joint density classifier with normal distributions where the discriminant function $h_{j}(y_{it},x_{it}|\theta_j,\psi_j) = f_{j}(y_{it}|x_{it};\theta_j)p_j(x_{it}|\psi_j)$. To make sure that the C-ML function never decreases between two consecutive iterations of the C-EM algorithm, I use $l^{CM}(\theta,\psi,\mathbf{z}) = \sum_{i=1}^N\sum_{t=1}^T\sum_{g=1}^G z^D_{itg}(\theta,\psi) \log(f_{j}(y_{it}|x_{it};\theta_j)p_j(x_{it}|\psi_j))$ as the objective function (see the proof of Lemma (ref) for a justification). All models are correctly specified, including the use of normal distributions for the random-unit effect $\alpha_{ij}$, the error term $\epsilon_{it}$, and the density of $x_{it}$. To ease computational burden, 250 replications were performed with only 25 distinct sets of randomly-generated initial parameter values per replication. Comparison of the finite-sample performance between the EM and the C-EM algorithms was performed by using the same set of initial values for both algorithms. All covariates, true parameters, and random terms are allowed to vary across replications according to their respective distribution or generating process. More details regarding the two estimation procedures for this second simulation exercise are given in Appendix (ref). Since the random part $\alpha_{iz^0_{it}} + \epsilon_{it}$ is normally distributed, the penalized likelihood function $\tilde{l}^{ML}(\xi)$ is used instead of the regular mixture likelihood $l^{ML}(\xi)$ to ensure consistency. Non-singularity of the estimated $\hat{\Sigma}^{(k)}_g$ is obtained by ensuring that all eigenvalues of $\hat{\Sigma}^{(k)}_g$ remain above 1e-8 for every $g \in \mathbb{G}$ at each iteration. The maximum number of iterations for every replication is set to 100. Convergence of the EM algorithm is assumed to be reached when the relative change between two consecutive values of $\tilde{l}^{ML}(\xi)$ is less than 0.01%. Convergence of the C-EM algorithm is achieved when the estimated group memberships do not change between two consecutive iterations. All simulations are performed using Python 3.11.10 with NumPy version 2.2.0. Replication codes for the two simulation exercises are available on GitHub \href{https://github.com/raph1177/Finite_mixture}{here}. Additional simulation results for mixtures of binary choice models (Probit) are also available from the author upon request.

Simulation Results

Finite-Sample Bias of Standard MLE

\setstretch{1.3} Figure (ref) shows the results of the first simulation exercise for the two-component mixture of normal distributions with true mean values $\boldsymbol{\mu}^0 = (0.25,-0.25)$, true standard deviations $\sigma^0 = (1,1)$, and true mixing weight $\pi^0_1 = 0.5$. Panel (a) of Figure (ref) shows that the estimation bias for both the mean values and the standard deviations, which corresponds to the distance between the blue dots and the red dashed line, tends to increase with the sample size. Note that the estimated mixing weights are unbiased in this case, but the estimation variance (due to the unbiasedness) is very high because some estimated mixing weights are very close to the boundary of the parameter space $\Pi$. Note also that the upper graphs of Figure (ref) confirm the formal intuition described in Example (ref), which is that $\hat{\mu}_1$ and $\hat{\mu}_2$ are positively and negatively biased in finite samples, respectively, when the distance between the true mean values is small and ${\mu}^0_1 > {\mu}^0_2$. The results shown in panel (b) of Figure (ref) also confirm that, for several simulated samples, one of the two sets of estimated parameters that maximize $\tilde{l}^{ML}(\xi)$ is strongly influenced by a small number of outliers. Indeed, panel (b) shows that, when $N=10,000$, several replications feature a set of maximizers $\hat{\xi}$ that are heavily influenced by outliers, as indicated by the strange estimated parameter values across several replications. Only a few number of estimated parameters feature a very small level of bias compared to the others, indicating that not all simulated datasets generated sufficiently large outliers (or a sufficiently large number of outliers) to drastically change the location of the global maximum of $\tilde{l}^{ML}(\xi)$.

\setstretch{1.3} Figure (ref) shows the results for the two-component mixture of normal distributions with true parameter values that are similar to those used to produce Figure (ref), except for the mean values where $\boldsymbol{\mu}^0 = (0.75,-0.75)$. Increasing the distance between the true mean values leads to a reduction in the degree of overlap between the two component densities, which translates into a smaller influence of outliers on the estimated parameters. Figure (ref) also shows that increasing the sample size from 1,000 to 10,000 leads to a sharp reduction in estimation biases and/or mean squared errors for all parameters. Specifically, the results portrayed in Figure (ref) confirm the intuition behind Figure (ref), which is that there exists a threshold value $N^*$ such that when $N> N^*$, the location of $\hat{\xi}^{ML}$ then becomes much closer to the true parameter values $\xi^0$. In the case illustrated by Figure (ref), it is clear that this threshold lies somewhere between $1,000$ and $10,000$ for most simulated datasets. Determining the location of this threshold as a function of $\xi^0$ and characteristics of the data goes beyond the scope of this paper but is an interesting avenue for future research. Other simulation results with a mixture of normal distributions where standard deviations and mixing weights are allowed to vary are shown in Appendix (ref). For instance, Figure(ref) shows that the value of the threshold $N^*$ is, on average, lower when the true standard deviations differ across components, and that this value may differ across groups. Indeed, Figure (ref) shows that the reductions in estimation biases and mean squared errors are captured primarily by the first group of parameters given the larger value of the true mixing weight $\pi_1^0$. However, the differences in standard deviations and mixing weights across groups tend to matter less as the distance between the true mean values increases, as shown in Figure (ref). Figure (ref) also shows that large differences in standard deviations alone are sufficient to reduce estimation biases and mean squared errors with a moderately large sample size. In fact, this is the only scenario where no component density is ever influenced by any outlier when $N = 10,000$ (see Panel (b) of Figure (ref)). Finally, Figures (ref) to (ref) of Appendix (ref) show similar results to the previous figures, but for mixtures of Poisson and exponential distributions. Since these distributions are both governed by a single parameter, the estimation biases and mean squared errors are less salient than for mixtures of normal distributions. Nonetheless, Figures (ref) and (ref) show that both the estimated mean and mixing weight parameters from a mixture of Poisson or exponential distributions can be heavily biased in small- to medium-sized samples when the true mean values are too close to each other.

\setstretch{1.3}

Latent Group Linear Panel Structure

The four graphs of Figure (ref) show the estimation biases of the mean parameters (i.e., $\beta_g$, $\gamma_g$, and $\delta_{tg}$ for all $g \in \mathbb{G}$ and $t \in \{1,2,...,T\}$), variance parameters, and mixing weights for each group and each algorithm when $N=500$, $T=5$, $G=2$, and the number of variables $p$ is equal to one. This implies that the joint density classifier is based only on two variables : the outcome $y_{it}$ (conditional on the linear predictor $X_{it,1}\tilde{\beta}_g$) and the scalar covariate $x_{it,1}$. The two upper graphs of Figure (ref) show that the C-EM algorithm yields mean parameter estimates that have similar biases -- except for $\beta_g$ -- but much smaller mean squared errors compared to those produced by the EM algorithm. The two other graphs of Figure (ref) also show that the variance terms and the mixing weights estimated by the C-EM algorithm feature less variability and smaller or equal bias than those produced by the EM algorithm. This is because the C-EM algorithm, even with $p$ as low as 1, leads to a misclassification rate that is, on average, lower than the one produced by the EM algorithm. Note that the misclassification rate for the EM algorithm is based on the maximum posterior probability at convergence, and it should be equal to or close to zero in order to limit estimation biases and variability.

\setstretch{1.3} Figures (ref) and (ref) show similar results, but when the number of covariates $p$ is equal to 5 and 10, respectively. When $p$ goes from 1 to 5, there is a substantial reduction in mean squared errors for all parameters estimated by the C-EM algorithm, including the variance terms and the mixing weights. The misclassification rate is also close to zero for the vast majority ($\ge 95\%$) of the replications. As expected, the overall performance of the EM algorithm is not influenced by the additional covariates that are not included in the outcome model shown in Eq. ((ref)). It is worth noting that introducing additional covariates in the outcome model will not improve the general performance of the EM algorithm unless the associated coefficients largely differ across groups, but this will also improve the performance of the C-EM algorithm in finite samples. This intuition is confirmed by additional simulation results that are available from the author upon request. When the number of covariates is equal to 10, the results illustrated in Figure (ref) show that the C-EM performs even better than with $p=5$. However, the marginal reduction in biases and mean squared errors diminishes as the number of covariates increases, which is due to the C-EM reaching the oracle efficiency lower bound for this sample size. The lower graph of Figure (ref) indeed shows that the misclassification rate produced by the C-EM algorithm is equal to zero in (at least) 19 replications out of 20. In this case, the C-EM algorithm produces mean parameter and mixing weights estimates with very low biases and mean squared errors, while the estimates produced by the EM algorithm remain slightly biased but highly variable. Note also that the bias and precision of the variance terms estimated by the C-EM algorithm are similar between $p=5$ and $p=10$ due to the short length of the panel.

\setstretch{1.3} Figures (ref) to (ref) in Appendix (ref) show the set of results obtained from the same simulation scenarios as the ones used to produce Figures (ref) to (ref), but with $N=750$ and $T=8$. When $p=1$, Figure (ref) shows that increasing the sample size does not necessarily lead to less biased and/or less variable estimates for either algorithm. However, when $p$ increases from 1 to 10, the C-EM algorithm produces less biased estimates with lower mean squared errors than with the smaller sample size, which is true for all parameters, including the variance terms. Figures (ref) to (ref) also show that there is no clear improvement in the performance of the EM algorithm after increasing the sample size from $NT = 2,500$ to $NT = 6,000$. Such a poor finite-sample performance can be explained by an insufficient increase in the sample size, or perhaps by the relatively small number of initial parameter values that were used to identify the global maximum of each objective function. In any case, the results obtained with the C-EM algorithm demonstrate that it is not always necessary to use a larger number of initial parameter values to yield estimated parameters that are close to true values. Finally, Figures (ref) to (ref) in Appendix (ref) show the set of results associated with similar simulation scenarios as the three shown above, but with an additional group (i.e., $N=500$, $T=5$, and $G=3$). For each value of $p$, increasing the number of groups leads to higher misclassification rates for the C-EM algorithm compared to when $G=2$, which then leads to higher mean squared errors for most estimated parameters. These results also show that increasing the number of groups $G$ entails similar effects on the parameters estimated by the C-EM algorithm as if the number of covariates $p$ had been reduced with $G$ fixed. This validates the conclusion of Corollary (ref), which states that the ratio $p/G$ has to go to infinity for the classification procedure to be consistent. Therefore, increasing the number of groups $G$ while keeping $p$ constant increases the misclassification rate and thus leads to larger biases and/or mean squared errors overall. If the general form of the joint density $f(y_{it},x_{it}|\theta,\psi)$ is known, the results presented in this section show that the joint density classifier generally outperforms the results produced by the EM algorithm Assumption (ref) is satisfied, even when the numbers of covariates employed to improve classification is as small as 1. To note, the EM algorithm produces values of $\hat{\beta}_g$ that are always positively biased across all simulation scenarios, while the C-EM algorithm always yields mean parameters with very low biases and mean squared errors when $p \ge 5$, although this may be an artifact of the employed DGPs. If the general form of the joint density is unknown and no assumption on this general form is being made, then the Mahalanobis distance classifier can be used as a nonparametric substitute for every classification step of the C-EM algorithm. However, the Mahalanobis distance classifier typically leads to higher misclassification rates compared to the joint density classifier (when the parametric assumptions are correct) and might result in a decrease of the C-ML function between two consecutive iterations. Using the Mahalanobis distance classifier with a slightly larger number of covariates will however lead to similar performances as with the joint density classifier in the context of the current simulation exercise (results not shown).

Empirical Application

Objectives

The main objective of the empirical analysis is to group within the same component all the observations that share the same latent, unobserved individual characteristics. In this context, it is assumed that the unobserved characteristics explain a non-null proportion of the observed variation in individual healthcare expenditures (HCE), and that the unobserved characteristics can also influence the relationships between the outcome and the observed individual characteristics. For instance, individuals who are more frail than others, which is unobserved, are more likely to develop adverse health outcomes and become more dependent on the healthcare system after experiencing a minor injury. Such poor underlying characteristics will be reflected in the data by higher individual HCE in the periods following the minor injury. This approach is similar to health state modeling in the literature on hidden Markov models (HMMs) luo_bayesian_2021,komariah_health_2019 but using a two-step approach where the first consistently estimates the group membership (i.e., the “health state”) and the second step models the dynamic behavior of the group membership. Unlike this literature, I will refer to health groups rather than “health state” given that the unobserved individual characteristics explaining HCE can be related to non-health factors (e.g. accessibility to healthcare services, peer effects, etc.). The second objective of the empirical analysis is to compare the results coming from both algorithms (C-EM and EM) when using real-world data with missing information. The predictive performance of the results obtained from each algorithm is compared using cross-validation, as detailed in Section (ref).

Data Sources and Characteristics

The employed dataset is the Quebec portion of the Canadian Emergency departments Team Initiative (CETI). The CETI research program on mobility and aging is a Canadian clinical program “which aims to improve emergency department (ED) care for older adults with minor injuries” provencher_decline_2015. In the province of Quebec, the CETI research program followed 1,391 patients with a medical consultation at the ED of one of three hospitals (Hôpital du Sacré-Coeur de Montréal, Hôpital de l’Enfant-Jésus, and Hôpital du Saint-Sacrement) after a minor injury that occurred between 2009 and 2015. Individuals were included in the cohort if they were aged 65 and older and suffered from a minor injury that did not lead to hospitalization. A detailed description of the eligibility criteria and monitoring schedules can be found in provencher_decline_2015. The health administrative data for each participant has been extracted from Quebec's physician claims database, which is managed by the Régie de l’Assurance-Maladie du Quebec (RAMQ). Each claim included a unique patient identification number, billing codes corresponding to the services rendered, the diagnostic code according to the International Classification of Diseases, 9th Revision (ICD-9-CM), dates and locations of the services provided, an identification code for the physician making the claim, and the amount paid to the physician. The period covered by the dataset goes from January 1, 2008, to June 30, 2016. Note that physician claims do not include all the HCE that an individual could incur in a given period. However, physician claims are strongly correlated with total HCE while being usually more precise than other cost measures (e.g., hospitalization costs, costs of treatment episodes, etc.). An unbalanced panel dataset containing 1,330 individuals and seven three-month periods per individual was created using the information described above. Some individuals died before the end of the last period, hence leading to 64 observations with no information. A weight of zero was attributed to those 64 observations at every classification/expectation step of both algorithms. Consequently, the total number of observations with non-null weights in the dataset is 9,246. The initial ED visit occurred at the end of the second period, while two follow-up visits occurred at the end of the third and fourth period,s respectively. Price effects (inflation) were accounted for by adjusting all costs with Quebec’s all-item consumer price index (CPI) between 2008 and 2016.

Estimation Strategy

The employed model is a mixture of two-part models, where the first part corresponds to a Probit binary choice model and the second part uses a lognormal distribution. The lognormal distribution accounts for the heavy-tailed distribution that is usually observed in individual HCE manning_estimating_2001. Conceptually speaking, the first part models the “decision” of the $i^{th}$ participant to incur a strictly positive amount of HCE at the $t^{th}$ period, whereas the second part models the total amount, if any, that was incurred at this period by the participant. Two-part models have been extensively used in econometrics and health economics to model the decision-making of agents that lead to the generation of a strictly positive, continuous outcome depending on the initial decision of the agent norton_specification_2008, neelon_bayesian_2011. The specification used within each part of the model is similar to the one used for the second simulation exercise, as formulated by Eq. ((ref)). For simplicity, it is assumed that all covariates follow a normal distribution with means and covariance matrices varying across components. The general form of the $g^{th}$ joint component's density for the $it^{th}$ observation is written as follows

with $\eta^b_{itg} = X_{it}^{b}\tilde{\beta}^{b}_{g}$ and $\eta_{itg} = \frac{X_{it}\tilde{\beta}_g-y^c_{it}}{\sigma^{2}_{\alpha+\epsilon,g}}$, where $f_N(\cdot|\boldsymbol{\mu}_g,\Sigma_g)$ corresponds to the multivariate normal pdf with mean $\boldsymbol{\mu}_g$ and variance-covariance matrix $\Sigma_g$, $\Phi(\cdot)$ and $\phi(\cdot)$ are defined as in Example (ref), $y^c_{it} = \log(y_{it})$ stands as the log value of the HCE for the ${it}^{th}$ observation (and is set to zero if $y_{it} = 0$), $d_{it} = \mathbbm{1}[y_{it}>0]$ is equal to one if $y_{it}>0$ and zero otherwise, and where $\theta_g = (\tilde{\beta}_g,\tilde{\beta}^{b}_g, \sigma^{2}_{\alpha+\epsilon,g})$ and $\psi_g=(\boldsymbol{\mu}^b_g,\Sigma^b_g,\boldsymbol{\mu}_g,\Sigma_g)$ for any $g \in \mathbb{G}$. All other coefficients have the same interpretation as in Section (ref), with $X_{it}$ including all covariates and fixed effects, and with the $b$ superscript denoting the estimates of the binary part. Note that all covariates and parameters associated with the binary part need not be equal to those of the continuous part, even if the true distributions of the common elements in $x^b_{it}$ and $x_{it}$ are identical. More details on the covariates included in each part are provided in the next subsection. Estimation of the model is carried out using similar approaches to the ones used in the second simulation exercise. Two differences are however worth noting. First, each M-step for both binary parts (one for each algorithm) corresponds to a single Newton-Raphson step to update the estimates $\tilde{\beta}_g^b$. Second, the C-EM algorithm uses the joint density classifier $z^{D}_{itg}(\theta,\psi)$ at each C-step of the algorithm, where $h_{j}(y^c_{it},x_{it}|\theta,\psi) = f_j(y^c_{it},x_{it}|\theta_j,\psi_j)$, whereas the EM algorithm uses the posterior probabilities $\tau_{itg}(\theta,\psi)$ at each E-step of the algorithm, where $\tau_{itg}(\theta,\psi) = \frac{f_g(y^c_{it},x_{it}|\theta_g,\psi_g)}{\sum_{j=1}^{G} f_{j}(y^c_{it},x_{it}|\theta_j,\psi_j)}$. Such posterior probabilities avoid the need to choose between the binary and the continuous density at each E-step of the EM algorithm, while leading to values of $\tau_{itg}(\theta,\psi)$ that are closer to either 0 or 1 on average. Not introducing any mixing weight in the RHS of the posterior probability $\tau_{itg}(\theta,\psi)$ also makes the two algorithms more comparable to each other. Such a formulation also suggests that the covariate densities $f_N(x^b_{it}|\boldsymbol{\mu}^b_g,\Sigma^b_g)$ and $f_N(x_{it}|\boldsymbol{\mu}_g,\Sigma_g)$ are used as prior grouping information for computing $\tau_{itg}(\theta,\psi)$. For simplicity, the unit-random effects of the binary and the continuous parts are assumed to be independent of each other. It has been shown that this assumption is likely to introduce bias in the continuous part of the model su_bias_2009. Note that this bias is different from the sample selection bias often encountered in econometrics heckman_sample_1979. If the probability for every individual of generating a non-null amount of HCE is never equal to zero over time, then the sample selection bias will vanish as $T$ increases (provided a representative sample). On the other hand, the bias introduced by the independence assumption between the two parts of the model may be maintained even in the absence of any sample selection bias. However, the finite mixture setup naturally accounts for cross-part correlation that might be caused by correlated group-specific intercepts/time-fixed effects. Introducing correlated random effects is also possible using Bayesian estimation methods or simulated maximum likelihood, but using such methods goes beyond the scope of this paper. Inference on the estimated parameters is performed using a cluster-robust variance estimator that is described in Appendix (ref). If the independence assumption described above is satisfied, then consistent estimation of the parameters can be done by estimating each part of the model separately. This implies that only the observations with a strictly positive outcome value are used to estimate the continuous parts of the model. The estimation of the binary parts is performed using an iterative weighted Newton-Raphson procedure that uses all the observations in the sample, where the weights correspond to the joint density classifier (for the C-EM) and the posterior probabilities (for the EM) described above. Convergence of the whole estimation procedure is assumed to be achieved when the relative change between two consecutive log likelihood values is less than 0.01%. A maximum number of 100 iterations is enforced for each iterative procedure.

Covariates

The list of covariates included in $X^b_{it}$ and $X_{it}$ is shown in Appendix (ref). The included covariates are the same across all component densities. To proxy for frailty, I use the Elder's Risk Assessment (ERA) index, which predicts the hospitalization and health risks among elders crane_use_2010. Frailty is known to be an important predictor of medical resource consumption and individual HCE over time sirven_cost_2017. The composition of the ERA index is shown in Appendix (ref) and slightly differs from the original index due to data availability issues. In addition to frailty, the global amount of comorbidity is also known to be an important predictor of individual HCE charlson_charlson_2008. Consequently, the Charlson index has been used to control for the overall burden of comorbidities in both parts of the model. The composition of the Charlson index is shown in Appendix (ref). Given that the ERA and Charlson indices have comorbidities in common, the overlapping covariates were removed from the Charlson index to limit collinearity issues. Note that only the time-averaged Charlson index was used in every specification due to strong collinearity between the time-varying and the time-averaged Charlson indices. Finally, the Continuity of Care Index (COCI) was also introduced in the model to account for the peculiar relationships between each patient and the healthcare system. Continuity of care is defined as how one patient experiences care over time as coherent and linked bice_quantitative_1977, haggerty_experienced_2013. Values of the index range from zero to one, with a zero value referring to a total absence of continuity of care (i.e., each visit is associated with a different provider), whereas a value of one represents perfect continuity of care (i.e., all visits are associated with the same unique provider). Several studies have shown that the COCI is strongly associated with individual HCE and that even modest variations in the index value are associated with large variations in medical costs chu_continuity_2012,hussey_continuity_2014. The computation of the COCI relies on the identification of the billing physician at each visit, which cannot be done if $y_{it}=0$ for a given pair $(i,t) \in [N]\times \{1,2,...,T\}$. Although the missing COCIs could be proxied from other information in the dataset, I include the time-varying COCI only in the continuous part of the model for simplicity.\footnote{To maintain comparability with the other coefficients, the COCI has been rescaled from 0 to 10.}

Selection of the Initial Parameter Values and Number of Groups

Between 300 and 2,000 different sets of initial parameters were assessed to estimate the model for each algorithm and each value of $G$.\footnote{The number of initial parameter values is different for each algorithm since it is more frequent for the C-EM algorithm to experience convergence issues due to the presence of empty component(s). The exact number of initial parameter values used for each value of $G$ is indicated in the Python code provided by the author upon request.} The selection of the “optimal” initial parameter values and number of groups $G$ was performed using the maximum value of the corresponding likelihood function and cross-validation (CV). If some observations are misclassified, then maximizing the likelihood might select parameter estimates that may be quite far from the true parameters. This is why CV was also used to select the optimal set of initial parameter values and the number of groups. The goal of the CV procedure is to estimate the out-of-sample prediction error associated with the estimated parameters and the chosen number of groups. Following the recommendations of zhang_cross-validation_2015, a 2-fold CV procedure with 10 different data splittings was performed for each value of $G \in \{2,6\}$ and each set of estimates associated with the 15 highest likelihood values.\footnote{The 15 highest likelihood values were chosen to reduce the computational burden of the CV procedure for the selection of the initial parameter values.} All data splittings randomly allocated each unit to the training set or the test set, which maintains the correlation structure among units in the training set. To limit the probability of being “trapped” in a local maximum of the objective function, the estimated parameter values of the complete dataset were used as the initial values for each repetition during the CV procedure. The predicted values of the test set for both algorithms are computed as follows

where $w^{(k)}_{itg}$ is defined as in Section (ref), and where $\hat{\tilde{\beta}}^{b,(k)}_{g}$ and $\hat{\tilde{\beta}}^{(k)}_g$ are the respective estimated analogs of $\tilde{\beta}^{b}_g$ and $\tilde{\beta}_g$ after convergence of the algorithm. The cross-validated RMSEs were computed by combining the prediction errors of all observations for each test set and each data splitting. Note that all $\hat{\tilde{\beta}}^{(k)}_g$ correspond to semi-elasticities and that retransformation to the original scale is not necessary to perform the CV procedure.

Results

Selection of the Initial Parameter Values and Number of Groups

Figure (ref) shows the relative, cross-validated RMSE values obtained by the CV procedure described in Section (ref) for each algorithm, each value of $G$, and each set of initial parameter values that produced the 15 highest likelihood values. The values shown in Figure (ref) are relative to the RMSE computed when $G=1$, which corresponds to a value of 2.05. The results illustrated in Figure (ref) show that the EM algorithm minimizes the cross-validated RMSE when $G=4$ while the C-EM algorithm minimizes the cross-validated RMSE when $G=5$. Several repetitions of the CV procedure have been performed to ensure the stability of the lowest RMSE values for each algorithm (results not shown). The lowest RMSE value obtained from the C-EM algorithm is equal to 0.89, which is 17.6% smaller than the lowest RMSE value obtained from the EM algorithm. This also corresponds to a reduction of 56.6% compared to the standard, single-component two-part model (2.05 vs. 0.89). Figure (ref) also shows that the out-of-sample prediction errors obtained from the C-EM algorithm are generally lower and less variable than those obtained from the EM algorithm, especially when $G=2$. For parsimony, the subsequent subsections will focus on the set of parameter estimates that yielded the lowest RMSE among all values.

Groups' Analysis

Table (ref) shows the estimated mean and variance values of each covariate for each group of the optimal set of estimates according to the CV procedure. Note that the distributions of the covariates in each group are likely to overlap with each other, as denoted by the large within-group variance values. Although each group contains observations that are homogeneous with respect to their unobserved characteristics, this means that all groups might contain observations that appear very different from each other based on their observed characteristics. Because of this observed heterogeneity, it is hard to broadly characterize the observations contained within each group. This is why the expression “mostly feature” is used below for the description of each group.

The first group contains observations that mostly feature a moderately high number of comorbidities, as indicated by columns (4) and (7), and low continuity of care. The second group contains observations that mostly feature a larger number of comorbidities compared to the first group, but also a higher contemporaneous level of continuity of care, as indicated by column (5). The third group contains observations that mostly feature very low continuity of care and a moderately high number of comorbidities. Note that this group is the largest in size and also has the largest proportion of males. The fourth group contains observations that mostly feature a very low number of comorbidities and a relatively high continuity of care. This group also features the smallest proportion of males. Finally, the fifth group contains exclusively observations that feature perfect continuity of care. It is also the smallest group in size, and observations contained in this group mostly feature a relatively low number of comorbidities. Given that one covariate is constant within this group, this creates a perfect collinearity issue that can be dealt with by removing the group's intercept. Note that there exist large differences between the time-varying and the time-averaged COCI in most groups. For instance, the first and third groups contain observations that mostly feature low continuity of care at the current period, but high continuity of care on average over time. This is the opposite for the fourth and fifth groups. Therefore, it is not unreasonable to think that the level of continuity of care, both contemporaneously and over time, is an important factor in determining group memberships in the sample.

Figure (ref) shows the estimated time-fixed effects for each group and each part of the model that are associated with the same set of estimates. Confidence intervals on the time-fixed effects (and all other coefficients) are computed using Eq. ((ref)). The first period, ranging from 6 months to 3 months before the initial ED visit, is set as the reference value for the other time-fixed effects. The shaded areas correspond to the 95% cluster-robust confidence intervals. Given that the continuous part of the model is in the log scale, each mean value in the graphs on the right of Figure (ref) corresponds to the average, relative increase in individual HCE (in %) for each period. For instance, the observations in Group 1 feature an average increase of 70% in their individual HCE at $t=5$ compared to $t=1$, which is explained by the low levels of continuity of care and the poor health of the observations contained within this group. The marginal effects of the time-fixed effects in the binary parts can be obtained using the estimated coefficient values shown in Appendix (ref) and the covariates' average values presented in Table (ref). Nonetheless, the sign of the mean values of the time-fixed effects on the LHS of Figure (ref) still indicates the direction of each effect on the probability of consuming medical resources. The very large time-fixed effects depicted in Group 4 at $t=3$ and $t=4$ indicate that every or almost every observation in this group consumed a strictly positive amount of medical resources during those two periods. Note that the initial ED visit after the minor injury did not necessarily lead to a visit to the physician, thus explaining why some time-fixed effects at $t=2$ are significantly negative in both parts of the model.

The analogous set of time-fixed effects that are associated with the “best” set of parameters generated by the EM algorithm is shown in Appendix (ref). Trying to make formal connections between the estimates obtained from each algorithm is, according to me, hazardous given that the two methods do not treat the observed information similarly, and can produce results that are substantially different from each other. Such differences can also be appreciated by examining all other parameters estimated by each algorithm.

Transition Between Groups

Using the results from Table (ref) and Figure (ref), it is possible to broadly characterize each group. This is shown below in Figure (ref), along with the proportions of the total number of observations contained in each group at each period. This “naming” exercise is highly subjective and does not rely on a comprehensive analysis of all observed characteristics and results. The proposed “names” for the groups rely exclusively on frailty and continuity of care, which are known to be good predictors of individual HCE sirven_cost_2017, chu_continuity_2012, hussey_continuity_2014.\footnote{More precisely, it was assumed that a group is composed of frail individual-periods if at least one of the time-fixed effects of the continuous part is significantly positive after the initial injury, and if the group mean value for the ERA or the Charlson index is larger than its global average. It was assumed to be composed of robust individual-periods if no time-fixed effect of the continuous part is significantly different from zero after the initial injury. Overall continuity of care was based exclusively on the average value of the time-varying COCI covariate, as shown in column (6) of Table (ref), where the attribution of the “type” of continuity of care (i.e., low, moderate, and perfect) is obvious. Note that the overlap between the distributions of the COCI covariate for each type of continuity of care is very small compared to the other covariates, which advocates for the use of this covariate to broadly define the groups.} Note that the group numbers have been previously arranged in descending order so that the first group corresponds to the “high-cost” group, whereas the fifth group corresponds to the “low-cost” group. Figure (ref) shows that the initial ED visit substantially reduced the number of units in the group with perfect continuity of care at $t=2$. This is not surprising since the initial ED visit was associated with an unexpected minor trauma, therefore increasing the chance for the patient's regular healthcare provider to be unavailable at this precise moment. On the other hand, the initial ED visit seemed to have substantially increased the number of patients in Group 3. This can be explained by the fact that (robust) individuals who rarely consume medical resources can easily switch from perfect continuity to low continuity of care if they see a different healthcare provider each time they consume medical resources. Figure (ref) also shows that Group 4 experienced the largest reduction in size with a decrease of 3.4 percentage points between period 1 and period 7, while Group 2 is the group whose size increased the most with 2.3 additional percentage points between period 1 and period 7. This can be easily explained by a natural transition from robustness to frailty as time passes after the minor injury.

Figure (ref) depicts the estimated transition matrix of the group memberships over time. This transition matrix shows that robustness and frailty are persistent characteristics since transitions from a robust group to another robust group account for 74% of all transitions out of a robust group on average, whereas transitions from a frail group to another frail group account for 65% of all transitions out of a frail group on average. Therefore, robustness appears more persistent than frailty, which is in line with the fact that frailty is reversible provided appropriate care and healthy lifestyle habits kolle_reversing_2023. Such a transition matrix can also be used to predict the group membership at period $t+1$ given the group membership at period $t$. Using the Bayes' rule of categorical assignment and the probabilities shown in Figure (ref), group membership at period $t+1$ is correctly predicted 3,226 times out of 7,980 total memberships, which corresponds to a “success” rate of 40.4%. This success rate could be easily improved upon using more sophisticated approaches, such as dynamic multinomial logit models, namely by including additional observed factors that might affect future group memberships.

Conclusion

This paper showed that maximizing the likelihood function of a mixture density, as defined generally by Eq. ((ref)), leads to estimates that can be severely biased in finite samples under weak regularity conditions. It is worth emphasizing that the resulting biases do not depend on the chosen algorithm to maximize the mixture likelihood if the global maximum of the mixture likelihood is indeed reached. The motivation behind the widespread use of the EM algorithm is that it maximizes the mixture likelihood function just as any other numerical optimization method, but is easier to implement in practice. Still, it is true that for the same initial parameter values, different algorithms might converge to different local maxima. This, however, does not invalidate any of the claims made in this paper. Rather than relying on the mixture likelihood function, this paper showed that maximizing the classification-mixture likelihood function combined with a consistent classifier leads to (much) less biased and consistent estimates of all parameters in the mixture. As shown in Section (ref), consistency of the classifier is essential to obtain consistent estimates of the component parameters and the mixing weights. Although uniform consistency of the chosen classifier is a desirable property, Theorem (ref) showed that it is not necessary to obtain consistent estimates of all parameters in the mixture. Nonetheless, Corollary (ref) showed that the Mahalanobis distance classifier is uniformly consistent if the number of covariates goes to infinity at a faster rate than the sample size for a fixed number of groups, which is however problematic for estimation in large samples. Consistency of the Mahalanobis distance classifier (and, by extension, of the joint density classifier under normality) heavily hinges on the validity of Assumptions (ref)$(iv)-(v)$, which assumes that all the elements in the outcome vector $y_i$ and/or the covariates vector $x_i$ share the same group structure. To what extent the violation of those assumptions reduces the performance, and the usefulness of the proposed estimation strategy is an empirical question : stronger violations of those assumptions are likely to lead to worse performance of the C-EM algorithm and vice versa. A strong violation would be a case where the group structure of the conditional outcome density is completely different from that of the covariates' density. Because it is not difficult to rule out such a possibility before analyzing the data, it is a good practice to use both the EM and the C-EM algorithms to compare the predictive performance of each approach, as done in Section (ref) of this paper. Contrary to the recent papers on latent group panel structures, the estimation strategy proposed in this paper has the benefit of leaving group membership completely unrestricted. This has important implications from a policy perspective given that transitions between groups over time allow for the prediction of future group memberships, which can then be used to improve decision-making in healthcare. Achieving such an objective implies that group membership has to be consistently estimated in the first step. The dynamic behavior of the group membership can then be modeled in the second step using a consistent estimator for categorical outcomes. Such a two-step procedure is advocated by bonhomme_distributional_2019 in the context of matched employer-employee panel data analysis, while being more in line with the literature on finite mixtures, where no restriction is imposed on group memberships to estimate the mixture parameters. A general recommendation concerning the estimation of any kind of finite mixture models (e.g., Gaussian mixture models, mixture of experts, latent group panel structures, etc.) with unrestricted group membership can thus be formulated : instead of using the EM algorithm, one should always use the C-EM algorithm combined with a consistent classifier with as many covariates as possible. If the large number of covariates makes the computational burden too heavy, then one could remove from the conditional mixture density $f(y_i|x_i;\xi)$ (if it is the density of interest) the covariates that are deemed to be irrelevant for predicting the outcome $y_i$. If this is not sufficient, then a simple and effective procedure for variable selection at every classification step would be to remove covariates that feature the lowest in-sample variances. Indeed, covariates that feature a low variance do not contain a substantial amount of information in order to guide the classification process. Therefore, removing these covariates from all classification steps should not significantly increase the misclassification rate while easing computational burden. Results from the simulation exercises and a real-world application showed that the estimation procedure leads to less biased and more stable estimates than those obtained by standard MLE procedures such as the EM algorithm. Results from the empirical analysis also show that the proposed estimation strategy identified five heterogeneous health groups that account for a large part of the unobserved heterogeneity in the sample. The use of the proposed estimation strategy reduced the out-of-sample prediction error by more than 55% compared to the single-component model and by 17.6% compared to the best results obtained from standard, widely used MLE procedures. Compared to other algorithms, the proposed estimation strategy contrasts with HMMs, where the dynamic behavior of the latent variable is estimated simultaneously with the parameters of each component density. Although this method has been proven to yield consistent estimates douc_asymptotics_2001, the two-step estimation procedure proposed in this paper seems to be less biased and more efficient than HMMs since the latter is often estimated using special cases of the EM algorithm such as the Baum-Welch or the Viterbi algorithm bishop_pattern_2006. Confirming this intuition goes beyond the scope of this paper, but it is an interesting avenue for future research. Introduction of endogeneity in the covariates of the conditional mixture density, namely through feedback effects, is also another avenue for future research chamberlain_feedback_2022. \setstretch{1.15}