Assessing Heterogeneity of Two-Part Model via Bayesian Model-Based Clustering with Its Application to Cocaine Use Data

Abstract

The purpose of this chapter is to provide an introduction to the model-based clustering within the Bayesian framework and apply it to asses the heterogeneity of fractional data via finite mixture two-part regression model. The problems related to the number of clusters and the configuration of observations are addressed via Markov Chains Monte Carlo (MCMC) sampling method. Gibbs sampler is implemented to draw observations from the related full conditionals. As a concrete example, the cocaine use data are analyzed to illustrate the merits of the proposed methodology.

1. Introduction

A recurring theme in the statistical analysis is to separate the unstructured data into groups to detect the similarity or discrepancy within or between groups. This is especially true in the fields, e.g., discriminant analysis [1, 2, 3], pattern recognition [4, 5], gene expression [6, 7, 8], machine learning [9], and artificial intelligence [10]. In the literature, the clustering problem is often formulated within the cluster analysisframework, which is generally categorized into two classes: the non-probabilistic framework and the probabilistic framework. The non-probabilistic clustering method, including the $K$-means method [9, 11, 12] and the hierarchical/agglomerative clustering algorithms [13, 14, 15], is based on the distancebetween any two observations or groups. It clusters data by merging or removing observations according to the “closeness” specified by the distance. This method is more general since it does not impose any distributional assumptions on data, hence having greater flexibility in the real applications. Instead, the non-probabilistic clustering algorithm, also termed the model-based clustering, groups data by positing a probability model on data and then clustering data via configuration function related to the model. Compared with the non-probabilistic framework, the model-based methods enable us to assess the statistical properties of the solutions, e.g., how many clusters are there, how well the configuration function works, and how robust the method is against the model deviation and so on. There is rich literature on this issue. Among them, finite mixture model (FMM, [16, 17, 18]) perhaps is the most popular choice and has often been proposed and studied in the context of clustering (see a short review in Fraley and Raftery [2]). FMM assumes that each cluster is identified with a probability distribution indexed by the cluster-specific parameter(s), and each observation is related to clusters via configuration or membership function. The statistical task is the inference about the number of clusters, the estimation of the unknown parameters, and the allocation of observations.

In this chapter, we pursue a Bayesian model-based method to address the heterogeneity of fraction data. Fractional data are very common in the social and economical surveys. A distinguished feature of fractional responses is that its measurements are responded on a scale in the unity interval [0,1] but suffer from excessive zeros and unities on the boundaries. In understanding such type of data, the commonly used method is to separate the whole data into three parts: two corresponding to the zeros and unities respectively, and one corresponding to the continuously positive values. Two separative logistic models are suggested to model two discrete value parts respectively while single normal linear regression model is formulated for the continuous value part. This method, though more appealing, ignores the instinct association across different parts and readily leads to inconsistence of the occurrence probabilities on each part. Instead, we propose a three-category multinomial model for the occurrence variable, in which the usual separated models can be considered as the marginal models of our proposal. Such modeling always ensures the probabilities on each part to be proper, thus avoiding parameter constraints, see for example, [19]. To assess the heterogeneity underlying data, we formulate the problem into a finite mixture analysis of which each component is specified by two-part regression model. In view of the model complexity, we implement Markov Chains Monte Carlo sampling method to implement posterior analysis. Block Gibbs sampler is implemented to draw observations from the target distributions. The posterior inference including parameters estimates, model selection, and the configuration determination of observations are obtained based on the simulated observations.

The chapter is organized as follows. Section 2 introduces a general model-based clustering method to address the heterogeneity of regression model within the Bayesian framework. In Section 3, we apply the proposed method to the fractional data. Section 4 presents a cocaine use study. And Section 5 concludes the chapter.

2. Method description

2.1 General framework

Suppose that for $i=1,2,\cdots ,n$, $y_i$is an observed response, each associated with an $m$dimensional fixed covariates ${\mathbf{x}}_i=\left(x_{i1},\cdots ,x_{\mathit{im}}\right)$. In the context of regression analysis, the interest mainly focuses on exploring the pattern of the influence of ${\mathbf{x}}_i$on $y_i$and predicting the mean of a future response $y$in terms of a new $\mathbf{x}$. This is usually achieved by formulating $\left\{{\mathbf{x}}_i,y_i\right\}$as $\mathbb{E}\left(y_i|{\mathbf{x}}_i\right)=m\left({\mathbf{x}}_i\right)$for some mean function $m\left(\cdot \right)$. In the parametric fitting framework, the function $m\left(\mathbf{x}\right)$is assumed to be related to $\mathbf{x}$via linking function as the form of

$$m\left(\mathbf{x}\right)=h\left({\mathbf{x}}^T\mathit{\beta }\right)$$E1

which induces the so-called generalized linear model [20] for $\left\{{\mathbf{x}}_i,y_i\right\}$, where $\mathit{\beta }$is the regression coefficients used to quantify the uncertainty about $m$, and $h\left(\cdot \right)$is the known linking function used to link the mean and the predictors.

More often, the single relationship such as Eq. (1) may not be sufficient when the patterns among the subjects take on the heterogeneity such as clustering. The heterogeneous data occur when the observations are generated from the different populations of which the number of populations and the membership of each observation to the population are unknown. The main objective is to separate data into different clusters to detect the possible similarity within clusters or the discrepancy between clusters. This is generally accomplished by defining a cluster’s membership/configuration function $\mathcal{K}:\left\{\left({\mathbf{x}}_1,y_1\right),\cdots ,\left({\mathbf{x}}_n,y_n\right)\right\}\mapsto \left\{1,\cdots ,K\right\}$such that $K_i=\mathcal{K}\left(\left({\mathbf{x}}_i,y_i\right)\right)=k$if $\left({\mathbf{x}}_i,y_i\right)$belongs to the cluster $k$, where $K$is assumed to be less than $n$. The discrepancy between any two clusters is characterized by the cluster-specific parameters such as intercepters, regression coefficients, and/or disperse parameters.

The model-based clustering assumes that given the clusters membership $K_i$, $\left({\mathbf{x}}_i,y_i\right)$within the cluster $k$has the following sampling density

$$\left(y_i|K_i=k,{\mathbf{x}}_i\right)\stackrel{\mathit{ind}.}{\sim }{f}_k\left(y_i|{\mathbf{x}}_i^T{\mathit{\beta }}_k,{\tau }_k\right)$$E2

