Bayesian Feature Extraction for Two-Part Latent Variable Model with Polytomous Manifestations

1. Introduction

Semi-continuous data, characterized by excessive zeros, are very common in the fields of social science and economics. A typical example is given by [8] in the analysis of medical expenditures, in which the zeros correspond to a subpopulation of patients who do not use health services, while the positive values describe the actual levels of expenditures among users. In understanding such type of data structure, two-part model [7] is a widely appreciated statistical method. The basic assumption on two-part model is that the overall model is consisted of two processes: one binary process (Part one) and one continuous positive-valued process (Part two). The binary process, usually formulated within the logistic or probit regression model, is used to indicate whether the items are responded or not, while the continuous process, conditioning on the binary process, is used to describe what the actual levels of responses are (see, e.g., [32]). By combining two processes into one, two-part model provides a unified and flexible way in describing various relevance underlying semi-continuous data. Now, two-part model has been widely used in the health service [30,43,44], the medical expenditures [8,29,40,48], the household finance [3], the substance use study [33,50] and the genome analysis [4].

The traditional two-part model usually formulates the exogenous explanatory factors as fixed and observed. However, in the real applications especially in the social survey, many unobserved/latent and random factors also have important impacts on the outcome variable. This fact is revealed by [25] in the study of children’s aggressive behavior. [25] noted that two factors, the propensity to engage in aggressive behavior and the propensity to have high aggressive activity levels, had significant influence on the children’s aggressive behavior. They incorporated such two latent factors into analysis and established a two-component-two-part mixture model to identify the heterogeneity of population; [12] noticed that in China, the finance literacy of a family had a nonignorable influence on the desire to holding finance debts, and also affected the amount of finance debt being held. They suggested conducting a joint analysis of latent factor and observed covariates in two-part regression model. The latent factor is further manifested by the multiple binary measurements via factor analysis model. [49] incorporated two-part regression model into the general latent variable model framework and analyzed the internal relationships between multiple factors longitudinally. These methods have brought significant attention to the two-part model in behavioral science, economics, psychology, and medicine in recent years, see for example, [4,20,51] and references therein for further developments of two-part model.

In the analysis of semi-continuous data, an important issue is to determine which explanatory factor is helpful in improving model fit. This issue is especially true when the number of exogenous factors is large since the commonly used forward and backward regression procedure is extremely time-consuming. Now, lasso and its extensions [11,14,21,22,26,47,54,55] have been the most commonly used methods for the feature extraction. A typical feature of these methods is to put some suitable penalties on the coefficients and shrink many coefficients to zero, thus performing variable selection. Recently, these penalization/regularizatio approaches have been applied widely for prediction and prognosis (see for example, [52,53]). Though more appealing, the lasso-type regularization also suffers some limitations, e.g., most contributions are developed within the frequency framework and the performance heavily depends on the large sample theory. It also readily leads to computational difficulty in the analysis of mixed data. An alternative for the variable selection is conducted within the Bayesian framework. Statisticians have introduced hierarchical models with mixture spike-and-slab priors that can adaptively determine the amount of shrinkage [15,16]. The spike-and-slab prior is the fundamental basis for most Bayesian variable selection approaches, and has proved remarkably successful [6,15,16,23,24,31,37]. Recently, Bayesian spike-and-slab priors have been applied to predictive modeling and variable selection in large scale genomic studies, see [45] for a simple review. Nevertheless, model selection has never been considered in the two-part regression model with latent variables. In this study, we introduce spike and slab model and Bayesian lasso into two-part latent variable model, which is the first attempt for this model.

Our formulation is more along with the lines of spike and slab bimodal prior in [24] and Bayesian Lasso in [34]. We formulate the problem by specifying a normal distribution with mean zero to the regression coefficient or factor loading of interest. The probability of a related variable being excluded or included is governed by the variance. To model the shrinkage of coefficients properly, we consider two schemes for the variance parameter: one is the two-point mixture model with one component located at the point close to zero and the other component situated at the point far away zero. The proportion is governed by a beta-distribution with suitable hyperparameters. Another scheme is along with the Bayesian lasso in which the variance is specified via a gamma distribution scaled by the penalty parameters. Two schemes are unified into a hierarchical mixture model. Within the Bayes paradigm, we developed a fully Bayesian selection procedure for the two-part latent variable model. We resort to the Markov Chains Monte Carlo sampling method. Gibbs sampler is used to draw observations from the posterior. We obtained all full conditionals. Posterior analysis is carried out based on the simulated observations. We investigate the performance of the proposed methods via simulation study and a real example. Our empirical results show that two schemes results in the similar results in the variable selection but SS with suitable hyperparameters slightly outperforms over BaLsso in the correct rate.

The remainder of this paper is organized as follows. Section 2 introduces the proposed model for the semi-continuous data with latent variables. Section 3 develops the MCMC algorithm for the proposed model. Bayesian inference procedures including parameters estimation and model assessment are also presented in this section. In Section 4, we illustrate the use of our proposed model with household finance debt data and present the results of a simulation study to assess the performance of our proposed methodology. Section 5 concludes the paper with a discussion. Some technical details are given in the Appendix.

2. Model description

In Section 2.1, a basic formulation for analyzing semi-continuous data with latent variables is presented. Section 2.2 presents a Bayesian procedure for the feature extraction.

2.1. Two-part latent variable model

Suppose that for $i=1,2,\dots ,n$, $s_i$ is a semicontinuous outcome variable which takes value in $[0,\infty )$; ${\mathbf{x}}_i$ is a generic vector consisted of r fixed covariates representing the collection of observed explanatory factors of interest. We assume that each $x_{ij}$ in ${\mathbf{x}}_i$ is standardized in the sense ${\sum }_{i=1}^n{x}_{ij}=0$ and ${\sum }_{i=1}^n{x}_{ij}^2=1$ for $j=1,\dots ,r$. Moreover, we include m letent/unobserved variables ${\mathit{\omega }}_i={({\omega }_{i1},\dots ,{\omega }_{im})}^T$ into analysis to account for unobserved heterogeneity of responses. Conceptually, these latent variables can be the covariates that are not directly observed or the synthesization of some highly correlated explanatory items suffering from the noisy. Inclusion of latent variables can improve model fits and strengthen the power of model interpretation, see [39] for more discussions on the latent variables in a general setting. To deal with the spike of $s_i$ at zero, we follow the common routine in literature (see for example, [33,48]) and identify $s_i$ with two surrogate variables: $u_i=I\{s_i>0\}$ and $z_i=log(s_i^{+})$, where $I(A)$ denotes the indicator function of set A and $a^{+}$ represents the positive part of a. That is, we separate the whole dataset into two parts: one part is the binary dataset which corresponds to the response-to-nonresponse indicators of subject and the other part is the logarithm of positive values. Our interest focuses on the exploration of effects of exogenous factors on two parts.

We assume that $u_i$ and $z_i$ satisfy the following sampling models:

$$\begin{array}{ccc}& & p(u_i|{\mathbf{x}}_i,{\mathit{\omega }}_i)=\frac{exp(u_i{\eta }_i^u)}{1+exp({\eta }_i^u)},\hfill \\ & & {\eta }_i^u=\alpha +{\mathit{\beta }}_x^T{\mathbf{x}}_i+{\mathit{\beta }}_{\omega }^T{\mathit{\omega }}_i,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}& & p(z_i|u_i=1,{\mathit{\omega }}_i)=N({\eta }_i^z,{\sigma }^2),\hfill \\ & & {\eta }_i^z=\gamma +{\mathit{\psi }}_x^T{\mathbf{x}}_i+{\mathit{\psi }}_{\omega }^T{\mathit{\omega }}_i,\hfill \end{array}$$

