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

Abstract

Semi-continuous data are very common in social sciences and economics. In this paper, a Bayesian variable selection procedure is developed to assess the influence of observed and/or unobserved exogenous factors on semi-continuous data. Our formulation is based on a two-part latent variable model with polytomous responses. We consider two schemes for the penalties of regression coefficients and factor loadings: a Bayesian spike and slab bimodal prior and a Bayesian lasso prior. Within the Bayesian framework, we implement a Markov chain Monte Carlo sampling method to conduct posterior inference. To facilitate posterior sampling, we recast the logistic model from Part One as a norm-type mixture model. A Gibbs sampler is designed to draw observations from the posterior. Our empirical results show that with suitable values of hyperparameters, the spike and slab bimodal method slightly outperforms Bayesian lasso in the current analysis. Finally, a real example related to the Chinese Household Financial Survey is analyzed to illustrate application of the methodology.

1. Introduction

Semi-continuous data, which are characterized by excessive zeros, are very common in the fields of social sciences and economics. A typical example is given by [1] in an 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. For understanding such a type of data structure, a two-part model [2] is a widely appreciated statistical method. The basic assumption for a two-part model is that the overall model consists of two processes: one binary process (Part One) and one continuous positive-valued process (Part Two). The binary process, usually formulated within a logistic or probit regression model, is used to indicate whether the items have been responded to or not, while the continuous process, conditioning of the binary process, is used to describe what the actual levels of the responses are (see, e.g., [3]). By combining two processes into one, a two-part model provides a unified and flexible way to describe various relationships underlying semi-continuous data. Now, two-part models have been widely used for health services [4,5,6], medical expenditures [1,7,8,9,10], household finances [11], substance use studies [12,13], and genome analysis [14].

A traditional two-part model usually formulates exogenous explanatory factors as fixed and observed. However, in many real applications, especially for socials survey, many unobserved/latent and random factors also have important impacts on the outcome variable(s). This fact is revealed by [15] in a study of children’s aggressive behavior. Ref. [15] noted that two factors, the propensity to engage in aggressive behavior and the propensity to have highly aggressive activity levels, had significant influence on children’s aggressive behavior. The authors incorporated two such latent factors into their analysis and established a two-component–two-part mixture model to identify the heterogeneity of the population. Ref. [16] noticed that in China, the financial literacy of a family had a non-ignorable effect on the desire to hold finance debts and also affected the amount of finance debt being held. They suggested conducting a joint analysis of latent factors and observed covariates in a two-part regression model. Latent factors are further manifested by multiple binary measurements via a factor analysis model. Ref. [17] incorporated a two-part regression model into a general latent variable model framework and analyzed the internal relationships between multiple factors longitudinally. These methods have brought significant attention to two-part models in behavioral science, economics, psychology, and medicine in recent years: see, for example [14,18,19], and references therein for further developments of two-part models.

In an analysis of semi-continuous data, an important issue is to determine which explanatory factors are helpful for 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, lassos and their extensions [20,21,22,23,24,25,26,27] have been the most commonly used methods for 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/regularization approaches have been applied widely for prediction and prognosis (see, for example, [28,29]). Though more appealing, lasso-type regularization methods also suffer some limitations. For example, most contributions are developed within the frequentist framework, and their performance heavily depends on the large sample theory (see, for example [20,21,26,27], and references therein). This also readily leads to computational difficulty in the analysis of mixed data. An alternative to variable selection is conducted within the Bayesian framework. Statisticians have introduced hierarchical models with mixed spike-and-slab priors that can adaptively determine the amount of shrinkage [30,31]. The spike and slab prior is the fundamental basis for most Bayesian variable selection approaches and has proved remarkably successful [30,31,32,33,34,35,36]. Recently, Bayesian spike and slab priors have been applied to predictive modeling and variable selection in large-scale genomic studies: see [37] for a simple review. Nevertheless, model selection has never been considered in a two-part regression model with latent variables. In this study, we introduce a spike and slab model and Bayesian lasso that have been combined into a two-part latent variable model, which is a first attempt for this model.

Our formulation is more along the lines of the spike and slab bimodal prior in [34] and the Bayesian lasso in [38]. 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 a two-point mixture model with one component located at a point close to zero and the other component situated at a point far away zero. The mixing proportion is governed by a beta-distribution with suitable hyperparameters. The other scheme uses a Bayesian lasso for which the variance is specified via a gamma distribution that is scaled by the penalty parameters. The 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 chain Monte Carlo sampling method. A Gibbs sampler is used to draw observations from the posterior. We obtain all full conditionals. Posterior analysis is carried out based on the simulated observations. We investigate the performance of the proposed methods via a simulation study and a real example. Our empirical results show that the two schemes result in similar results for variable selection, but the spike and slab bimodal prior with suitable hyperparameters slightly outperforms the Bayesian lasso in terms of the correct rate.

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

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 feature extraction.