while $K_i$is specified by

$$\mathbb{P}\left(K_i=k\right)={\pi }_k$$E3

where $f_k$, maybe independent of $k$, is the probability density function, ${\mathit{\beta }}_k$and ${\mathit{\tau }}_k$are the cluster-specific regression coefficients and the disperse parameters, respectively, and ${\pi }_k$is the mixing proportion identifying the proportion of the component $k$over the entire population. It is assumed that ${\pi }_k\ge 0$and ${\sum }_{k=1}^K{\pi }_k=1.0$.

Two important issues arise when formulating data clustering problem as Eqs. (2) and (3). One is related to the number of clusters, and the other is pertained to the determination of configurations. Within the Bayesian framework, several methods have been proposed for the first issue. One can, for example, follow [21] and treat $K$to be random and assign a prior to it. The reversible jump MCMC method (RJMCMC, [21, 22]) can be implemented to conduct the joint analysis of $K$with other random quantities. Another method is along the lines with the hypothesis test procedure and routinely to estimate $K$via model comparison/selection procedure. This perhaps is the most popular choice in the model-based clustering context, in which various measures such as the Akaike information criterion (AIC) [23], the corrected AIC (AICc) [24, 25], the Bayesian information criterion (BIC) [26], the integrated completed likelihood (ICL) [27], and Bayes factor (BF, [28, 29]) can be adopted to select a suitable model. It is worth pointing out that the deviance information criterion (DIC) [30] may not be appropriate for the mixture model comparison. The well-known software WinBUGS® [31] for Bayesian analysis does not provide DIC results for mixture analysis. In addition, many authors suggested modeling heterogeneous data into the mixture of Dirichlet process (MDP, [32, 33]). However, as discussed in Ishwaran and James [34], DP fitting often overestimates the number of clusters and readily leads to model over fitting.

For the second issue, the complexity of problem depends on the methods adopted in the analysis. In the frequency framework, for example, the configuration of observation $i$is often achieved by maximizing $\mathbb{P}\left(K_i=k|\mathbf{Y},\widehat{\mathit{\pi }},\widehat{\mathbf{\Xi }}\right)$over $k=1,\cdots ,K$, where $\widehat{\mathit{\pi }}$and $\widehat{\mathbf{\Xi }}$are the maximum likelihood estimates (MLE) obtained via, e.g., the expectation-maximization algorithm (EM, [35]). In the next section, we will present a Bayesian procedure for determining $\mathcal{K}$. Compared with the frequency approach, the nice feature of the Bayesian approach is its flexibility to utilize prior information for achieving better results. Also, the sampling-based Bayesian methods depend less on the asymptotic theory and hence have the potential to produce reliable results even with small sample size.

Let $\mathbf{Y}$be the set of all observed responses and $\mathbf{X}$be the set of fixed covariates; Write $\mathbf{\Xi }$as the collection of ${\mathit{\beta }}_k$and ${\mathit{\tau }}_k$. Integrating over $K_i$produces a $K$-component mixture model for $y_i$, which is given by

$$p\left(y_i|\mathit{\pi },\mathbf{\Xi },{\mathbf{x}}_i\right)=\sum _{k=1}^K{\pi }_k{f}_k\left(y_i|{\mathbf{x}}_i^T{\mathit{\beta }}_k,{\mathit{\tau }}_k\right).$$E4

The log-likelihood of the observed data conditional on $K$is given by

$$\mathcal{L}\left(\pi ,\mathbf{\Xi }|\mathbf{Y},\mathbf{X}\right)=\sum _{k=1}^{n}log\left(\sum _{k=1}^K{\pi }_k{f}_k\left(y_i|{\mathbf{x}}_i^T{\mathit{\beta }}_k,{\mathit{\tau }}_k\right)\right).$$E5

As an illustration, Figure 1 presents a three-component normal linear mixture regression model with one covariate. It can be seen clearly that the density function illustrates strong heterogeneity. The regression line is obviously different from those of components, which indicates that single model is unappreciate in fitting such data. In what follows, we suppress $\mathbf{X}$for notational simplicity.

Figure 1.

Plot of the three-component normal mixture model$0.3N\left(-4-2x,1\right)+0.5N\left(0.5+0.5x,1\right)+0.2N\left(4.5+3x,1\right)$. Left panel: Plot of the density functions of the mixture as well as their three weighted components ; right panel: plots of regression lines. Mixture model: solid line “$-$” component one: dotted lines “$\cdots$” component two: dashed lines “$--$” and component three: dotted-dashed lines “$-\cdot$”

2.2 Bayesian model-based clustering via MCMC

Bayesian analysis for analyzing Eqs. (2) and (3) especially $\mathcal{K}$requires the specification of a prior distribution $p\left(\mathit{\pi },\mathbf{\Xi }\right)$for the parameters of the mixture model. By model convention, it is naturally to assume that $\mathit{\pi }$and $\mathbf{\Xi }$are independent, and the components among $\mathbf{\Xi }$are also independent. In particular,

$${\mathit{\beta }}_k\stackrel{\mathit{iid}.}{\sim }{N}_m\left({\mathit{\beta }}_0,{\mathbf{\Sigma }}_0\right),{\mathit{\tau }}_k^{-1}\stackrel{\mathit{iid}.}{\sim }W\left({\rho }_0,{\mathit{R}}_0\right)$$E6

in which $W\left({\rho }_0,{\mathit{R}}_0^{-1}\right)$is the Wishart distribution with the degrees of freedom ${\rho }_0$and the scale matrix ${\mathit{R}}_0$, and reduces to the scaled Chi-square distribution when ${\mathit{\tau }}_k$is a univariate; ${\mathit{\beta }}_0$, ${\mathbf{\Sigma }}_0$, ${\rho }_0$and ${\mathit{R}}_0$are the hyper-parameters, which are treated to fixed and known. In the real applications, if no extra information can be available, the values of these hyper-parameters are often taken to ensure ${\mathit{\beta }}_k$and ${\mathit{\tau }}_k$to be dispersed enough. For example, one can set ${\mathbf{\Sigma }}_0={\lambda }_0\mathbf{I}$with large ${\lambda }_0$(Throughout, we use $\mathbf{I}$to signify an identify matrix). In this case, the values of ${\mathit{\beta }}_0$are not really important and can be set to any values, e.g., zeros. Note that for the mixture models, Diebolt and Robert [36] (see also, for example, [37]) showed that using fully non-informative prior distributions may lead to improper posterior distributions and hence is strictly prohibitive.