(2)

in which $\alpha$ and $\gamma$ are the intercept parameters, ${\mathit{\beta }}_x$ and ${\mathit{\psi }}_x$ are the vectors of regression coefficients, and ${\mathit{\beta }}_{\omega }$ and ${\mathit{\psi }}_{\omega }$ are the vectors of factor loadings; ${\sigma }^2$ is the scale and ‘T’ is the transpose operator of vector or matrix; For compactness, we write $\mathit{\beta }={({\mathit{\beta }}_x^T,{\mathit{\beta }}_{\omega }^T)}^T$ and $\mathit{\psi }={({\mathit{\psi }}_x^T,{\mathit{\psi }}_{\omega }^T)}^T$ and treat ${\mathbf{w}}_i={({\mathbf{x}}_i^T,{\mathit{\omega }}_i^T)}^T$ as the complete explanatory variables.

The involvement of latent variables apparently complicates the model. It readily results in model identification problem [2,27]. This is especially true when the dimension of ${\mathit{\omega }}_i$ is high. In this case, any auxiliary information is required to manifest ${\mathit{\omega }}_i$ further. Among various-easy-constructs, we consider latent variable (LV, [2,27]) approach. A basic assumption on LV approach is that there exist, say p manifestations ${\mathbf{y}}_i={(y_{i1},\cdots ,y_{ip})}^T$, of which each $y_{ij}$ may be continuous, counted or categorical, and assuming that they satisfy the following link equation

$$\begin{array}{c}\hfill F({\mathbf{y}}_i,{\mathit{\omega }}_i,{\mathit{ϵ}}_i,\mathit{\phi })=0,\end{array}$$

(3)

where F is a known and fixed link function, ${\mathit{ϵ}}_i$ is the vector of errors used to identity the idiosyncratic part of ${\mathbf{y}}_i$ that can not be explained by ${\mathit{\omega }}_i$, and $\mathit{\phi }$ is the vector of unknown parameters used to quantify the uncertainty of model. The information about ${\mathit{\omega }}_i$ is manifested by ${\mathbf{y}}_i$ via F. In this paper, in view of the real applications, we consider p ordered categorical variable ${\mathbf{y}}_i={(y_{i1},\cdots ,y_{ip})}^T$, of which $y_{ij}$ takes value in $\{0,1,\dots ,c_j\}(c_j>1)$ and satisfies the following link model:

$$\begin{array}{c}\hfill y_{ij}={\ell }_j\mathrm{if}{\delta }_{j{\ell }_j}<y_{ij}^{\ast }\le {\delta }_{j,{\ell }_j+1},\end{array}$$

(4)

where ${\delta }_{j0}<{\delta }_{j1}<\cdots <{\delta }_{j{c}_j}<{\delta }_{j,c_j+1}$ are the threshold parameters with ${\delta }_{j0}=-\infty$ and ${\delta }_{j,c_j+1}=+\infty$, and ${\mathbf{y}}_i^{\ast }={(y_{i1}^{\ast },\cdots ,y_{ip}^{\ast })}^T$ is the vector of latent responses satisfying the factor analytic model:

$$\begin{array}{cc}\hfill & {\mathbf{y}}_i^{\ast }=\mathit{\mu }+\mathbf{\Lambda }{\mathit{\omega }}_i+{\mathit{ϵ}}_i,\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\mathit{\omega }}_i\stackrel{iid.}{\sim }{N}_m[\mathbf{0},\mathsf{\Phi }],{\mathit{ϵ}}_i\sim N_p[\mathbf{0},{\mathbf{I}}_p],\mathrm{and}{\mathit{\omega }}_i\perp {\mathit{ϵ}}_i,\hfill \end{array}$$

(6)

where $\mathit{\mu }$ is a p-dimensional intercept vector, $\mathbf{\Lambda }$ is the $p\times m$-dimensional factor loading matrix, and ${\mathbf{I}}_p$ is the identity matrix of order p. We assume that conditional upon ${\mathit{\omega }}_i$, $s_i$ and ${\mathbf{y}}_i$ are independent.

We refer to the model specified by (1), (2) and (4) associated with (5) as the two-part latent variable model with polytomous responses. It provides a unified framework to explore the dependence of binary, continuous and categorical data simultaneously. The dependence between them results from the share of common factors or latent variables. If ${\mathit{\omega }}_i$ is degenerated at zeros or the factor loadings are taken as zeros, the dependence among them disappears, and the overall model reduces to the traditional two-part model and ordinal regression model.

To facilitate the efficient calculation, motivated by the key identity in [35] (see squation (2) in their seminar paper), we express the logistics model (1) as the mixture model of form

$$\begin{array}{ccc}\hfill & & \frac{exp(u_i(\alpha +{\mathit{\beta }}^T{\mathbf{w}}_i))}{1+exp(\alpha +{\mathit{\beta }}^T{\mathbf{w}}_i)}\hfill \\ & & =2^{-1}exp({\kappa }_i(\alpha +{\mathit{\beta }}^T{\mathbf{w}}_i)){\int }_0^{\infty }exp\left\{-\frac{u_i^{\ast }}{2}{(\alpha +{\mathit{\beta }}^T{\mathbf{w}}_i)}^2\right\}{p}_{PG}(u_i^{\ast })\mathrm{d}{u}_i^{\ast },\hfill \end{array}$$

(7)

where ${\kappa }_i=u_i-0.5$, and $p_{PG}(u)$ is the standard Pólya-Gamma probability density function. Assuming that we introduce auxiliary variables $u_i^{\ast }$ and augment them with $u_i$, then equation (1) can be considered as the marginal density of the joint distribution

$$\begin{array}{c}\hfill p(u_i,u_i^{\ast }\mid {\mathbf{x}}_i,{\mathit{\omega }}_i)=2^{-1}exp\left\{{\kappa }_i{\eta }_i^u-\frac{u_i^{\ast }}{2}{\eta }_i^{u2}\right\}{p}_{PG}(u_i^{\ast }).\end{array}$$

(8)

Note that the exponential part in the brackets is the kernel of normal density function with respect to ${\eta }_i^u$. Hence, it admits conjugate full-conditional distributions for all regression coefficients, factor loadings and factor variables, leading to a straightforward Bayesian computation.

Let $\mathbf{U}={\left\{u_i\right\}}_{i=1}^n$, $\mathbf{Z}={\left\{z_i\right\}}_{i=1}^n$, and $\mathbf{Y}={\left\{{\mathbf{y}}_i\right\}}_{i=1}^n$ be the sets of observed variables; We write $\mathsf{\Omega }={\left\{{\mathit{\omega }}_i\right\}}_{i=1}^n$ for the collect of factor variables, and write ${\mathbf{U}}^{\ast }={\left\{u_i^{\ast }\right\}}_{i=1}^n$, ${\mathbf{V}}^{\ast }={\left\{v_i^{\ast }\right\}}_{i=1}^n$, ${\mathbf{Y}}^{\ast }={\left\{{\mathbf{y}}_i^{\ast }\right\}}_{i=1}^n$ for the sets of latent response variables. The complete-data likelihood is given by