2.1. Two-Part Latent Variable Model

Suppose that for $i=1,\dots ,n$, $s_i$ is a semicontinuous outcome variable that takes a value in $[0,\infty )$; ${\mathbf{x}}_i$ is a generic vector with 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 latent/unobserved variables ${\mathit{\omega }}_i={({\omega }_{i1},\dots ,{\omega }_{im})}^T$ into the analysis to account for the 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 with the noise. Inclusion of latent variables can improve model fit and strengthen the power of model interpretation: see [39] for more discussion of the latent variables in a general setting. To deal with the spike of $s_i$ at zero, we follow the common routine in the literature (see, for example [10,12]) and identify $s_i$ with two surrogate variables: $u_i=I\{s_i>0\}$ and $z_i=log\left(s_i\right|s_i>0)$, where $I\left(A\right)$ denotes the indicator function of set A. That is, we separate the whole dataset into two parts: one part is the binary dataset that corresponds to the response-to-nonresponse indicators of the subject and the other part is the set of logarithms of positive values. Our interest focuses on the exploration of the effects of exogenous factors on the two parts.

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

$$\begin{array}{ccc}& & log\left(\frac{P(u_i=1)}{1-P(u_i=1)}\right)={\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\left(z_i\right|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 scalars of the intercept parameters, ${\mathit{\beta }}_x$ and ${\mathit{\psi }}_x$ are the vectors of the regression coefficients, and ${\mathit{\beta }}_{\omega }$ and ${\mathit{\psi }}_{\omega }$ are the vectors of the factor loadings; ${\sigma }^2$ is the scale and ‘T’ is the transpose operator of the 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. Note that Equation (1) can be represented as

$$\begin{array}{c}\hfill p\left(u_i\right|{\mathbf{x}}_i,{\mathit{\omega }}_i)=exp\left(u_i{\eta }_i^u\right)/(1+exp\left({\eta }_i^u\right)),u_i=0,1,\end{array}$$

(3)

and we refer to it as the logistic model.

The involvement of latent variables apparently complicates the model. It readily leads to model identification problems [40,41]. 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 a latent variable (LV) [40,41] approach. A basic assumption of the LV approach is that there exists, say, p manifestations ${\mathbf{y}}_i={(y_{i1},\cdots ,y_{ip})}^T$, of which each $y_{ij}$ may be continuous, counted, or categorical, and we assume 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}$$

(4)

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 cannot be explained by ${\mathit{\omega }}_i$, and $\mathit{\phi }$ is a vector of unknown parameters used to quantify the uncertainty of the model. The information about ${\mathit{\omega }}_i$ is manifested by ${\mathbf{y}}_i$ via F. In this paper, in view of real applications, we consider p ordered categorical variables ${\mathbf{y}}_i={(y_{i1},\cdots ,y_{ip})}^T$, for which $y_{ij}$ takes a 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}$$

(5)

where ${\delta }_{j,0}<{\delta }_{j,1}<\cdots <{\delta }_{j,c_j}<{\delta }_{j,c_j+1}$ are the threshold parameters satisfying ${\delta }_{j,0}=-\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 analysis model:

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

(6)

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

(7)

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

We refer to the model specified by (1), (2), and (5) associated with (6) 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 sharing 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 a traditional two-part model and ordinal regression model.

To facilitate efficient calculation and motivated by the key identity in [42] (see Equation (2) in their seminal paper), we express model (3) as a mixture model with form

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

(8)

where ${\kappa }_i=u_i-0.5$, and $p_{PG}\left(u\right)$ 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 (3) can be considered to be 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}\left(u_i^{\ast }\right).\end{array}$$

(9)