We assign a symmetric Dirichlet distribution to $\mathit{\pi }$as follows

$$\mathit{\pi }\mid \alpha \sim D_K\left(\alpha ,\cdots ,\alpha \right)$$E7

in which $\alpha \left(>0\right)$is the hyper-parameter, which is treated to fixed and unknown. In the applications, we can take sensitive analysis by setting smaller and larger values for $\alpha$. See section 4 for more details.

Let $\mathbf{K}=\left\{K_1,\cdots ,K_n\right\}$be the collection of all configurations. A Bayesian procedure for model-based clustering mainly focuses on exploring the behavior of the posterior of $\mathbf{K}$given data, which is given by

$$p\left(\mathbf{K}|\mathbf{Y}\right)\propto p\left(\mathbf{Y}|\mathbf{K}\right)p\left(\mathbf{K}\right)$$E8

where $p\left(\mathbf{Y}|\mathbf{K}\right)$is the marginal distribution of $p\left(\mathbf{Y},\pi ,\mathbf{\Xi }|\mathbf{K}\right)$with $\mathit{\pi }$and $\mathbf{\Xi }$being integrated out. Generally, no closed form can be available for this target distribution. Markov Chain Monte Carlo [38, 39] sampling method can be used to conduct posterior analysis. In particular, one can follow the routine in Tanner and Wong [40] and treat the latent quantities $\left\{\mathit{\pi },\mathbf{K},\mathbf{\Xi }\right\}$as the missing data and augment them with the observed data. Posterior analysis is carried out based on the joint distribution $p\left(\mathit{\pi },\mathbf{K},\mathbf{\Xi }|\mathbf{Y}\right)$. In this case, block Gibbs sampler [41, 42] can be implemented to draw observations from such target distribution. The Gibbs sampler is iteratively implemented by drawing: (i) $\mathrm{\Xi }$from $p\left(\mathbf{\Xi }|\mathit{\pi },\mathbf{K},\mathbf{Y}\right)$; (ii) $\mathit{\pi }$from $p\left(\mathit{\pi }|\mathbf{K},\mathbf{\Xi },\mathbf{Y}\right)$and $\mathbf{K}$from $p\left(\mathbf{K}|\mathit{\pi },\mathbf{\Xi },\mathbf{Y}\right)$till convergence. The convergence can be monitored by the “estimated potential scale reduction” (EPSR) values [43] or by plotting the traces of estimates against iterations under different starting values. Note that except for (i), all full conditionals involved in the Gibbs sampler are standard. However, drawing $\mathbf{\Xi }$in (i) depends on the specific form of the density function $f_k$and sometimes requires implementing Metropolis-Hastings algorithm (MH, [44, 45]) or rejection sampling [46].

2.3 Label switching

Formulating the model-based clustering problem into mixture model Eq. (2) faces the model identification. A statistical model is said to be identified if the observed likelihood is uniquely determined by unknown parameters. A less identified model may be problematic and will distort the estimates of unknown parameters. It is easily showed that the observed likelihood of data is only determined up to the permutation of the component labels. As a matter of fact, suppose that there are the pair $\left\{{\mathit{\pi }}^{\left(1\right)},{\mathbf{\Xi }}^{\left(1\right)}\right\}$and $\left\{{\mathit{\pi }}^{\left(2\right)},{\mathbf{\Xi }}^{\left(2\right)}\right\}$such that

$$p\left(y|{\mathit{\pi }}^{\left(1\right)},{\mathbf{\Xi }}^{\left(1\right)}\right)=p\left(y|{\mathit{\pi }}^{\left(2\right)},{\mathbf{\Xi }}^{\left(2\right)}\right)$$E9