$$\begin{array}{cc}\hfill & p(\mathbf{U},{\mathbf{U}}^{\ast },\mathbf{Z},{\mathbf{V}}^{\ast },\mathbf{Y},{\mathbf{Y}}^{\ast },\mathsf{\Omega }|\mathit{\theta })\hfill \\ \hfill =& p(\mathbf{U},{\mathbf{U}}^{\ast }|\mathsf{\Omega },\alpha ,\mathit{\beta })p(\mathbf{Z},{\mathbf{V}}^{\ast }|\mathbf{U},\mathsf{\Omega },\gamma ,\mathit{\psi },{\sigma }^2)p(\mathbf{Y}|{\mathbf{Y}}^{\ast },\mathit{\delta })p({\mathbf{Y}}^{\ast }|\mathsf{\Omega },\mathit{\mu },\mathbf{\Lambda })p(\mathsf{\Omega }|\mathsf{\Phi })\hfill \\ \hfill =& \prod _{i=1}^{n}exp\left\{{\kappa }_i{\eta }_i^u-\frac{1}{2}{u}_i^{\ast }{({\eta }_i^u)}^2\right\}{p}_{PG}(u_i^{\ast }|1,0)\hfill \\ \hfill & \times \prod _{i\in I}\frac{1}{\sqrt{2\pi }\sigma }exp\left\{-\frac{1}{2{\sigma }^2}{(z_i-{\eta }_i^z)}^2\right\}\hfill \\ \hfill & \times \prod _{i=1}^n\prod _{j=1}^p\sum _{\ell =0}^{c_j}I\{{\delta }_{j\ell }<y_{ij}^{\ast }\le {\delta }_{j,\ell +1},y_{ij}=\ell \}\hfill \\ \hfill & \times \prod _{i=1}^n\prod _{j=1}^p\frac{1}{\sqrt{2\pi }}exp\left\{-\frac{1}{2}{(y_{ij}^{\ast }-{\mu }_j-{\mathbf{\Lambda }}_j^T{\mathit{\omega }}_i)}^2\right\}\hfill \\ \hfill & \times \prod _{i=1}^n\frac{1}{{(\sqrt{2\pi })}^m{\left|\mathsf{\Phi }\right|}^{1/2}}exp\left\{-\frac{1}{2}tr\left[{\mathsf{\Phi }}^{-1}{\mathit{\omega }}_i{\mathit{\omega }}_i^T\right]\right\}.\hfill \end{array}$$

where $I=\{i:u_i=1\}$ is the set of indices, $\mathit{\delta }=\left\{{\delta }_{j\ell }\right\}$ is the set of threshold parametes, and $\mathit{\theta }=\{\alpha ,\mathit{\beta },\gamma ,\mathit{\psi },{\sigma }^2,\mathit{\mu },\mathbf{\Lambda },\mathsf{\Phi },\mathit{\delta }\}$ is the vector of unknown parameters. For the moment, we assume ${\theta }_j$ in $\mathit{\theta }$ all are free.

2.2. Bayesian feature selection

Generally speaking, regression variables ${\mathbf{x}}_i$ and factor variables ${\mathit{\omega }}_i$ may have not impacts on the $u_i$ and $z_i$ simultaneously, and some redundant variables may exist. The presence of redundant variables not only decreases the model fit but also weakens thepower of model interpretation. Therefore, it is necessary to determine which regression coefficient or factor loading is significantly away from zero. In the context of frequency statistics, this issue is generally tackled out via stepwise regression, in which each variable is decided to be exclude or included according to the model fit. However, the situation becomes complex when the number of independent variables is large. In this paper, we pursuit a bayesian variable selection procedure. To this end, we follow [34] and assume

$$\begin{array}{cc}\hfill & \mathit{\beta }\sim N_q(\mathbf{0},diag\left\{{\gamma }_{\beta k}^2\right\}),\mathit{\psi }\sim N_q(\mathbf{0},{\sigma }^{2}diag\left\{{\gamma }_{\psi k}^2\right\}),\hfill \end{array}$$

(9)

in which we use $diag\left\{a_k\right\}$ to represent a diagonal matrix with the $k^{th}$ diagonal element $a_k$. That is, we assume that each ${\beta }_k$ in $\mathit{\beta }$ (${\psi }_k$ is similar) is centered at zero (or equivalently each $w_{ik}$ is excluded from ${\mathbf{w}}_i$) but with the probability governed by the variance ${\gamma }_{\beta k}^2$. If ${\gamma }_{\beta k}^2$ is close to zero, then the probability of ${\beta }_k$ taking zero increases, and $w_{ik}$ tends to be excluded; conversely, if ${\gamma }_{\beta k}^2$ is large, then the probability of ${\beta }_k$ being zero is small and $w_{ik}$ tends to be maintained. As a result, the value of ${\gamma }_{\beta k}^2$ plays a key role in determining whether $w_k$ is relevant to be selected in Part one. With this in mind, a reasonable assumption on ${\gamma }_{\beta k}^2$ and ${\gamma }_{\psi k}^2$ is that:

$$\begin{array}{ccc}& & {\gamma }_{\beta k}^2\stackrel{ind.}{\sim }(1-w_{\beta }){\delta }_{{\nu }_{\beta 0}{\eta }_{\beta k}^2}(·)+w_{\beta }{\delta }_{{\eta }_{\beta k}^2}(·),\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& & {\gamma }_{\psi k}^2\stackrel{ind.}{\sim }(1-w_{\psi }){\delta }_{{\nu }_{\psi 0}{\eta }_{\psi k}^2}(·)+w_{\psi }{\delta }_{{\eta }_{\psi k}^2}(·),\hfill \end{array}$$

(11)

where ${\delta }_a(·)$ is the Dirac measure concentrated at point a, $w_{\beta }$ is the random weight used to measure the similarity between ${\gamma }_{\beta k}^2$ and ${\eta }_{\beta k}^2$, and ${\eta }_{\beta k}^2$ is the hyperparameter used to represent how far ${\beta }_k$ is away from zero or slab; ${\nu }_{\beta 0}$ is a previously specified small positive values used to identity the ‘spike’ of ${\beta }_k$ at zero. In other words, every ${\gamma }_{\beta k}^2$ is assumed to be equal to ${\eta }_{\beta k}^2$ with probability $w_{\beta }$ and equal to ${\nu }_{\beta 0}{\eta }_{{\beta }_k}^2$ with probability $1-w_{\beta }$. This is also true for $w_{\psi }$, ${\eta }_{\psi k}$ and ${\nu }_{\psi 0}$. To model $w_{\beta }$ and $w_{\psi }$ properly, we assign the following beta distributions to them

$$\begin{array}{c}\hfill p(w_{\beta }|a_{\beta },b_{\beta })=Beta(a_{\beta },b_{\beta }),p(w_{\psi }|a_{\psi },b_{\psi })=Beta(a_{\psi },b_{\psi }),\end{array}$$

(12)