Note that the exponential part is the kernel of the 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 $\mathbf{\Omega }={\left\{{\mathit{\omega }}_i\right\}}_{i=1}^n$ for the collection of factor variables and ${\mathbf{U}}^{\ast }={\left\{u_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{Y},{\mathbf{Y}}^{\ast },\mathbf{\Omega }|\mathit{\theta })\hfill \\ \hfill =& p(\mathbf{U},{\mathbf{U}}^{\ast }|\mathbf{\Omega },\alpha ,\mathit{\beta })p\left(\mathbf{Z}\right|\mathbf{U},\mathbf{\Omega },\gamma ,\mathit{\psi },{\sigma }^2)p\left(\mathbf{Y}\right|{\mathbf{Y}}^{\ast },\mathit{\delta })p\left({\mathbf{Y}}^{\ast }\right|\mathbf{\Omega },\mathit{\mu },\mathbf{\Lambda })p(\mathbf{\Omega }|\mathbf{\Phi })\hfill \\ \hfill =& \prod _{i=1}^{n}exp\left\{{\kappa }_i{\eta }_i^u-\frac{1}{2}{u}_i^{\ast }{\left({\eta }_i^u\right)}^2\right\}{p}_{PG}\left(u_i^{\ast }\right|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}{{\left(\sqrt{2\pi }\right)}^m{|\mathbf{\Phi }|}^{1/2}}exp\left\{-\frac{1}{2}tr\left[{\mathbf{\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 parameters, and $\mathit{\theta }=\{\alpha ,\mathit{\beta },\gamma ,\mathit{\psi },{\sigma }^2,\mathit{\mu },\mathbf{\Lambda },\mathbf{\Phi },\mathit{\delta }\}$ is the vector of unknown parameters. For the moment, we assume ${\theta }_j$ in $\mathit{\theta }$ are all free.

2.2. Bayesian Feature Selection

Generally speaking, regression variables ${\mathbf{x}}_i$ and factor variables ${\mathit{\omega }}_i$ may not impact $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 the power of model interpretation. Therefore, it is necessary to determine which regression coefficient or factor loading is significantly far from zero. In the context of frequentist inference, this issue is generally tackled via stepwise regression, during which for each variable it is decided to be excluded or included according to the model fit. However, the situation becomes complex when the number of independent variables is large. In this paper, we pursue a Bayesian variable selection procedure. To this end, we follow [38] 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}$$

(10)

in which we use $diag\left\{a_k\right\}$ to represent a diagonal matrix with the $k^{th}$ diagonal element $a_k$ and let $q=r+m$. 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 the probability is 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 about ${\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}$$

(11)

$$\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}$$

(12)

where ${\delta }_a(·)$ is the Dirac measure concentrated at point a, $w_{\beta }$ is the random weight used to measure the degree of 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 the slab; ${\nu }_{\beta 0}$ is a previously specified small positive value 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\left(w_{\beta }\right|a_{\beta },b_{\beta })=Beta(a_{\beta },b_{\beta }),p\left(w_{\psi }\right|a_{\psi },b_{\psi })=Beta(a_{\psi },b_{\psi }),\end{array}$$

(13)

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

Finally, to measure the magnitudes of the ‘slab’ 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}$$

(14)

where we use ‘$IG(\tau ,\zeta )$’ to signify the inverse gamma distribution with mean $\zeta /(\tau -1)$ for $\tau >1$ and variance ${\zeta }^2/\left({(\tau -1)}^2(\tau -2)\right)$ for $\tau >2$; ${\tau }_{\beta 0}$, ${\zeta }_{\beta 0}$, ${\tau }_{\psi 0}$, and ${\zeta }_{\psi 0}$ are the hyperparameters, which are treated as fixed and known. Similarly, one can assign values to them to ensure that (14) is dispersed enough. For example, we can follow the routines in [34] for 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 (11) and (12) can be formulated as a 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}$$

(15)

$$\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}$$

(16)

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

It is instructive to compare the proposed method to the Bayesian lasso [38], in which the variance parameters ${\gamma }_{\beta k}^2$ and ${\gamma }_{\psi k}^2$ in Equation (10) 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}$$

(17)

$$\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}$$

(18)

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/penalty 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 (17) and (18) leads to marginal distributions of ${\beta }_k$ and ${\psi }_k$ as Laplace distributions with location zero and scale ${\lambda }_k$. The penalty parameters ${\lambda }_{\beta k}^2$ and ${\lambda }_{\psi k}^2$ are rather crucial for determining the amount of shrinkage of parameters. Figure 1 presents the plots of densities of Laplace distribution $LA\left(\lambda \right)(\lambda >0)$ across various choices of $\lambda$. It can be seen that the larger the value of $\lambda$ is, the more kurtosis the density has, indicating more penalties on the regression coefficient.

Due to the key roles in Equations (17) and (18) of ${\lambda }_{\beta }^2$ and ${\lambda }_{\psi }^2$, we assign the following gamma priors to them, i.e.,

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

(19)

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

(20)

where ‘$Ga(\nu ,\lambda )$’ denotes a gamma distribution with mean $\nu /\lambda$. As previously discussed, 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 [43], one can set $a_{k0}=c_{k0}=1$ and $b_{k0}=d_{k0}=0.05$ to enhance the robustness of inference. This routine is 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 the 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 the spike and slab bimodal prior to SS and the 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 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 $\mathbf{\Phi }$, we consider the following conjugate priors:

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

(21)

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

(22)

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

(23)