then there exists a permutation $\nu :\left\{1,2,\cdots ,K\right\}\mapsto \left\{1,2,\cdots ,K\right\}$such that ${\pi }_k^{\left(1\right)}={\pi }_{\nu \left(k\right)}^{\left(2\right)}$, ${\mathit{\beta }}_k^{\left(1\right)}={\mathit{\beta }}_{\nu \left(k\right)}^{\left(2\right)}$and ${\mathit{\tau }}_k^{\left(1\right)}={\mathit{\tau }}_{\nu \left(k\right)}^{\left(2\right)}$. In this setting, we can not distinguish $\mathcal{K}$and $\nu \circ \mathcal{K}$in terms of data (“$\circ$” denotes the operator of function composition). With this in mind, any exchangeable priors on $\mathit{\pi }$and $\mathbf{\Xi }$like Eqs. (6) and (7) produces symmetric and multi-modal posterior distributions with up to $K!$copies of each “genuine” mode, which induces the so-called label switching problem on Bayesian estimate. Traditional approaches to eliminating such exchangeability is to impose identifiability constraints on the parameter space. However, as pointed out by Frühwirth-Schnatter [18], an unappropriate identifiability constraint may not be able to eliminate label switching. Many efforts have been devoted to coping with this issue, see Chapter 11 in Lee [47] for a review. Among them, the relabeling algorithm [48] is more appealing due to its simplicity and flexibility. The relabeling sampling procedure takes a decision-theoretical approach and requires specifying an appropriate loss function to measure the loss in terms of the classification probability. The model identification problem is addressed via postprocessing the MCMC output to minimize the posterior expected loss. Specifically, let $\mathit{\theta }$be the collection of $\mathbf{\Xi }$and $\mathit{\pi }$, and write $\mathit{Q}=\{q_{\mathit{ik}}\left(\mathit{\theta }\right)$as the matrix of allocation probabilities of order $n\times K$with $q_{\mathit{ik}}\left(\theta \right)=\mathbb{P}\left(K_i=k|\mathbf{Y},\mathit{\theta }\right)$. In the context of clustering, the loss function can be defined on the cluster label $\mathcal{K}$as follows

$${\mathcal{L}}_0\left(\mathcal{K};\mathit{\theta }\right)=-\sum _{i=1}^{n}log{q}_{{\mathit{iK}}_i}\left(\mathit{\theta }\right).$$E10

Given that ${\mathit{\theta }}^{\left(1\right)},\cdots ,{\mathit{\theta }}^{\left(M\right)}$are the sampled parameters and let ${\nu }_1,\cdots ,{\nu }_M$be the permutation applied to them. The relabeling algorithm proceeds by selecting initial values for the ${\nu }_m$s, which are generally taken to be the identity permutations, then iterating the following steps until a fixed point is reached.

1. Choose $\widehat{\mathcal{K}}$to minimize ${\sum }_{m=1}^M{\mathcal{L}}_0(\mathcal{K},{\nu }_m\left({\mathit{\theta }}^{\left(m\right)}\right)$;

2. For $m=1,2,\cdots ,M$, choose ${\nu }_m$to minimize ${\mathcal{L}}_0(\widehat{\mathcal{K}},v_m\left({\mathit{\theta }}^{\left(m\right)}\right)$.

2.4 Posterior inference

Once the label switching is taken care of, the MCMC samples can be used to draw posterior inference. For example, the joint Bayesian estimate of $\mathit{\theta }$can be obtained easily via the corresponding sample means of the generated observations via ergodic average as follows:

$${\widehat{\mathit{\beta }}}_k=M^{-1}\sum _{m=1}^M{\mathit{\beta }}_k^{\left(m\right)},{\widehat{\mathit{\tau }}}_k=M^{-1}\sum _{m=1}^M{\mathit{\tau }}_k^{\left(m\right)},\text{and}{\widehat{\pi }}_k=M^{-1}\sum _{m=1}^M{\pi }_k^{\left(m\right)}$$E11

The consistent estimates of the covariance matrix of estimates can be obtained via sample covariance matrix.

Given the observations $\left\{{\mathbf{K}}^{\left(m\right)}:m=1,2,\cdots ,M\right\}$drawn from the posterior $p\left(\mathbf{K}|\mathbf{Y}\right)$via MCMC sampling, serval methods can be available for arriving at a point estimate of the clustering using draws from the posterior clustering distribution. The simplest method, known as the maximum a posteriori (MAP) clustering, is to select the observed clustering that maximizes the density of the posterior clustering distribution, i.e.,

$$\widehat{\mathcal{K}}:{\widehat{K}}_i={\text{argmax}}_{k=1,\cdots ,K}\mathbb{P}\left(K_i=k|\mathbf{Y}\right)$$E12

in which $\mathbb{P}\left(K_i=k|\mathbf{Y}\right)$can be approximated by

$$\mathbb{P}\left(K_i=k|\mathbf{Y}\right)\approx M^{-1}\sum _{m=1}^{M}I\left\{K_i^{\left(m\right)}=k\right\}.$$E13

A more appreciate alternative to MAP is based on the pairwise probability matrix, an $n\times n$association matrix $\delta \left(\mathcal{K}\right)$with the $\left(i,j\right)$th element formed by the indicator of whether the subject $i$is clustered with subject $j$. Element-wise averaging of these association matrices yields the pairwise probability matrix of clustering, denoted $\widehat{\psi }$. Medvedovic and Sivaganesan [49] and Medvedovic et al. [50] suggested a clustering estimate of $\mathcal{K}$by using the pairwise probability matrix $\widehat{\psi }$as a distance matrix in hierarchical/agglomerative clustering. However, as augured by Dahl [51], such routine seems counterintuitive to apply an ad hoc clustering method on top of a model which itself produces clusterings. In the context of Dirichlet process mixture-based clustering, Dahl [51] proposed a least-squares model-based clustering method by using draws from a posterior clustering distribution. Specifically, the least-squares clustering ${\mathcal{K}}_{\mathit{LS}}$is the observed clustering ${\mathcal{K}}_{\mathit{LS}}$, which minimizes the sum of squared deviations of its association matrix $\left(\mathcal{K}\right)$from the pairwise probability matrix:

$${\widehat{\mathcal{K}}}_{\mathit{LS}}={\text{argmin}}_{\mathcal{K}\in \left\{{\mathbf{K}}^{\left(1\right)},\cdots ,{\mathbf{K}}^{\left(m\right)}\right\}}\sum _{i=1}^n\sum _{j=1}^n{\left(\delta \left(i,j\right)\left(\mathcal{K}\right)-\widehat{\psi }\left(i,j\right)\right)}^2.$$E14

Dahl [51] showed that the least-squares clustering has the advantage over those in Medvedovic and Sivaganesan [49] since it utilizes the information from all the clusterings and is intuitively appealing for the “average” clustering instead of forming a clustering via an external, ad hoc clustering algorithm.

3. Assessing heterogeneity of two-part model

In this section, we first proposed a two-part regression model for the fractional data especially for the U shaped fractional data and then extend the method discussed above to the current situation to address the possible heterogeneity of the population underlying data.

3.1 Two-part model for U shaped fractional data

Suppose that for subject/individual $i\left(=1,\cdots ,n\right)$, $y_i$is an univariate fractional response taking values in $\left[0,1\right]$; ${\mathbf{x}}_i$is an $m\times 1$fixed covariate vector denoting various explanatory factors under consideration. Usually, $y_i$suffers from excess zeros and ones on the boundaries, and the whole data set takes on the U shape. In modeling such data, we introduce a three-category indicator variable $d_i$and a continuous intensity variable $z_i$such that

$$d_i=\left\{\begin{array}{lll}1& \text{if}& y_i=0\\ 2& \text{if}& y_i=1\\ 3& \text{if}& 0<y_i<1\end{array}\right.\text{and}{z}_i=\left\{\begin{array}{ccc}h\left(y_i\right)& \text{if}& 0<y_i<1\\ \text{irrelevant}& \text{if}& y_i=0,1\end{array}\right.$$E15

where $h\left(\cdot \right)$is any monotone increasing function such that $z_i\in \left(-\infty ,+\infty \right)$. That is, we break the data set into three parts: two parts corresponding to zeros and ones respectively and one part corresponding to the continuous values between 0 and 1. We formulate a two-part model for $y_i$by first specifying a baseline-category logits model [52] for $d_i$and then a conditional continuous model for $z_i$. The baseline-category logits model is assumed that conditional upon ${\mathbf{x}}_i$, $d_i$s are independent satisfying the following logits models simultaneously: for $j=1,2$,

$$log\frac{\mathbb{P}\left(d_i=j|{\mathbf{x}}_i\right)}{\mathbb{P}\left(d_i=3|{\mathbf{x}}_i\right)}={\mathbf{x}}_i^T{\mathit{\alpha }}_j$$E16

where ${\mathit{\alpha }}_j$is an $m\times 1$regression coefficients vector. We use category $d_i=3$as the reference for the ease of parameters interpretation. For example, the magnitude of ${\alpha }_{j\mathrm{\ell }}$in ${\mathit{\alpha }}_j$indicates that the increase of one unit in $x_{i\mathrm{\ell }}$will increase $e^{{\alpha }_{j\mathrm{\ell }}}$times chance of $d_i=j$over that of $d_i=3$.

The conditional continuous model for $z_i$is given by

$$p\left(z_i|d_i=3,{\mathbf{x}}_i\right)=p^z\left(z_i|{\mathbf{x}}_i^T\mathit{\gamma },\tau \right)$$E17

or equivalently

$$p\left(y_i|0<y_i<1,{\mathbf{x}}_i\right)=p^z\left(h\left(y_i\right)|{\mathbf{x}}_i^T\mathit{\gamma },\psi \right)\mid \dot{h}\left(y_i\right)\mid$$E18

where $\dot{h}\left(s\right)=\mathit{dh}/\mathit{ds}$, $p^z\left(u|a,\tau \right)$is the normal density with mean $a$and variance $\tau >0$, and $\mathit{\gamma }$like that in Eq. (16), is the regression coefficient vector. Although the identical covariates are taken in Eqs. (16) and (17), this is not necessary in practice. Each equation can own their covariates. This can be achieved by imposing particular structure on the regression coefficients. For example, we can exclude $x_{i1}$from Eq. (17) by restricting ${\gamma }_1$in $\mathit{\gamma }$to be zero.

It follows from Eqs. (16) and (17) that marginal distribution of $y_i$is given by

$$p\left(y_i|{\mathbf{x}}_i,\mathit{\beta },\tau \right)=q_{i1}{\delta }_0+q_{i2}{\delta }_1+\left(1-q_{i1}-q_{i2}\right)p\left(y_i|0<y_i<1,{\mathbf{x}}_i,\mathit{\gamma },\tau \right)$$E19

where $q_{\mathit{ij}}=\mathbb{P}\left(d_i=j|{\mathbf{x}}_i,{\mathit{\alpha }}_j\right)\left(j=1,2\right)$is the response probability specified by Eq. (16) and $\mathit{\beta }$is the regression parameters constituted by ${\mathit{\alpha }}_1,{\mathit{\alpha }}_2$and $\mathit{\gamma }$.

3.2 Assessing heterogeneity of two-part model

To detect the possible heterogeneity among $y_i$, we extend the model Eq. (18) to the mixture case by assuming that conditional upon $K_i=k$, $d_i$and $z_i$satisfy Eqs. (16) and (17) with ${\mathit{\alpha }}_j$replaced by ${\mathit{\alpha }}_{\mathit{jk}}$and $\left(\mathit{\gamma },\tau \right)$by $\left({\mathit{\gamma }}_k,{\tau }_k\right)$respectively. This indicates that the mixture component $f_k$in Eq. (1) in Section 2 is given by Eq. (19) with $\mathit{\beta }={\mathit{\beta }}_k$and $\tau ={\tau }_k$.

For the Bayesian analysis, the general forms of full conditionals involved in the model-based clustering have been given in Section 2. We here only focus on the technical details of the conditional distribution of $\mathbf{\Xi }$in (i) in the Gibbs sampler.

We assume that the prior of ${\tau }_k$is the same as that in Eq. (6), while the priors of ${\mathit{\beta }}_k$are taken as $p\left({\mathit{\beta }}_k\right)=p\left({\mathit{\alpha }}_{k1}\right)p\left({\mathit{\alpha }}_{k2}\right)p\left({\mathit{\gamma }}_k\right)$, in which

$$p\left({\mathit{\alpha }}_{k\mathrm{\ell }}\right)\stackrel{D}{=}{N}_m\left({\mathit{\alpha }}_{\mathrm{\ell }0},{\mathbf{\Sigma }}_{\alpha \mathrm{\ell }0}\right)\left(\mathrm{\ell }=1,2\right),p\left({\mathit{\gamma }}_k\right)\stackrel{D}{=}{N}_m\left({\mathit{\gamma }}_0,{\mathbf{\Sigma }}_{\gamma 0}\right).$$E20

where ${\mathit{\alpha }}_{\mathrm{\ell }0}$, ${\mathit{\gamma }}_0$, ${\mathbf{\Sigma }}_{\alpha \mathrm{\ell }0}$and ${\mathbf{\Sigma }}_{\gamma 0}$are the hyper-parameters treated to be known.

Gibbs sampling $\mathbf{\Xi }$now becomes drawing ${\mathit{\alpha }}_k$, ${\mathit{\gamma }}_k$and ${\tau }_k$alternatively from the full conditional distributions $p\left({\mathit{\alpha }}_k|\mathbf{K},\mathbf{Y}\right)$, $p\left({\mathit{\gamma }}_k|{\tau }_k,\mathbf{K},\mathbf{Y}\right)$and $p\left({\tau }_k|{\mathit{\gamma }}_k,\mathbf{K},\mathbf{Y}\right)$respectively. By some algebras, it can be shown that

$$\begin{array}{l}p\left({\mathit{\alpha }}_k|\mathbf{K},\mathbf{Y}\right)\propto p\left({\mathit{\alpha }}_k\right)\prod _{K_i=k}p\left(d_i|{\mathbf{x}}_i,{\mathit{\alpha }}_k\right),\\ p\left({\mathit{\gamma }}_k|{\tau }_k,\mathbf{K},\mathbf{Y}\right)\propto p\left({\mathit{\gamma }}_k\right)\prod _{K_i=k}p\left(z_i|d_i=3,{\mathbf{x}}_i^T{\mathit{\gamma }}_k,{\tau }_k\right),\\ p\left({\tau }_k|{\mathit{\gamma }}_k,\mathbf{K},\mathbf{Y}\right)\propto p\left({\tau }_k\right)\prod _{K_i=k}p\left(z_i|d_i=3,{\mathbf{x}}_i^T{\mathit{\gamma }}_k,{\tau }_k\right)\end{array}$$E21

in which the full conditionals of ${\gamma }_k$and ${\tau }_k$are easily obtained and given by

$$p\left({\mathit{\gamma }}_k|{\tau }_k,\mathbf{K},\mathbf{Y}\right)\stackrel{D}{=}N\left({\widehat{\mathit{\gamma }}}_k,{\widehat{\mathbf{\Sigma }}}_{\mathit{\gamma k}}\right)$$E22

$$p\left({\tau }_k^{-1}|{\mathit{\gamma }}_k,\mathbf{K},\mathbf{Y}\right)\stackrel{D}{=}\mathit{\text{Gamma}}\left({\widehat{\alpha }}_k,{\widehat{\beta }}_k\right)$$E23

in which

$$\begin{array}{l}{\widehat{\mathbf{\Sigma }}}_{\mathit{\gamma k}}={\left(\sum _{K_i=k:d_i=2}{\mathbf{x}}_i{\mathbf{x}}_i^T/{\tau }_k+{\mathbf{\Sigma }}_{\gamma 0}^{-1}\right)}^{-1},\\ {\widehat{\mathit{\gamma }}}_k={\widehat{\mathbf{\Sigma }}}_k\left({\mathbf{\Sigma }}_{\gamma 0}^{-1}{\gamma }_0+\sum _{K_i=k,d_i=3}{\mathbf{x}}_i{z}_i/{\tau }_k\right),\\ {\widehat{\alpha }}_k={\alpha }_0+n_k/2,\\ {\widehat{\beta }}_k={\beta }_0+\sum _{K_i=k,d_i=3}{\left(z_i-{\mathbf{x}}_i^T{\mathit{\gamma }}_k\right)}^2/2\end{array}$$E24

and $n_k=\#\left\{K_i=k,d_i=3\right\}$.

However, drawing ${\mathit{\alpha }}_{k\mathrm{\ell }}$is more tedious since its distribution loses the standard form. We first note that

$$p\left({\mathit{\alpha }}_{k\mathrm{\ell }}|{\mathit{\alpha }}_{k,-\mathrm{\ell }},\mathbf{K},\mathbf{Y}\right)\propto p\left({\mathit{\alpha }}_{k\mathrm{\ell }}\right)\prod _{K_i=k}^n\frac{exp\left({\tilde{d}}_{i\mathrm{\ell }}\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)\right)}{1+exp\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)}$$E25

where ${\tilde{d}}_{i\mathrm{\ell }}=I\left\{d_i=\mathrm{\ell }\right\}$and $C_{\mathit{ik}\mathrm{\ell }}=log\left(1.0+exp\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k,-\mathrm{\ell }}\right)\right)$; ${\mathit{\alpha }}_{k,-\mathrm{\ell }}$denotes the set ${\mathit{\alpha }}_k$with ${\mathit{\alpha }}_{k\mathrm{\ell }}$removed. Following the similar routine in Polson, Scott, and Windle [53], we recast the logistic function Eq. (25) as follows

$$\begin{array}{l}\frac{exp\left({\tilde{d}}_{i\mathrm{\ell }}\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)\right)}{1+exp\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k1}-C_{\mathit{ik}\mathrm{\ell }}\right)}=2^{-1}exp\left\{{\kappa }_{i1}\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)\right\}\\ \times {\int }_0^{\infty }exp\left\{-\frac{1}{2}{\omega }_{i\mathrm{\ell }}{\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)}^2\right\}{p}_{\mathrm{PG}}\left({\omega }_{i\mathrm{\ell }}\right)\mathrm{d}{\omega }_{i\mathrm{\ell }}\end{array}$$E26