where $a_{\beta }$, $a_{\psi }$, $b_{\beta }$ and $b_{\psi }$ are the hyperparameters used to control the shape of beta density, that is, to determine the magnitude of weights in $(0,1)$. For example, if $a_{\beta 1}$ in equation (12) is small and $b_{\beta 1}$ is large, then equation (12) encourages $w_{\beta }$ to take small value with high probability. In contrast, it follows from $1-Beta(a_{\beta },b_{\beta })=Beta(b_{\beta },a_{\beta })$ that large $a_{\beta 1}$ and small $b_{\beta 1}$ encourage $w_{\beta }$ to take large value in $(0,1)$. In the case that $a_{\beta 1}=b_{\beta 1}=1.0$, equation (12) reduces to the uniform distributions on $(0,1)$. In this case, every value in $(0,1)$ is possible for $w_{\beta }$ with the same probability. In the real applications, if no informative can be available, one can assign the values to them to ensure the beta distribution to be inflated enough.

Finally, to measure the magnitudes of ‘slap’ in the distributions of ${\beta }_k$ and ${\psi }_k$, we specify gamma distributions for ${\eta }_{\beta k}^{-2}$ and ${\eta }_{\psi k}^{-2}$, or equivalently,

$$\begin{array}{c}\hfill {\eta }_{\beta k}^2|a_{\beta 1},a_{\beta 2}\stackrel{iid.}{\sim }IG({\tau }_{\beta 0},{\zeta }_{\beta 0}),{\eta }_{\psi k}^2|a_{\psi 1},a_{\psi 2}\stackrel{iid.}{\sim }IG({\tau }_{\psi 0},{\zeta }_{\psi 0}),\end{array}$$

(13)

where ‘$IG(a,b)$’ denotes the inverse-gamma distribution with mean $b/(a-1)$ for $a>1$ and variance $b^2/({(a-1)}^2(a-2))$ for $a>2$; ${\tau }_{\beta 0}$, ${\zeta }_{\beta 0}$, ${\tau }_{\psi 0}$ and ${\zeta }_{\psi 0}$ are the hyperparameters which are treated to fixed and known. Similarly, one can set the values to them to ensure the inverse gamma distribution inflated enough. For example, we can follow the routine in [24] in the ordinary regression analysis, and set ${\tau }_{\beta 0}={\tau }_{\psi 0}=1.0$ and ${\zeta }_{\beta 0}={\zeta }_{\psi 0}=0.05$ to obtain dispersed priors.

Note that equations (10) and (11) can be formulated as hierarchy as follows: for $k=1,\dots ,q$,

$$\begin{array}{ccc}& & {\gamma }_{{\beta }_k}^2=f_{\beta k}{\eta }_{\beta k}^2,{\gamma }_{\psi k}^2=f_{\psi k}{\eta }_{\beta k}^2,\hfill \end{array}$$

$$\begin{array}{ccc}& & f_{\beta k}|{\nu }_{\beta 0},{\omega }_{\beta }\stackrel{iid.}{\sim }(1-w_{\beta }){\delta }_{v_{\beta 0}}(·)+w_{\beta }{\delta }_1(·),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & f_{\psi k}|{\nu }_{\psi 0},w_{\psi }\stackrel{iid.}{\sim }(1-w_{\psi }){\delta }_{v_{\psi 0}}(·)+w_{\psi }{\delta }_1(·),\hfill \end{array}$$

(15)

where $f_{\beta k}$ and $f_{\psi k}$ are the latent binary variables respectively. Such a formulation aims to separate ${\eta }_{\beta k}^2$ and ${\eta }_{\psi k}^2$ from the distributions (10) and (11) to facilitate posterior sampling.

It is instructive to compare the proposed method to the Bayesian lasso [34], in which the variance parameters ${\gamma }_{\beta k}^2$ and ${\gamma }_{\psi k}^2$ in equation (9) are specified via exponential distributions as follows:

$$\begin{array}{cc}\hfill & p({\mathit{\gamma }}_{\beta }^2\mid {\mathit{\lambda }}_{\beta }^2)=\prod _{k=1}^q\frac{{\lambda }_{\beta k}^2}{2}exp(-{\lambda }_{\beta k}^2{\gamma }_{\beta k}^2/2),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill & p({\mathit{\gamma }}_{\psi }^2\mid {\mathit{\lambda }}_{\psi }^2)=\prod _{k=1}^q\frac{{\lambda }_{\psi k}^2}{2}exp(-{\lambda }_{\psi k}^2{\gamma }_{\psi k}^2/2),\hfill \end{array}$$

(17)

where ${\mathit{\lambda }}_{\beta }^2={({\lambda }_{\beta 1}^2,\dots ,{\lambda }_{\beta q}^2)}^T$ and ${\mathit{\lambda }}_{\psi }^2=({\lambda }_{\psi 1}^2,\dots ,$${\lambda }_{\psi q}^2{)}^T$${\lambda }_{\psi k}^2$ are the shrinkage/penality parameters used to control the amount of shrinkage of ${\beta }_k$ and ${\psi }_k$ toward zero.

Modeling ${\gamma }_{\beta k}^2$ and ${\gamma }_{\psi k}^2$ like equations (16) and (17) lead to marginal distributions of ${\beta }_k$ and ${\psi }_k$ as the laplace distributions with location zero and scale ${\lambda }_k$. The penalty parameters ${\lambda }_{\beta k}^2$ and ${\lambda }_{\psi k}^2$ are rather crucial in determining the amount of shrinkage of parameters. Figure 1 presents the plots of densities of Laplace distribution $LA(\lambda )(\lambda >0)$ across various choices of $\lambda$. It can be seen that the larger the value of $\gamma$, the more kurtosis the density, indicating much penalty on the regression coefficient.

Due to their key role in equations (16) and (17), for ${\lambda }_{\beta }^2$ and ${\lambda }_{\psi }^2$, we assign gamma priors to them, i.e.,

$$\begin{array}{c}\hfill p({\mathit{\lambda }}_{\beta }^2)=\prod _{k=1}^{q}p({\lambda }_{\beta k}^2)=\prod _{k=1}^{q}Ga(a_{k0},b_{k0}),\end{array}$$

(18)

$$\begin{array}{c}\hfill p({\mathit{\lambda }}_{\psi }^2)=\prod _{k=1}^{q}p({\lambda }_{\psi k}^2)=\prod _{k=1}^{q}Ga(c_{k0},d_{k0}),\end{array}$$

(19)

where ‘$Ga(\nu ,\lambda )$’ denotes the gamma distribution with mean $\nu /\lambda$. As the previous discussions, the values of $a_{k0}$, $b_{k0}$, $c_{k0}$ and $d_{k0}$ should be selected with care since they relate the shrinkages directly. Similar to that in (13), we set $a_{k0}=c_{k0}=1$, $b_{k0}=d_{k0}=0.05$ to enhance the robustness of inference. Their procedure is also followed in our empirical study.

Let ${\mathbf{F}}_{\beta }^{\ast }=\left\{f_{\beta k}\right\}$, ${\mathbf{F}}_{\psi }^{\ast }=\left\{f_{\psi k}\right\}$, ${\mathit{\gamma }}_{\beta }^2=\left\{{\gamma }_{\beta k}^2\right\}$, ${\mathit{\gamma }}_{\psi }^2=\left\{{\gamma }_{\psi k}^2\right\}$, ${\mathit{\eta }}_{\beta }^2=\left\{{\eta }_{\beta k}^2\right\}$, ${\mathit{\eta }}_{\psi }^2=\left\{{\eta }_{\psi k}^2\right\}$. We treat ${\nu }_{\beta 0}$ and ${\nu }_{\psi 0}$ as known hyperparameters. Note that ${\mathit{\gamma }}_{\beta }^2$ and ${\mathit{\gamma }}_{\psi }^2$ are totally determined by ${\mathbf{F}}_{\beta }^{\ast }$, ${\mathbf{F}}_{\psi }^{\ast }$ and ${\mathit{\eta }}_{\beta }^2$, ${\mathit{\eta }}_{\psi }^2$. In the following, we abbreviate spike and slab bimodal prior to SS and Bayesian lasso to BaLsso.