where ‘$IW(\rho ,\mathbf{R})$’ denotes an inverse Wishart distribution with degrees of freedom $\rho$ and scale matrix $\mathbf{R}$ [44]; ${\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 as fixed and known.

Secondly, for $\alpha$, $\gamma$, and ${\sigma }^2$ in Parts One and Two, we assume they are mutually independent and satisfy

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

(24)

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

Lastly, for 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 c. Moreover, we assume that $p\left(\mathit{\delta }\right)={\prod }_{j=1}^{p}p\left({\mathit{\delta }}_j\right)$, where ${\mathit{\delta }}_j=\left({\delta }_{jk}\right)$ 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=\Phi (·/{\tau }_0)$ for some ${\tau }_0>0$ or the distribution function of Student’s t-distribution with degrees of freedom ${\nu }_0$, where $\Phi (·)$ is the standard normal distribution function. To specify a prior for $\mathit{\delta }$, we follow [45] and let $p_j=F_0\left({\delta }_j\right)-F_0\left({\delta }_{j-1}\right)$ for $j=1,\cdots ,c$. It is easily shown that this transformation is invertible with Jacobi determination unity. We first consider the following Dirichlet distribution for $p={(p_1,\cdots ,p_c)}^T$:

$$\begin{array}{c}\hfill \pi \left(p\right)=\frac{1}{\mathcal{B}({\eta }_1,\dots ,{\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,\dots ,{\eta }_{c+1})={\prod }_{j=1}^{c+1}\Gamma \left({\eta }_j\right)/\Gamma \left({\sum }_{j=1}^{c+1}{\eta }_j\right)$ is the multivariate beta function evaluated at ${\eta }_1,\dots ,{\eta }_{c+1}$, and ${\eta }_j>0$. Then, by the formula of inverse transformation, the joint distribution of $\mathit{\delta }$ is given by

$$\begin{array}{c}\hfill \pi \left(\mathit{\delta }\right)=\frac{1}{\mathcal{B}({\eta }_1,\dots ,{\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\left({\delta }_j\right)I\{{\delta }_1<\cdots <{\delta }_c\},\end{array}$$

(25)

where $f_0\left(x\right)$ is the derivative of $F_0\left(x\right)$ with respect to x. We call (25) the transformed Dirichlet prior and use it as the prior of $\mathit{\delta }$. An advantage of working with (25) 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\left({\delta }_j\right)-F_0\left({\delta }_{j-1}\right)}{F_0\left({\delta }_{j+1}\right)-F_0\left({\delta }_{j-1}\right)}|({\delta }_{j-1},{\delta }_{j+1})\sim Beta({\eta }_j,{\eta }_{j+1}),(j=1,\dots ,c).$$

(26)

3.2. MCMC Sampling

With the prior given above, the inference about $\mathit{\theta }$ is based on the posterior $p^o\left(\mathit{\theta }\right|\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 augment them to the observed data to form the complete data. The statistical inference is carried out based on the complete data likelihood. To this end, apart from $\mathbf{\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 }\}$ under SS and $\{{\mathit{\lambda }}_{\beta }^2$, ${\mathit{\lambda }}_{\psi }^2\}$ under BaLsso. Rather than working with the posterior $p^o$ directly, we consider the following joint distribution

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

(27)

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

• Draw $\mathbf{\Omega }$ from $p(\mathbf{\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 \mathbf{\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 \mathbf{\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 \mathbf{\Omega },{\mathbf{U}}^{\ast },{\mathbf{Y}}^{\ast },\mathit{\theta },\mathbf{U},\mathbf{Z},\mathbf{V})$;
• Draw $\mathit{\theta }$ from $p(\mathit{\theta }\mid \mathbf{\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 the simulated observations. The convergence of the algorithm can be monitored by plotting the traces of estimates for different starting values or observing the values of EPSR [49] of unknown parameters. The technical details for implementing MCMC sampling are given in Appendix A.

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 }}^{\left(m\right)},\end{array}$$

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

The main purpose of introducing SS and BaLsso is to screen the variables in ${\mathbf{w}}_i$. Unlike that in the frequentist inference, Bayesian variable selection does not produce estimates $\widehat{\mathit{\beta }}$ and $\widehat{\mathit{\psi }}$ that are 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 (PCIs) of ${\beta }_j$ and ${\psi }_j$, which are given by

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

(28)

where $\alpha$ is any previously specified value in $(0,1)$. Calculation of the PCI can be achieved via the Monte Carlo method. For example, let $\{{\beta }_j^{\left(k\right)}: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 probabilities 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}^{\left(k\right)}=1\},{\widehat{f}}_{\psi j}=\frac{1}{K}\sum _{k=1}^{K}I\{f_{\psi j}^{\left(k\right)}=1\},\end{array}$$

(29)

where $f_{\beta j}^{\left(k\right)}$ and $f_{\psi j}^{\left(k\right)}(k=1,\dots ,K)$ are the k observations drawn from the posterior distribution via the Gibbs sampler. The variable $w_j$ is selected in Parts 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 (5) associated with (6), respectively, in which the number of fixed covariates is set to five. We generate $x_{i1}$ and $x_{i2}$ from the standard normal distribution, $x_{i3}$ and $x_{i4}$ from a binomial distribution with a probability of success of 0.3, and $x_5$ from a uniform distribution on $(0,1)$. All covariates were standardized 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 the 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$, where ${\mathbf{1}}_6$ is a $6\times 1$ vector for which the elements are in unity. The factor loading matrix $\mathbf{\Lambda }$ and covariance matrix $\mathbf{\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],\mathbf{\Phi }=\left[\begin{array}{cc}1.0& 0.3\\ 0.3& 1.0\end{array}\right],\end{array}$$

(30)

in which the ones and zeros in $\mathbf{\Lambda }$ are treated as fixed to identify the 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,\mathbf{\Phi })$ and then drawing latent responses ${\mathbf{Y}}^{\ast }$ from (6). The the indicator responses $\mathbf{U}$, the intensity responses $\mathbf{Z}$ and the ordered categorical responses $\mathbf{Y}$, are sequentially generated from (1), (2) and (5). 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 the 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,\dots ,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 the intercept parameters $\alpha$, $\gamma$, and scale ${\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 that the corresponding priors are inflated enough; hence, we can expect this to enhance the robustness of the inference. In addition, we set ${\nu }_{\beta 0}={\nu }_{\psi 0}=0.001$ in Equations (11) and (12) 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 are conducted as a pilot 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 values 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 the burn-in phase and further collect 3000 observations for calculating the bias (BIAS), root mean square (RMS), and standard deviation (SD) of the estimate across 100 replications. The BIAS and RMS of the j-th component ${\widehat{\theta }}_j$ in $\widehat{\mathit{\theta }}$ are respectively defined as follows:

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

(31)

where ${\theta }_j^0$ is the j-th element of the population parameters ${\mathit{\theta }}^0$. The summaries of the estimates of the main parameters for the two scenarios are reported in

Table 1. Summary of the estimates of unknown parameters under SS and BaLsso: simulation study and $n=400$.

	SS			BaLsso
PAR	BIAS	RMS	SD	BIAS	RMS	SD
$\alpha =0.7$	−0.015	0.097	0.129	0.028	0.150	0.134
${\beta }_1=0.7$	−0.056	0.143	0.142	−0.152	0.217	0.136
${\beta }_2=0.0$	−0.001	0.021	0.061	−0.019	0.042	0.079
${\beta }_3=0.7$	−0.144	0.216	0.145	−0.122	0.251	0.148
${\beta }_4=0.0$	0.005	0.030	0.064	−0.008	0.040	0.078
${\beta }_5=0.7$	−0.091	0.147	0.137	−0.045	0.135	0.137
${\beta }_6=0.0$	0.017	0.028	0.075	0.026	0.055	0.096
${\beta }_7=0.8$	−0.187	0.237	0.184	−0.126	0.209	0.184
$\gamma =0.7$	0.010	0.079	0.084	0.008	0.063	0.085
${\psi }_1=0.7$	−0.035	0.079	0.077	−0.011	0.065	0.074
${\psi }_2=0.0$	0.005	0.032	0.051	−0.018	0.031	0.054
${\psi }_3=0.7$	−0.007	0.061	0.070	−0.021	0.085	0.069
${\psi }_4=0.0$	−0.007	0.029	0.049	−0.003	0.031	0.053
${\psi }_5=0.7$	−0.070	0.093	0.077	−0.018	0.082	0.075
${\psi }_6=0.8$	−0.040	0.086	0.089	−0.020	0.069	0.088
${\psi }_7=0.0$	−0.011	0.033	0.062	0.014	0.036	0.069
${\sigma }^2=1.0$	0.085	0.129	0.117	0.038	0.082	0.111
${\lambda }_{21}=0.8$	0.042	0.078	0.073	0.058	0.098	0.071
${\lambda }_{31}=0.8$	0.030	0.072	0.071	0.034	0.063	0.072
${\lambda }_{52}=0.8$	0.058	0.079	0.072	0.052	0.090	0.073
${\lambda }_{62}=0.8$	0.031	0.060	0.072	0.037	0.064	0.073
${\Phi }_{12}=0.3$	0.014	0.041	0.074	0.018	0.058	0.076
Total	-	1.870	1.975	-	2.016	2.035

and

Table 2. Summary of the estimates of unknown parameters under SS and BaLsso: simulation study and $n=1000$.

	SS			BaLsso
PAR	BIAS	RMS	SD	BIAS	RMS	SD
$\alpha =0.7$	0.052	0.096	0.087	0.009	0.092	0.087
${\beta }_1=0.7$	0.005	0.069	0.089	0.055	0.117	0.090
${\beta }_2=0.0$	0.003	0.048	0.058	0.032	0.052	0.060
${\beta }_3=0.7$	0.007	0.086	0.093	−0.045	0.076	0.091
${\beta }_4=0.0$	0.004	0.015	0.049	−0.020	0.043	0.060
${\beta }_5=0.7$	0.010	0.071	0.086	0.013	0.074	0.085
${\beta }_6=0.0$	−0.003	0.029	0.059	0.032	0.064	0.077
${\beta }_7=0.8$	0.002	0.102	0.120	−0.042	0.108	0.114
$\gamma =0.7$	0.017	0.042	0.053	0.030	0.056	0.054
${\psi }_1=0.7$	−0.023	0.038	0.046	−0.016	0.039	0.047
${\psi }_2=0.0$	−0.007	0.019	0.033	−0.005	0.018	0.037
${\psi }_3=0.7$	−0.028	0.060	0.042	−0.014	0.026	0.043
${\psi }_4=0.0$	−0.007	0.023	0.033	0.000	0.018	0.036
${\psi }_5=0.7$	−0.005	0.035	0.046	0.003	0.043	0.047
${\psi }_6=0.8$	−0.031	0.058	0.053	−0.039	0.063	0.054
${\psi }_7=0.0$	−0.001	0.031	0.045	−0.025	0.081	0.053
${\sigma }^2=1.0$	0.018	0.049	0.068	0.041	0.053	0.071
${\lambda }_{21}=0.8$	0.021	0.041	0.045	0.033	0.038	0.045
${\lambda }_{31}=0.8$	0.016	0.049	0.045	0.028	0.038	0.045
${\lambda }_{52}=0.8$	0.032	0.049	0.045	0.054	0.057	0.045
${\lambda }_{62}=0.8$	0.043	0.059	0.046	0.043	0.054	0.046
${\Phi }_{12}$	0.016	0.043	0.049	0.005	0.037	0.048
Total	-	1.112	1.290	-	1.247	1.335

, where the sums of the SDs and RMSs across the estimates are presented in the last rows.

Examination of Table 1 and Table 2 gives the following findings: (i) Both methods produce satisfactory results. The performance of SS is slightly superior to that of BaLsso in terms of the sums of RMS and SD. 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) Except ${\psi }_2$ in Table 1 and ${\psi }_4$ in Table 2, for the regression coefficients and factor loadings with true values at zero, the RMSs and SDs produced by SS are uniformly smaller than those produced by BaLsso. This indicates that SS imposes more penalties on the regression coefficients and factor loadings than BaLsso to shrink them toward zero. (iii) For non-zero ${\beta }_j$ and ${\psi }_j$, the situation becomes fiendishly complex. For example, the RMSs and SDs of the estimates of ${\beta }_5$, ${\beta }_7$, ${\psi }_1$, ${\psi }_5$, and ${\psi }_6$ in Table 1 show that the SS does a weaker job than BaLsso. This is also true for ${\beta }_3$ and ${\psi }_3$ in Table 2. The underlying reason perhaps lies in that the penalties imposed on the non-zero regression coefficients and factor loadings by SS are similar to those of BaLsso, and one is not overwhelmingly superior to the other. (iv) As expected, increasing the sample size improves the accuracy of the estimates for both SS and BaLsso.

Another simulation is conducted to assess the performance of the proposed method for 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 zero and covariance matrix $\mathbf{\Sigma }(7\times 7)$ with ${\Sigma }_{jk}={\rho }^{|j-k|}$, where ${\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 represent, respectively, weak, mild, and 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)}^T$ and ${(1.0,0.0,1.0,0.0,1.0,1.0,0.0)}^T$, respectively, and the sample size is taken as $n=1000$. The remainder of the model is set up the same as before. We implement MCMC sampling and collect 3000 observations after removing the first 2000 observations for posterior inference. We follow [43] and treat a regression coefficient as zero if the absolute value of its estimate is less than 0.1.

Table 3. Number of correctly selected variables in the two-part model on the simulated datasets.

	SS			BaLsso
PAR	$\mathit{\rho }\mathbf{=}\mathbf{0.1}$	$\mathit{\rho }\mathbf{=}\mathbf{0.5}$	$\mathit{\rho }\mathbf{=}\mathbf{0.8}$	$\mathit{\rho }\mathbf{=}\mathbf{0.1}$	$\mathit{\rho }\mathbf{=}\mathbf{0.5}$	$\mathit{\rho }\mathbf{=}\mathbf{0.8}$
${\beta }_1=1.0$	100	100	100	100	100	100
${\beta }_2=0.0$	98	96	85	88	86	76
${\beta }_3=1.0$	100	100	100	100	100	100
${\beta }_4=0.0$	96	95	86	93	93	85
${\beta }_5=1.0$	100	100	100	100	100	100
${\beta }_6=0.0$	96	94	93	97	92	87
${\beta }_7=1.0$	100	100	100	100	100	100
${\psi }_1=1.0$	99	100	100	100	100	100
${\psi }_2=0.0$	100	99	95	100	98	93
${\psi }_3=1.0$	100	100	100	100	100	100
${\psi }_4=0.0$	100	100	97	98	100	91
${\psi }_5=1.0$	100	100	100	100	100	100
${\psi }_6=1.0$	100	100	100	100	100	100
${\psi }_7=0.0$	100	98	97	97	96	96

gives the summary of variable selection across 100 replications.

Based on Table 3, it can be found that (i) for nonzero regression coefficients, the two methods exhibit satisfactory performances, with both having 100% correct rates across all situations; (ii) for zero regression coefficients, there are differences between the two methods, and SS uniformly outperforms BaLsso. The underlying reason perhaps is that for SS, the variances of estimates are set to be small enough to ensure that the coefficients are close to zero, while for BaLsso, the variances of estimates are controlled by the shrinkage parameters, which may not be large enough to ensure this point; (iii) with an increase in the strength of dependence, the correct rates of the two methods decrease. As an illustration, Figure 3 gives the plots of the correct rates of the selected variables for three scenarios; the values on the x-axis are the true values of ${\beta }_1$ to ${\psi }_7$ in the order given in Table 3.

5. Chinese Household Financial 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 Chinese Household Financial Survey (CHFS), a non-profit institute organized by the Southeast University of Finance and Economics in Chengdu, Sichuan Province, China. The survey covers a series of questions that touch on information about various aspects of a household’s financial situation. In this study, we only focus on the measurement ‘gross debts per household (DEB)’: the amount of secured debt and unsecured debt for the household under investigation. We extracted the data from a survey of Zhejiang Province in 2013. Due to some uncertain factors, some measurements of 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 DEB measurements 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 4 presents the histogram of DEB as well as the logarithms of positive values. It can be seen clearly that the dataset illustrates strong heterogeneity. The skewness and kurtosis of DEB are 1.1042 and 2.3361, respectively, which indicates that a single parametric model for DEB may be unappreciated.

We include the following measurements as the potential explanatory factors to interpret the variability of 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. Descriptive statistics of explanatory variables: CHFS data.

Variable	Description	Mean	Max	Min	SD
Gender ($x_1$)	=1, male; =0, otherwise	0.756	1	0	0.430
Age ($x_2$)		51.81	91	19	14.931
Marital status ($x_3$)	=1, married; 0, otherwise	0.863	1	0	0.344
Health condition ($x_4)$	=1, good; 0, otherwise	0.833	1	0	0.373
Educational experience $\left(x_5\right)$	=1, high school or above;
	=0, otherwise	0.352	1	0	0.478
Employment ($x_6$)	=1, yes; 0, otherwise	0.092	1	0	0.290
No. of adults ($x_7$)		3.002	3	0	1.301
Annual Income (CYN) $\left(x_8\right)$ *		${9.376}^4$	${8.060}^5$	0	${4.249}^4$

* Note: Superscripts are used to indicate values raised to the power of 10 (thus, $a^b=a\times {10}^b$). The measurement is taken as the middle value of the range in the questionnaire.

gives the descriptive summary of these measurements. To unify the scale, all covariates were standardized.

Besides the observed factors mentioned above, we also include family culture $\eta$ as a latent factor into the current analysis. It is well-known that China is an ancient civilization with a long history, and Confucian culture is deeply rooted in social development. Economic activity or social development cannot be independent of cultural development. Hence, it is of practical interest to investigate how family culture affects behavior related to household finance debt. Based on the design of the questionnaire, we select the following three measurements as manifestations for $\eta$: (i) The preference for boys (BP, $y_1$): this is a three-category measurement coded as 0, 1, and 2, which correspond to the attitudes ‘opposed’, ‘doesn’t matter’, and ‘strongly support’; (ii) The attitude toward having a single child (SC), coded by 0, 1 and 2, according to the level of support. (iii) The importance of family in one’s life: this measurement was originally based on a six-point scale (0 to 5) according to the support level. However, in view of the frequencies in the last three groups being small, we grouped them into three categories and recoded them as 0 (does not matter), 1 (important), and 2 (very important). In addition, as some manifestations are missing, we treat missing data as missing randomly and ignorable [50] and ignore the specific mechanics that result 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 collection 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 (6) and assume that ${\eta }_i,iid.\sim N(0,1)$. The inputs of the hyperparameters in the priors are taken as follows: ${\Lambda }_{j0}=0.0$, $H_{j0}=1$, and ${\eta }_{j1}={\eta }_{j2}={\eta }_{j3}=2.0$. The values of other hyperparameters are taken to be the same as those in the simulation study. To implement the MCMC sampling algorithm, we need to impute the missing data in $\mathbf{Y}$. This is done by drawing $y_{ij,mis}$ from the conditional distribution $p\left(y_{ij,mis}\right|\mathit{\theta },{\mathbf{Y}}_{obs})=N({\mu }_{j,mis}+{\Lambda }_{j,mis}{\eta }_i,1)$, where ${\mu }_{j,mis}$ and ${\Lambda }_{j,mis}$ are the components of $\mathit{\mu }$ and $\mathbf{\Lambda }$, respectively, that correspond to the missing entries $y_{ij,mis}$ in ${\mathbf{y}}_i$. In addition, to identify the model and scale the factor, we set ${\Lambda }_1=1$. We also adopt the method in [51] in the context of a latent variable model with polytomous data and fix ${\delta }_{j1}$ at ${\Phi }^{-1}(f_{j1}/n_j)$, where $n_j$ is the size of $y_{obs,ij}$ that is 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 5). 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 their standard deviations.

Table 5. Estimates and standard deviation estimates of unknown parameters under SS and BaLsso: CHFS data.

	SS		BaLsso		SS		BaLsso
Par	Est.	SD	Est.	SD	Par	Est.	SD	Est.	SD
$\alpha$	−0.835	0.078	−0.838	0.080	$\gamma$	9.782	0.152	9.670	0.125
${\beta }_1$	0.050	0.063	0.076	0.070	${\psi }_1$	−0.137	0.103	−0.107	0.088
${\beta }_2$	−0.750	0.099	−0.757	0.102	${\psi }_2$	−0.147	0.141	−0.015	0.081
${\beta }_3$	0.107	0.085	0.147	0.088	${\psi }_3$	−0.022	0.065	−0.006	0.075
${\beta }_4$	0.428	0.062	0.072	0.070	${\psi }_4$	−0.019	0.060	−0.029	0.069
${\beta }_5$	0.577	0.070	0.082	0.081	${\psi }_5$	0.259	0.123	0.322	0.107
${\beta }_6$	0.004	0.040	0.005	0.052	${\psi }_6$	0.035	0.058	0.053	0.067
${\beta }_7$	0.118	0.079	0.130	0.079	${\psi }_7$	0.043	0.072	0.281	0.113
${\beta }_8$	0.747	0.073	0.092	0.077	${\psi }_8$	0.384	0.132	0.188	0.118
${\beta }_{\eta }$	−0.059	0.112	−0.039	0.092	${\psi }_{\eta }$	1.205	0.106	1.910	0.104
${\sigma }^2$	0.312	0.150	0.300	0.152
${\lambda }_{21}$	−0.791	0.062	−0.714	0.057
${\lambda }_{31}$	−0.865	0.067	−0.625	0.068

gives the summary of two estimates of unknown parameters in the two-part model and factor analysis model. Examination of Table 5 shows that most of the estimates for both models are very close, but there exist 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$, and $0.747$ with standard deviations of $0.062$, $0.070$, and $0.072$, respectively, while they equal $0.072$, $0.082$, and $0.092$ with standard deviations of $0.07$, $0.081$, and $0.092$ under BaLsso. These differences reflect the fact that the two methods impose different penalties on the regression coefficients during variable selection.

To see more clearly,

Table 6. The selected variables in the CHFS data—0: excluded and 1: included.

	Part One		Part Two
VAR	SS	BaLsso	SS	BaLsso
Gender	0	0	1	1
Age	1	1	1	0
Marital status	1	1	0	0
Health condition	1	0	0	0
Education	1	0	1	1
Employment	0	0	0	0
No. of adults	1	1	0	1
Income	1	0	1	1
Family culture	0	0	1	1

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’, ‘marital status’, ‘employment’, ‘number of adults’, and ‘family culture’. The two methods favor ‘age’,‘marital status’, and ‘number of adults’ as helpful for improving model fit, while ‘gender’ and ‘family culture’ have less influence on the probability of holding household finance debt. However, there exist contradictory conclusions when selecting ‘health condition’, ‘education’, and ‘income’. (ii) For Part Two, except for the factors ‘age’ and ‘number of adults’, the two methods give the same results. In particular, both methods support that ‘family culture’ is relevant to the amount of household finance debt being held. This fact is also revealed by [18] in an analysis of the CHFS using a two-part nonlinear latent variable model. Further interpretation is omitted in order to save space.

6. Discussion

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

Although Bayesian variable selection has its unique advantages, there are still some limitations that need to be considered with care. Firstly, 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 datasets. 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 do not conform to the model’s convention, the performance of the method is poor. Therefore, these issues need to be carefully considered in real applications to ensure that the Bayesian SS method can be effectively applied to the specific dataset.

The proposed method can be applied to more general latent variable models, including multilevel SEMs [51] and longitudinal dynamic variable models with discrete variables [17,52]. These extensions are left for further study.