in which ${\kappa }_{i\mathrm{\ell }}={\tilde{d}}_{i\mathrm{\ell }}-1/2$and $p_{\mathrm{PG}}$is the well-known PG$\left(1,0\right)$density function [53]. If one introduces $n$independent Pólya-Gamma variables ${\omega }_{i\mathrm{\ell }}$into the current analysis, then,

$$p\left({\omega }_{i\mathrm{\ell }}|{\mathit{\alpha }}_{k\mathrm{\ell }},{\mathit{\alpha }}_{k,-\mathrm{\ell }},\mathbf{K},\mathbf{Y}\right)\stackrel{D}{=}\mathrm{PG}\left(1,\left({\mathbf{x}}_i^T{\mathit{\alpha }}_{k\mathrm{\ell }}-C_{\mathit{ik}\mathrm{\ell }}\right)\right)$$E27

$$p\left({\mathit{\alpha }}_{k\mathrm{\ell }}|{\mathit{\alpha }}_{k,-\mathrm{\ell }},\mathrm{\Omega },\mathbf{K},\mathbf{Y}\right)\stackrel{D}{=}N\left({\widehat{\mathit{\alpha }}}_{k\mathrm{\ell }},{\widehat{\mathbf{\Sigma }}}_{\mathit{\alpha k}\mathrm{\ell }}\right)$$E28