3. Bayesian inference

3.1. Prior specification and MCMC sampling

In view of the model complexity, we consider Bayesian inference. Some priors are required to specify for unknown parameters to complete Bayesian model specification. Based on the model convention, it is natural to assume that the parameters involved in the different models are independent.

Firstly, for $\mathit{\mu }$, $\mathbf{\Lambda }$ and $\mathsf{\Phi }$, we consider the following conjugate priors:

$$\begin{array}{ccc}& & p(\mathit{\mu })=N_p({\mathit{\mu }}_0,{\mathbf{\Sigma }}_0),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & p(\mathbf{\Lambda })=\prod _{k=1}^{p}p({\mathbf{\Lambda }}_k)=\prod _{k=1}^p{N}_m({\mathbf{\Lambda }}_{0k},{\mathbf{H}}_{0k}),\hfill \end{array}$$

(21)

$$\begin{array}{ccc}& & p(\mathsf{\Phi })=IW({\rho }_0,{\mathbf{R}}_0^{-1}),\hfill \end{array}$$

(22)

where ‘$IW(\rho ,\mathbf{R})$’ denotes the inverse Wishart distribution with degrees of freedom $\rho$ and scale matrix $\mathbf{R}$ [1]; ${\mathbf{\Lambda }}_k^T$ is the $k^{th}$ row vector of $\mathbf{\Lambda }$; ${\mathit{\mu }}_0$, ${\mathbf{\Sigma }}_0(p\times p)>0$, ${\mathbf{\Lambda }}_{0k}$, ${\mathbf{H}}_{0k}(m\times m)>0$, ${\rho }_0>0$, and ${\mathbf{R}}_0(m\times m)>0$ are the hyperparameters which are treated to be fixed and known.

Secondly, for $\alpha$, $\gamma$, ${\sigma }^2$ in part one and two, we assume they are mutually independent and satisfy

$$\begin{array}{c}\hfill p(\alpha )=N({\alpha }_0,{\sigma }_{\alpha 0}^2),p(\gamma )=N({\gamma }_0,{\sigma }_{\gamma 0}^2),p({\sigma }^{-2})\sim Ga(a_0,b_0),\end{array}$$

(23)

where ${\alpha }_0$, ${\sigma }_{\alpha 0}^2$, ${\gamma }_0$, ${\sigma }_{\gamma 0}^2$ and $a_0$, $b_0$ are the fixed hyperparameters.

Finally, for the threshold parameter $\mathit{\delta }$, without loss of generality, we assume that $c_j$, the number of categories of $y_{ij}$, is invariant across the subscript j and equals to c. Moreover, we assume that $p(\mathit{\delta })={\prod }_{j=1}^{p}p({\mathit{\delta }}_j)$, where ${\mathit{\delta }}_j=({\delta }_{jk})$ is the $j^{th}$ row vector of $\mathit{\delta }$. In the following, we suppress the subscript j in ${\delta }_{jk}$ for notational simplicity and write $\mathit{\delta }$ for ${\mathit{\delta }}_j$.

Let $F_0(·)$ be any strictly monotonically increasing and differentiable function on $\mathbb{R}$ with $F_0(+\infty )=1$ and $F_0(-\infty )=0$. For example, one can take $F_0=\mathsf{\Phi }(·/{\tau }_0)$ for some ${\tau }_0>0$ or student distribution with degrees of freedom ${\nu }_0$, where $\mathsf{\Phi }(·)$ is the standard normal distribution function. To specify a prior for $\mathit{\delta }$, we follow [38] and let $p_j=F_0({\delta }_j)-F_0({\delta }_{j-1})$ for $j=1,\cdots ,c$. It is easily to show that the transformation is invertible with Jacobi determination unity. We first consider the following Dirichlet distribution for $\mathit{p}={(p_1,\cdots ,p_c)}^T$:

$$\begin{array}{c}\hfill \pi (\mathit{p})=\frac{1}{\mathcal{B}({\eta }_1,\cdots ,{\eta }_{c+1})}{p}_1^{{\eta }_1-1}\cdots p_c^{{\eta }_c-1}{(1-\sum _{\ell =1}^c{p}_{\ell })}^{{\eta }_{c+1}-1}\end{array}$$

where $\mathcal{B}({\eta }_1,\cdots ,{\eta }_{c+1})={\prod }_{j=1}^{c+1}\mathsf{\Gamma }({\eta }_j)/\mathsf{\Gamma }({\sum }_{j=1}^{c+1}{\eta }_j)$ is the multivariate beta function evaluated at ${\eta }_1,\cdots ,{\eta }_{c+1}$, and ${\eta }_j>0$. By the formula of reverse transformation, the joint distribution of $\mathit{\delta }$ is given by

$$\begin{array}{c}\hfill \pi (\mathit{\delta })=\frac{1}{\mathcal{B}({\eta }_1,\cdots ,{\eta }_{c+1})}{p}_1^{{\eta }_1-1}\cdots p_c^{{\eta }_c-1}{(1-\sum _{\ell =1}^c{p}_{\ell })}^{{\eta }_{c+1}-1}\prod _{j=1}^c{f}_0({\delta }_j)I\{{\delta }_1<\cdots <{\delta }_c\},\end{array}$$

(24)

where $f_0(x)$ is the derivative of $F_0(x)$ with respect to x. We call (24) the transformed Dirichlet prior and use it as the prior of $\mathit{\delta }$. An advantage of working with (24) is that conditional upon ${\delta }_{j-1}$ and ${\delta }_{j+1}$, the transformed distribution of ${\delta }_j$ has the beta distribution given by

$$\frac{F_0({\delta }_j)-F_0({\delta }_{j-1})}{F_0({\delta }_{j+1})-F_0({\delta }_{j-1})}|({\delta }_{j-1},{\delta }_{j+1})\sim Beta({\eta }_j,{\eta }_{j+1}),(j=1,\cdots ,c).$$

(25)

3.2. MCMC sampling

With the prior given above, the inference about $\mathit{\theta }$ is based on the posterior $p^o(\mathit{\theta }|\mathbf{U},\mathbf{Z},\mathbf{Y})$, which has no closed form. Motivated by the key idea in [46], we treat latent quantities as the missing data and argument them to the observed data to form the complete data. The statistical inference is carried out based on the complete-data likelihood. For this end, apart from $\mathsf{\Omega }$, ${\mathbf{U}}^{\ast }$ and ${\mathbf{Y}}^{\ast }$ mentioned before, we further let ${\mathbf{Q}}^{\ast }$ be the collection of latent quantities involved in the specifications of $\mathit{\beta }$ and $\mathit{\psi }$, i.e, ${\mathbf{Q}}^{\ast }=\{{\mathbf{F}}_{\beta }^{\ast }$, ${\mathbf{F}}_{\psi }^{\ast }$, ${\mathit{\eta }}_{\beta }^2$, ${\mathit{\eta }}_{\psi }^2$, $w_{\beta }$, $w_{\psi }\}$ for the SS and $\{{\mathit{\lambda }}_{\beta }^2$, ${\mathit{\lambda }}_{\psi }^2\}$ for the BaLsso. Rather than working with the posterior $p^o$ directly, we consider the following joint distribution

$$\begin{array}{c}\hfill p^{joint}(\mathsf{\Omega },{\mathbf{U}}^{\ast },{\mathbf{Y}}^{\ast },{\mathbf{Q}}^{\ast },\mathit{\theta }|\mathbf{U},\mathbf{Z},\mathbf{Y}),\end{array}$$

(26)

where $p^o$ can be considered as the marginal of $p^{joint}$. We use Markov chain Monte Carlo (MCMC, [17,19]) sampling method to simulate observations from this target distribution. In particular, Gibbs sampler is implemented to draw observations iteratively from the full conditional distributions as follows:

• draw $\mathsf{\Omega }$ from $p(\mathsf{\Omega }\mid {\mathbf{U}}^{\ast },{\mathbf{Q}}^{\ast },{\mathbf{Y}}^{\ast },\mathit{\theta },\mathbf{U},\mathbf{Z},\mathbf{V})$,
• draw ${\mathbf{U}}^{\ast }$ from $p({\mathbf{U}}^{\ast }\mid \mathsf{\Omega },{\mathbf{Y}}^{\ast },{\mathbf{Q}}^{\ast },\mathit{\theta },\mathbf{U},\mathbf{Z},\mathbf{V})$,
• draw ${\mathbf{Y}}^{\ast }$ from $p({\mathbf{Y}}^{\ast }\mid \mathsf{\Omega },{\mathbf{U}}^{\ast },{\mathbf{Q}}^{\ast },\mathit{\theta },\mathbf{U},\mathbf{Z},\mathbf{V})$,
• draw ${\mathbf{Q}}^{\ast }$ from $p({\mathbf{Q}}^{\ast }\mid \mathsf{\Omega },{\mathbf{U}}^{\ast },{\mathbf{Y}}^{\ast },\mathit{\theta },\mathbf{U},\mathbf{Z},\mathbf{V})$, and
• draw $\mathit{\theta }$ from $p(\mathit{\theta }\mid \mathsf{\Omega },{\mathbf{U}}^{\ast },{\mathbf{Y}}^{\ast },{\mathbf{Q}}^{\ast },\mathbf{U},\mathbf{Z},\mathbf{V})$.

Upon convergence, the posterior is approximated by the empirical distribution of simulated observations. The convergence of algorithm can be monitored by plotting the traces of estimates under different starting values or observing the values of EPSR [18] of unknown parameters. The technical details on implementing MCMC sampling are given in Appendix.

Simulated observations obtained from the blocked Gibbs sampler can be used for statistical inference via a straightforward analysis procedure. For example, the joint Bayesian estimates of unknown parameters can be obtained via sample averaging as follows:

$$\begin{array}{c}\hfill \widehat{\mathit{\theta }}=M^{-1}\sum _{m=1}^M{\mathit{\theta }}^{(m)},\end{array}$$

where $\{{\mathit{\theta }}^{(m)}:m=1,\cdots ,M\}$ are the simulated observations from the posterior. The consistent estimates of covariance matrices of estimates can be obtained via sample covariance matrices.

The main purpose of introducing SS and BaLsso is to screen the variable in ${\mathbf{w}}_i$. Unlike that in the frequency statistics, Bayesian variable selection does not produce the estimates $\widehat{\mathit{\beta }}$ and $\widehat{\mathit{\psi }}$ exactly equal to zero, and hence it is necessary to determine which component can be treated as zero. This can accomplished via posterior confidence intervals (PCI) of ${\beta }_j$ and ${\psi }_j$, given by

$$\begin{array}{c}\hfill P(|{\beta }_j|<c_{\alpha /2}|\mathbf{U},\mathbf{Z},\mathbf{Y})=1-\alpha ,P(|{\psi }_j|<d_{\alpha /2}|\mathbf{U},\mathbf{Z},\mathbf{Y})=1-\alpha \end{array}$$

(27)

where $\alpha$ is any previously specified value in $(0,1)$. The calculation of PCI can be achieved via Monte Carlo method. For example, let ${\beta }_j^{(k)}:k=1,\dots ,K$ be the K observations generated from the posterior distribution, then the PCI of ${\beta }_j$ with confidence level $100(1-\alpha )\%$ is given by $[{\beta }_{j,100(\alpha /2)}$, ${\beta }_{j,100(1-\alpha /2)}]$, where ${\beta }_{j,k}$ is the $k^{th}$ order statistics.

Another choice for variable determination in SS is based on the posterior probability of $f_{\beta j}=1$ and $f_{\psi j}=1$, which can be approximated by

$$\begin{array}{c}\hfill {\widehat{f}}_{\beta j}=\frac{1}{K}\sum _{k=1}^{K}I\{f_{\beta j}^{(k)}=1\},{\widehat{f}}_{\psi j}=\frac{1}{K}\sum _{k=1}^{K}I\{f_{\psi j}^{(k)}=1\},\end{array}$$

(28)

where $f_{\beta j}^{(k)}$ and $f_{\psi j}^{(k)}(k=1,\dots ,K)$ are the k observations drawn from the posterior distribution via Gibbs sampler. The variable $w_j$ is selected in part one and two if ${\widehat{f}}_{\beta j}>0.5$ and ${\widehat{f}}_{\psi j}>0.5$.

4. Simulation study

In this section, a simulation study is conducted to assess the performance of the proposed method. The main objective is to assess the accuracy of estimates and the correct rate of variable selection. We consider one semi-continuous variable $s_i$, two factor variables ${\omega }_{i1}$ and ${\omega }_{i2}$, and six categorical variables $y_{ij}(j=1,\dots ,6)$. We assume that $s_i$, ${\omega }_{ij}$ and $y_{ij}$ satisfy equations (1), (2) and (4) associated with (5), respectively, in which the number of fixed covariates is set at five. We generate $x_{i1}$ and $x_{i2}$ from the standard normal distribution, $x_{i3}$ and $x_{i4}$ from the binomial distribution with probability of success 0.3, and $x_5$ from the uniform distribution on $(0,1)$. All covariates were standardize to unify the scales. For ordered categorical variables, we take $c_j=c=4$, that is, each $y_{ij}$ belongs to $\{0,1,2,3,4\}$.

The true values of population parameters are set as follows: $\alpha =\gamma =0.7$, $\mathit{\beta }={(0.7,0.0,0.7,0.0,0.7,0.0,0.8)}^T$, $\gamma =0.7$, $\mathit{\psi }={(0.7,0.0,0.7,0.0,0.7,0.8,0.0)}^T$, ${\sigma }^2=1.0$, $\mathit{\mu }=0.7\times {\mathbf{1}}_6$, in which ${\mathbf{1}}_6$ is a $6\times 1$ vector with elements being unity. The factor loading matrix $\mathbf{\Lambda }$ and conviance matrix $\mathsf{\Phi }$ are taken as

$$\begin{array}{c}\hfill {\mathbf{\Lambda }}^T=\left[\begin{array}{cccccc}1.0& 0.8& 0.8& 0.0& 0.0& 0.0\\ 0.0& 0.0& 0.0& 1.0& 0.8& 0.8\end{array}\right],\mathsf{\Phi }=\left[\begin{array}{cc}1.0& 0.3\\ 0.3& 1.0\end{array}\right],\end{array}$$