where

$${\widehat{\mathbf{\Sigma }}}_{\mathit{\alpha k}\mathrm{\ell }}={\left(\sum _{K_i=k}{\mathbf{x}}_i{\mathbf{x}}_i^T{\omega }_{i\mathrm{\ell }}+{\mathbf{\Sigma }}_{\alpha \mathrm{\ell }0}^{-1}\right)}^{-1},{\widehat{\mathit{\alpha }}}_{k\mathrm{\ell }}={\widehat{\mathbf{\Sigma }}}_{\mathit{\alpha k}\mathrm{\ell }}\left({\mathbf{\Sigma }}_{\alpha \mathrm{\ell }0}^{-1}{\mathit{\alpha }}_{\mathrm{\ell }0}+\sum _{K_i=k}{\mathbf{x}}_i{\eta }_{\mathit{ik}\mathrm{\ell }}\right)$$E29

with ${\eta }_{\mathit{ik}\mathrm{\ell }}={\kappa }_{i\mathrm{\ell }}+C_{\mathit{ik}\mathrm{\ell }}{\omega }_{i\mathrm{\ell }}$. Consequently, drawing ${\mathit{\alpha }}_{k\mathrm{\ell }}$is accomplished by first drawing ${\omega }_{i\mathrm{\ell }}$from the Pólya gamma distribution and then drawing ${\mathit{\alpha }}_{k\mathrm{\ell }}$from the normal distribution. The draw of ${\omega }_{i\mathrm{\ell }}$is a little intractable since its density function involves the infinite sum. By taking advantage of series sampling method [54], Polson et al. [53] devised a rejection algorithm for generating observations from such type of distribution. Their method can be adapted to draw ${\omega }_{i\mathrm{\ell }}$, see also [55].