(29)

in which ones and zeros in $\mathbf{\Lambda }$ are treated as fixed to identify model; the thresholds are set as ${\mathit{\delta }}_k={(-{1.5}^{\ast },0.0,1.2,{2.5}^{\ast })}^T$ for $k=1,\dots ,6$, where the elements with an asterisk are treated as fixed for model identification. Based on these setups, we generate data by first drawing latent factors from $N_2(0,\mathsf{\Phi })$, and then drawing latent responses ${\mathbf{Y}}^{\ast }$ from (5). The ordered categorical responses $\mathbf{Y}$, the indicator responses $\mathbf{U}$ and the intensity responses $\mathbf{Z}$ are sequentially generated from (5), (1) and (2). To investigate the effect of sample size on the estimates, we take $n=400$ and 1000, respectively, which represent small and large levels of sample size.

For Bayesian analysis, we consider the following inputs for hyperparameters: for the parameters involved in the measurement model, we take ${\mathit{\mu }}_0={\mathbf{0}}_6$, and ${\mathbf{\Sigma }}_0=100.0\times {\mathbf{I}}_6$; the elements in ${\mathbf{\Lambda }}_0$ corresponding to the free parameters in $\mathbf{\Lambda }$ are set at zero, and ${\mathbf{H}}_{0k}={\mathbf{I}}_2$ for $k=1,\cdots ,6$; ${\rho }_0=10.0$, and ${\mathbf{R}}_0^{-1}=6.0\times {\mathbf{I}}_3$; for the threshold parameters $\mathit{\delta }$, we take ${\eta }_1=\cdots ={\eta }_5=1.0$, which denotes the uniform distribution of p on the simplex in ${\mathbb{R}}^5$; for intercept parameters $\alpha$, $\gamma$ and ${\sigma }^2$ in the two-part model, we set ${\alpha }_0={\gamma }_0=0$, ${\sigma }_{\alpha 0}^2={\sigma }_{\gamma 0}^2=100$, and $a_0=b_0=2.0$; the hyperparameters involved in the formulation of $\mathit{\beta }$ and $\mathit{\psi }$ are set as before. Note that these values can ensure the priors inflated enough, hence it could be expected to enhance the robust of inference. In addition, we set ${\nu }_{\beta 0}={\nu }_{\psi 0}=0.001$ in equations (10) and (11) to guarantee ${\beta }_k$ and ${\psi }_k$ clumping at zero sufficiently.

The MCMC algorithm described in Section 3 is implemented to obtain the estimates of unknown parameters $\mathit{\theta }$. Before the formal implementation, a few test runs were conducted as pilots to monitor the convergence of the Gibbs sampler. We plot the values of EPSR of unknown parameters against the number of iterations under three different starting values. For SS, Figure 2 presents the plots of EPSR of unknown parameters under three different starting value with sample size $n=400$.

It can be found that the convergence of estimates is fast and all values of EPSR are less than 1.2 in about 300 iterations. To be conservative, we remove the first 2000 observations as burn-in phrase and further collect 3000 observations for calculating the bias (BIAS), the root mean squares (RMS) and the standard deviation (SD) of the estimate across 100 replications. The BIAS and RMS of the j-th component ${\widehat{\theta }}_j$ in estimates are defined as follows:

$$\begin{array}{c}\hfill BIAS({\widehat{\theta }}_j)=({\overline{\theta }}_j-{\theta }_j^0),{\overline{\theta }}_j=\frac{1}{100}\sum _{\kappa =1}^{100}{\widehat{\theta }}_j^{(\kappa )},RMS({\widehat{\theta }}_j)=\sqrt{\frac{1}{100}\sum _{\kappa =1}^{100}{({\widehat{\theta }}_j^{(\kappa )}-{\theta }_j^0)}^2,}\end{array}$$

(30)

where ${\theta }_j^0$ is the j-th element of population parameters ${\mathit{\theta }}^0$. The summaries of estimates of main parameters under two scenarios are reported in Table 1 and Table 2, where the sums of SD and RMS across the estimates are presented in the last rows.

Examinations of Table 1 and Table 2 present the following findings: (i) Both methods produce satisfactory results and the performance of SS are slightly superior to that of BaLsso. For $n=400$, the total RMS and SD are 1.870 and 1.975 respectively under SS, and amount to 2.016 and 2.035 respectively under BaLsso; (ii) As expected, increasing the sample size improves the accuracy of the estimates both for SS and BaLsso.

Another simulation is conducted to assess the performance of the proposed method in the variable selection when the covariates and latent variables are correlated. In this setting, we generate covariates and latent factors jointly from the multivariate normal distribution with mean zeros and covariance matrix $\mathbf{\Sigma }(7\times 7)$ with ${\mathsf{\Sigma }}_{jk}={\rho }^{|j-k|}$, where ${\mathsf{\Sigma }}_{jk}$ is the ${(j,k)}^{th}$ entry of $\mathbf{\Sigma }$. We consider three scenarios for $\rho$: (i) $\rho =0.1$, (ii) $\rho =0.5$ and (iii) $\rho =0.8$, which represents respectively the weak, the mild and the strong dependence among them. The values of $\mathit{\beta }$ and $\mathit{\psi }$ are taken as $(1.0,0.0,1.0,0.0,1.0,0.0,1.0)$ and $(1.0,0.0,1.0,0.0,1.0,1.0,0.0)$ respectively, and the sample size is taken as $n=1000$. The other model setups are set as the same as before. We implement MCMC sampling and collect 3000 observations after removing first 2000 observations for posterior inference. We follow [13] and treat a regression coefficient to be zero if the absolute value of its estimate is less than 0.1. Table 3 gives the summary of variable selection across 100 replications.

Based on Table 3, it can be found that (i) for nonzero regression coefficients, two methods exhibit satisfactory performances, both with 100% correct rates across all situations; (ii) for zero regression coefficients, there exist difference between two methods, and SS are uniformly outperforms BaLsso. The underlying reason perhaps is that for SS, the variances of estimates are set to be small enough to ensure the coefficients close to zero while for BaLsso, the variance of estimates are controlled by the shrinkage parameters which may not be large enough to ensure this point; (iii) with the increase of the strength of dependence, the correct rates of two methods decreases.

5. China Household Finance Survey data

To demonstrate the usefulness of the proposed methodology, in this section a small portion of Chinese household finance debt data is analyzed. The dataset is collected from the of China household financial survey (CHFS), a non-profit institute organized by the Southeast University of Finance and Economics. The survey covers a series of questions which touch on the information about various aspects of the household’s financial situation. In this study, we only focus on the measurement ‘gross debts per household (DEB)’, the amount of the secured debt and unsecured debt of a household under investigation. We extracted from the survey of Zhejiang Province in 2013. Due to some uncertain factors, some measurements in DEB are missing. The missing proportion is about 2.7%. We remove the subjects with missing entries and the ultimate sample size is 884. A preliminary data analysis shows that the measurements DEB contain excessive zeros and the proportion of zeros is about 72.58%. Naturally, we treat this variable as the outcome variable $s_i$, and identify it with $u_i$ and $z_i$. Figure 3 presents the histogram of DEB as well as the logarithms of the positive values. It can be seen clearly that data set illustrates strong heterogeneity. The skewness and kurtosis of DEB are 1.1042 and 2.3361, respectively, which indicates that single parametric model for DEB may be unappreciated.

We include the following measurements as the potential explanatory factors to interpret the variability in DEB: gender ($x_1$), age ($x_2)$, marital status ($x_3$), health condition ($x_4$), educational experience ($x_5$), employment status of the household head ($x_6$), the number of family members (aged over 16, $x_7$), and the household annual income ($x_8$). Table 4 gives the descriptive summary of the measurements under consideration. To unify the scale, all covariates were standardized.

Besides the observed factors mentioned above, we also include family culture $\eta$, a latent factor into current analysis. It is well-known that China is an ancient civilization country with a long history, and Confucian culture has deeply rooted in the social development. Economic activity or social development can not be independent of cultural development. Hence, it is of practical interest to investigate how the family culture affect the behaviour of the household finance debt. Based on the design of the questionnaire, we select the following three measurements as manifestations for $\eta$: (i) boys’s preference (BP, $y_1$). This is a three-category measurement coded by 0, 1 and 2, which corresponds to the attitude ‘oppositive’, ‘does’t matter’ and ‘strongly support’; (ii) attitude toward to the single of children (SC), coded by 0, 1 and 2, according to the leve of support; (iii) importance of a household head in a family. This measurement is originally coded in point 0 to 5 according the support level. However, in view of that the frequencies in the last three groups are small, we group them into three categories and recode them by 0 (does not matter), 1 (important) and 2 (very important). In addition, due to that some manifestations are missing, we treat missing data as missing random and ignorable [28] and ignore the specific missing mechanic that results in missing data.

Let $\mathbf{U}=\left\{u_i\right\}$, $\mathbf{Z}=\left\{z_i\right\}$, and $\mathbf{Y}=\{{\mathbf{Y}}_{obs},{\mathbf{Y}}_{mis}\}$, where ${\mathbf{Y}}_{obs}$ is the collectionn of observed data and ${\mathbf{Y}}_{mis}$ is the set of missing data. We formulate $\mathbf{U}$, $\mathbf{Z}$ and $\mathbf{Y}$ within equations (1), (2) and (5), and assume that ${\eta }_i,iid.\sim N(0,1)$. The inputs of hyperparameters in the priors are taken as follows: ${\mathsf{\Lambda }}_{j0}=0.0$, $H_{j0}=1$ and ${\eta }_{j1}={\eta }_{j2}={\eta }_{j3}=2.0$. The values of other hyperparameters are taken as the same as those in the simulation study. To implement MCMC sampling algorithm, we need to impute the missing data in $\mathbf{Y}$. This is just to do by drawing $y_{ij,mis}$ from the conditional distribution $p(y_{ij,mis}|\mathit{\theta },{\mathbf{Y}}_{obs})=N({\mu }_{j,mis}+{\mathsf{\Lambda }}_{j,mis}{\eta }_i,1)$, where ${\mu }_{j,mis}$ and ${\mathsf{\Lambda }}_{j,mis}$ are the components of $\mathit{\mu }$ and $\mathbf{\Lambda }$ respectively which corresponds to the missing entry $y_{ij,mis}$ in ${\mathbf{y}}_i$. In addition, to identify the model and scale the factor, we set ${\mathsf{\Lambda }}_1=1$. We also adopt the method in [41] in the context of latent variable model with polytomous data and fix ${\delta }_{j1}$ at ${\mathsf{\Phi }}^{-1}(f_{j1}/n_j)$, where $n_j$ is the size of $y_{obs,ij}$ equal to 1, and $f_{j1}$ is the observed frequency of 0 in $y_{obs,ij}$. To assess the convergence of the algorithm, for SS, we plot the traces of estimates under three different initial values (see Figure 4). It can be seen that the algorithm converges at about 3000 iterations. To be conservative, we collect 6000 observations after deleting the initial 4000 observations for calculating the estimates and the standard deviations.

Table 5 gives the summary of the estimates of unknown parameters in two parts and factor loadings. Examinations of Table 5 show that most estimates are very close but there exists differences in the estimates of ${\beta }_4$, ${\beta }_5$, ${\beta }_7$, ${\beta }_8$, ${\psi }_2$, ${\psi }_7$ and ${\psi }_8$. For example, the estimates of ${\beta }_4$, ${\beta }_5$ and ${\beta }_7$ under SS are $0.428$, $0.577$, $0.747$ with standard deviations $0.062$, $0.070$ and $0.072$ respectively, while equal to $0.072$, $0.082$ and $0.092$ with standard deviations $0.07$, $0.081$ and $0.092$ under BaLsso. These differences reflect the fact that two methods impose different penalties on the regression coefficients in the variable selection.

To see more clearly, Table 4 gives the resulting selected variables according to SS and BaLsso. It can be seen that (i) for part one, both methods give the same results for the selection of factors ‘gender’, ‘age’, ‘material status’, ‘employment’, ‘number of adults’ and ‘family culture’. Two methods favor that ‘age’,‘material status’, and ‘number of adults’ can be helpful in improving model fits while ‘gender’ and ‘family culture’ have less impacts on the probability of being held finance debt. However, there exist contradictory conclusion in selecting ‘health condition’, ‘education’ and ‘income’; (ii) For part two, except the factors ‘age’ and ‘number of adults’, two methods give the same results. In particular, two methods support that ‘family culture’ has a impact on the amount of household finance debts being held. This fact is also revealed by [20] in the analysis of two-part nonlinear latent variable model. The further interpretation is omitted for saving spaces.

6. Discussion

Two-part latent variable model can be considered as an extension of traditional two-part model to the situations where the latent variables are included to identify the unobserved heterogeneity of population resulting from the absence of the observed covariates. When analyzing such a model, an important issue is to determine which factor is relevant to the outcome variable. This is especially true when the number of exogenous factors is high because the usual model selection/comparison procedure is extremely time-consuming. In this paper, we restor to the Bayesian variable selection method and developed a fully Bayesian variable selection procedure for the semi-continuous data. Our formulation is along the lines with the spike and slab bimodal prior and recast the distribution of regression coefficients and factor loadings as hierarchy of priors over the parameter and model space. The selected variables is identified with high posterior probability of occurrence. We also consider a adaptive Bayesian lasso (BaLsso) for reference. To facilitate the computation, we recast the logistic regression model in part one as the flavor of normal mixture model by introducing latent Polya-gamma variables. which admits the conjugate conjugate full-conditional distributions for all regression coefficients, factor loadings and factor variables.

Although the Bayesian variable selection has its unique advantage, there are still some limitations that need to be considered with care. First, its computational complexity is high. Bayes SSL requires Monte Carlo sampling to estimate the posterior distribution, which can lead to slower calculation speed, especially when working with high-dimensional data sets. Secondly, the method is sensitive to hyperparameter and data distribution assumptions. The selection of the hyperparameters of the prior distribution, such as the ratio of spike to slab, lasso penalty parameters, and data distribution assumptions, will have a great impact on the results. When the data does not conform to the model convention, the performance of the model is poor. Therefore, these issues need to be carefully considered in practical application to ensure that the Bayesian SS method can be effectively applied to specific data sets.

The existing applications of the proposed methodology can be applied to more general latent variable models that include the the multilevel SEMs [41] and longitudinal dynamic variable models [42,49] with discrete variables. These extensions are left for further study.