4. A real example

In this section, a small portion of cocaine use data is analyzed to illustrate the practical value of the proposed methodology. The data are obtained from the 322 cocaine use patients who were admitted in 1988–89 to the West Los Angeles Veterans Affairs Medical Center. The original data set is made up of 68 measurements in which 17 items were assessed at four unequally spanned time points. In this study, we mainly focus on the measurements 1 year after treatment and ignore the initial effects at the baseline. The measurements cover the information on the cocaine use, treatment received, psychological problems, social status, employments, and so on. Among them, the measurement “cocaine use per month” (denoted by CC) plays a critical role since it identifies the severity of cocaine use of patients and therefore is treated as the dependent response. The CC is originally measured by 0–30 points but suffered from small portion of fractions. We identify CC/30 as the fraction response in [0,1]. In view of that the missing data are presented, we delete the subjects with missing values. The total sample size is 228. A primary analysis shows that CC/30 has excessive zeros and ones. Figure 2 gives the histograms of CC/30 and their fractional values in (0,1) via logistic transformation. It can be seen clearly that there is a large number of zeros and unities accumulated on the boundaries. The proportions of zeros and unities are about 15 and 4%, respectively. Moreover, panel (b) in Figure 2 indicates that single parametric model may be unappreciate for fitting the continuous valued variable.

Figure 2.

Plots of CC in cocaine use data: (a) Histograms of CC/30; and (b) histograms of CC/30 on logistic transformation conditional on CC/30 in (0,1).

To explore the effects of exogenous factors on the cocaine use, the following measurements are selected as the explanatory variables: the occupational status of a patient ($x_1$). This is a binary indicator: 1 for employment and 0 for non-employment; the level of technical proficiency of patients engaged in work ($x_2$): scaled on 0–4 points and the patient’s lifestyle ($x_3$) with five-point scale. To unify the scales, all covariates are standardized. However, a preliminary analysis shows that there exists strong multiple collinearity among these covariates. The minimum eigenvalue of sample covariance matrix equals to 0.06284, which approaches zero. We remove such collinearity by implementing principle component analysis (PCA) and treat the scores of the first two components (still denoted by $x_1$and $x_2$) as our explanatory variables. These two principle components can be interpreted as the levels related to the patients’ occupation and their live life.

To formulate a two-part model for the observed responses, we identity CC_i/30 with $d_i$and $z_i$, where $d_i$is the three-category indicator indicating the state of cocaine use after one year treatments: quitting cocaine successfully (state 1), insisting on cocaine use every day in a month (state 2) and taking the cocaine occasionally (state 3); $z_i$is the intensity variable representing the numer of days of cocaine use in a month. We assess the effects of exogenous factors $x_1$and $x_2$on the cocaine use via Eqs. (16) and (17), respectively.

We proceed data analysis by first fitting data to the $K$-component mixture two-part models with $K=1,2,\cdots ,6$. The model fits are assessed via AIC, AICc, and BIC, which are defined as $-2logp\left(\mathbf{Y}|{\widehat{\mathit{\theta }}}_K\right)$penalized by $2d_K$, $2n\left(d_K+1\right)/\left(n-d_K-2\right)$_, and $d_{K}logn$respectively, where ${\widehat{\theta }}_K$is the MLE of ${\mathit{\theta }}_K$and $d_K$is the dimension of unknown parameters under the model $K$. In view of that the Bayesian estimates and the ML estimates are close to each other, we replace the ML estimates by their Bayesian counterparts in evaluating AIC, AICc, and BIC. For computation, we take $\alpha =n^{-1}$, $n^0$, $n^1$, and $n^2$in Eq. (7), which represents our knowledge about $\mathit{\pi }$a prior. Note that for large value of $\alpha$, the Dirichlet distribution places most of the mass on its center and the prior Eq. (7) tends to be informative. However, for small $\alpha$, the Dirichlet distribution concentrates the mass on the boundaries of sampling space and the distribution tends to be degenerated and sparse. As a result, some components in $\mathit{\pi }$reduces to zeros. When $\alpha =1$, $D_K\left(\alpha ,\cdots ,\alpha \right)$becomes an uniform distribution on the simplex ${\mathbb{S}}_K$. For the inputs of the hyper-parameters involved in the priors Eq. (20), we take ${\mathit{\alpha }}_{0\mathrm{\ell }}={\mathit{\gamma }}_0={\mathbf{0}}_3$, ${\mathbf{\Sigma }}_{\alpha \mathrm{\ell }0}={\mathbf{\Sigma }}_{\gamma 0}=100{\mathbf{I}}_3$, ${\alpha }_{\gamma 0}=2.0$and ${\beta }_{\gamma 0}=2.0$. These values ensure the priors Eq. (20) to be inflated enough and represent the weak information on the parameters.

The relabeling MCMC algorithm described in Section 2 is implemented to draw observations from the posterior. The convergence of algorithm is monitored by plotting the traces of estimates against iterations under three starting values. Figure 3 presents the values of EPSR of unknown parameters against the number of iterations under three different starting values with $K=2$. It shows that the convergence of the proposed algorithm is fast and the values of EPSR are less than 1.2 in less than 1000 iterations. Hence, 3000 observations, after removing the initial 2000 iterations, are collected for calculating AIC, AICc, and BIC. The resulting summary is given in Table 1.

Figure 3.

Plots of values of EPSR of estimates of unknown parameters against the number of iterations under three different staring values in the cocaine use example:$K=2$.

	Model	$\alpha =1/n$	$\alpha =n^0$	$\alpha =n$	$\alpha =n^2$
AIC	$K=1$	921.3887	–	–	–
	$K=2$	923.0580	907.4485	901.9474	901.4380
	$K=3$	929.0698	926.5423	956.4945	994.5039
	$K=4$	990.7506	949.5966	1014.5477	1006.4228
	$K=5$	989.2483	971.4688	1069.1561	1037.5005
	$K=6$	1097.6853	1007.4899	1091.8049	1086.8491
AICc	$K=1$	882.4025	–	–	–
	$K=2$	885.5434	869.9339	864.4329	863.9234
	$K=3$	875.9005	873.3730	903.3253	941.3347
	$K=4$	925.3159	884.1618	949.1130	940.9880
	$K=5$	915.5836	897.8041	995.4914	963.8357
	$K=6$	1020.6483	930.4529	1014.7679	1009.8120
BIC	$K=1$	995.6821	–	–	–
	$K=2$	995.0742	979.4647	973.9637	973.4542
	$K=3$	1038.8088	1036.2814	1066.2335	1104.2429
	$K=4$	1138.2125	1097.0585	1162.0096	1153.8846
	$K=5$	1174.4330	1156.6534	1254.3408	1222.6851
	$K=6$	1320.5928	1230.3973	1314.7124	1309.7565

Table 1.

Summary statistics of AIC, AICcc, and BIC for model selection in cocaine use data analysis.

Examination of Table 1 shows that all measures favor the model with $K=2$. This indicates that the proposed model with two groups seems to give a better fit to the data. It also indicates that large $\alpha$favors the model fit. Furthermore, we calculate the posterior predictive density estimate of $z_i$under the elected model. Results (not represented here for saving spaces) show that our method can be successful in capturing the skewness and modes of data. We also follow [56] to plot the estimated residuals ${\widehat{\delta }}_i=z_i-\widehat{\mathit{\gamma }}{\mathbf{x}}_i^T$and find that these plots lie within two parallel horizontal lines that are centered at zero, with nonlinear or quadratic trends detected. This roughly indicates that the proposed linear model Eq. (18) is adequate.

Table 2 presents the estimates of unknown parameters associated with corresponding standard deviation (SD) estimates under $K=2$. Based on Table 2, we can find the following facts: (i) for Part one, we observe that except for ${\widehat{\alpha }}_{23}$, the Bayesian estimates of unknown parameters within two clusters have the same signs but their magnitudes are more different. For example, the estimate of ${\alpha }_{11}$within Cluster one is −1.540 with SD 0.587 while equals to −0.732 with SD 0.481 within Cluster two. This indicates that the baselines of logits Eq. (16) exist obvious difference. For ${\alpha }_{23}$, the estimates between two clusters have the opposite signs. Recall that ${\alpha }_{23}$quantifies the magnitude of effects of live life on the probability $\mathbb{P}\left(d_i=2\right)$over $\mathbb{P}\left(d_i=3\right)$on log scale. This shows that increasing the level of live life will lead to an opposite effect among two clusters; (ii) for Part two, although all the estimates within two clusters have the same signs but the levels of effects among them are obviously different. The estimates of ${\gamma }_1$is −2.779 with SD 0.144 in the cluster one and attains −0.490 associated with SD 0.215 in the Cluster two. This indicates that the baseline of cocaine use in Cluster one is 50 times as much as that in Cluster two; and (iii) investigation of the estimate of $\tau$also indicates that there exists the different amount of the fluctuation among two clusters.

Para.	Component I		Component II
	Est.	SD.	Est.	SD
${\alpha }_{11}$	−1.540	0.587	−0.732	0.481
${\alpha }_{12}$	0.150	0.317	0.604	0.322
${\alpha }_{13}$	0.261	0.703	0.188	0.601
${\alpha }_{21}$	−1.337	0.480	−1.059	0.545
${\alpha }_{22}$	−0.166	0.355	−0.229	0.418
${\alpha }_{23}$	0.232	0.378	−0.184	0.411
${\gamma }_1$	−2.779	0.144	−0.490	0.215
${\gamma }_2$	−0.029	0.080	−0.011	0.154
${\gamma }_3$	0.087	0.144	0.179	0.240
$\tau$	0.674	0.150	0.924	0.234

Table 2.

Summary statistics for the Bayesian estimates of unknown parameters in the cocaine use data.

5. Discussion

This chapter introduces a general Bayesian model-based clustering procedure for the regression model and proposed a Bayesian method for assessing the heterogeneity of fractional data within the mixture of two-part regression model framework. The heterogeneous fractional data arise mainly from two resources: one is that the excessive zeros and ones are accumulated upon the boundaries, and the other is that the underlying population may consist more than one components. For the first issue, we propose a novel two-part model, in which a three-category multinomial regression is suggested to model the occurrence probabilities of each part, and a conditional normal linear regression is used to fit the continuous positive values on logit scale. Such formulation is more appealing since it can ensure the probabilities on each part to be consistent and and at the same time maintains the coherent association across parts. For the second problem, we resort to the finite mixture model in which the cluster-specific components are specified via two-part model. MCMC sampling method is adopted to carry out the posterior analysis. The number of clusters and the configuration of observations are addressed based on the simulated observations from the posterior. We illustrate the proposed methodology in the analysis of cocaine use data.

When interest is concentrated upon the estimates, model identification is surely an important issue since it involves whether or not the estimates of component-specific quantities are meaningful. For a finite mixture model, model identification mainly stems from the label switching, in which the likelihood and the posterior are invariant under label permutation. Many efforts have devoted to alleviating such indeterminacy. Among them, parameters’ constraints may be the most popular choice. However, an unappreciated constraint fails to deal with the label switching. In this case, one can follow the routine in Frühwirth-Schnatter [18] and implement random permutation sampling to find the suitable identifiability constraints. The random permutation sampler is similar to the unconstrained MCMC sampling but only at each sweep, the labels $\left\{1,\cdots ,K\right\}$are randomly permutated. The permutation aims to deliver a sample that explores the whole unconstrained parameter space and jumps between the various labeling subspaces in a balanced fashion. The output of such balanced sample can help us to find a suitable identifiability constraint. A more detailed discussion on model identification in the mixture context can be referred to, for example, [18, 57]. Instead, we resort to the relabeling algorithm for simplicity. Compared with the random permutation sampling, the relabeling method requires implementing MCMC samplng only once, thus saving the computation cost.

The methodology developed in this chapter can be extended to the case where latent factors are included to identify the unobserved heterogeneity due to some fixed convariates absent. Another possible extension is to establish a dynamic LVM, wherein model parameters vary across times. These issues may raise theoretical and computational challenges and therefore require further investigation.

Acknowledgments

The work presented here was fully supported by grant from the National Natural Science Foundation of China (NNSF 11471161). The authors are thankful to Professor Xin-Yuan Song, the Chinese University of Hong Kong, for using her cocaine use data in the real example.

Conflict of interest

The authors have no conflicts of interest to disclose.