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 ${\mathsf{\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 [38] 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, [9,11]) and identify $s_i$ with two surrogate variables: $u_i=I\{s_i>0\}$ and $z_i=log\left(s_i^{+}\right)$, where $I\left(A\right)$ 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: in which $\alpha$ and $\gamma$ are the intercept parameters, ${\mathsf{\beta }}_x$ and ${\mathsf{\psi }}_x$ are the vectors of regression coefficients, and ${\mathsf{\beta }}_{\omega }$ and ${\mathsf{\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 $\mathsf{\beta }={({\mathsf{\beta }}_x^T,{\mathsf{\beta }}_{\omega }^T)}^T$ and $\mathsf{\psi }={({\mathsf{\psi }}_x^T,{\mathsf{\psi }}_{\omega }^T)}^T$ and treat ${\mathbf{w}}_i={({\mathbf{x}}_i^T,{\mathsf{\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 [39,40]. This is especially true when the dimension of ${\mathsf{\omega }}_i$ is high. In this case, any auxiliary information is required to manifest ${\mathsf{\omega }}_i$ further. Among various-easy-constructs, we consider latent variable (LV, [39,40]) 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 where F is a known and fixed link function, ${\mathsf{ϵ}}_i$ is the vector of errors used to identity the idiosyncratic part of ${\mathbf{y}}_i$ that can not be explained by ${\mathsf{\omega }}_i$, and $\mathsf{\phi }$ is the vector of unknown parameters used to quantify the uncertainty of model. The information about ${\mathsf{\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: 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^{*}={(y_{i1}^{*},\cdots ,y_{ip}^{*})}^T$ is the vector of latent responses satisfying the factor analytic model: where $\mathsf{\mu }$ is a p-dimensional intercept vector, $\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 ${\mathsf{\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 ${\mathsf{\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 [41] (see squation (2) in their seminar paper), we express the logistics model (1) as the mixture model of form 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^{*}$ and augment them with $u_i$, then equation (1) can be considered as the marginal density of the joint distribution 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$\Omega ={\left\{{\mathsf{\omega }}_i\right\}}_{i=1}^n$ for the collect of factor variables, and write ${\mathbf{U}}^{*}={\left\{u_i^{*}\right\}}_{i=1}^n$, ${\mathbf{V}}^{*}={\left\{v_i^{*}\right\}}_{i=1}^n$, ${\mathbf{Y}}^{*}={\left\{{\mathbf{y}}_i^{*}\right\}}_{i=1}^n$ for the sets of latent response variables. The complete-data likelihood is given by where $I=\{i:u_i=1\}$ is the set of indices, $\mathsf{\delta }=\left\{{\delta }_{j\ell }\right\}$ is the set of threshold parametes, and $\mathsf{\theta }=\{\alpha ,\mathsf{\beta },\gamma ,\mathsf{\psi },{\sigma }^2,\mathsf{\mu },\Lambda ,\mathsf{\Phi },\mathsf{\delta }\}$ is the vector of unknown parameters. For the moment, we assume ${\theta }_j$ in $\mathsf{\theta }$ all are free.

Generally speaking, regression variables ${\mathbf{x}}_i$ and factor variables ${\mathsf{\omega }}_i$ may not have 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 the power 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 [37] and assume 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 $\mathsf{\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: 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 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 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 information 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, where `$IG(a,b)$’ denotes the inverse-gamma distribution with mean $b/(a-1)$ for $a>1$ and variance $b^2/\left({(a-1)}^2(a-2)\right)$ 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 assign values to them to ensure (13) to be dispered enough. For example, we can follow the routine in [33] 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$, 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 [37], in which the variance parameters ${\gamma }_{\beta k}^2$ and ${\gamma }_{\psi k}^2$ in equation (9) are specified via exponential distributions as follows: where ${\mathbf{\lambda }}_{\beta }^2={({\lambda }_{\beta 1}^2,\dots ,{\lambda }_{\beta q}^2)}^T$ and ${\mathbf{\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\left(\lambda \right)(\lambda >0)$ across various choices of $\lambda$. It can be seen that the larger the value of $\gamma$, the more kurtosis the density, indicating more penalties on the regression coefficient.

Due to their key role in equations (16) and (17), for ${\lambda }_{\beta }^2$ and ${\lambda }_{\psi }^2$, we assign the following gamma priors to them, i.e., 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), one can set $a_{k0}=c_{k0}=1$, $b_{k0}=d_{k0}=0.05$ to enhance the robustness of inference. This routine is followed in our empirical study.

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

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 naturally to assume that the parameters involved in the different models are independent.

Firstly, for $\mathsf{\mu }$, $\Lambda$ and $\mathsf{\Phi }$ , we consider the following conjugate priors: where `$IW(\rho ,\mathbf{R})$’ denotes the inverse Wishart distribution with degrees of freedom $\rho$ and scale matrix $\mathbf{R}$[42]; ${\Lambda }_k^T$ is the $k^{th}$ row vector of $\Lambda$; ${\mathsf{\mu }}_0$, ${\mathsf{\Sigma }}_0(p\times p)>0$, ${\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 where ${\alpha }_0$, ${\sigma }_{\alpha 0}^2$, ${\gamma }_0$, ${\sigma }_{\gamma 0}^2$ and $a_0$, $b_0$ are the fixed hyperparameters.

Lastly, for threshold parameter $\mathsf{\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\left(\mathsf{\delta }\right)={\prod }_{j=1}^{p}p\left({\mathsf{\delta }}_j\right)$, where ${\mathsf{\delta }}_j=\left({\delta }_{jk}\right)$ is the $j^{th}$ row vector of $\mathsf{\delta }$. In the following, we suppress the subscript j in ${\delta }_{jk}$ for notational simplicity and write $\mathsf{\delta }$ for ${\mathsf{\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 $\mathsf{\delta }$, we follow [43] 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 to show that the transformation is invertible with Jacobi determination unity. We first consider the following Dirichlet distribution for $p={(p_1,\cdots ,p_c)}^T$: where $\mathcal{B}({\eta }_1,\cdots ,{\eta }_{c+1})={\prod }_{j=1}^{c+1}\mathsf{\Gamma }\left({\eta }_j\right)/\mathsf{\Gamma }\left({\sum }_{j=1}^{c+1}{\eta }_j\right)$ is the multivariate beta function evaluated at ${\eta }_1,\cdots ,{\eta }_{c+1}$, and ${\eta }_j>0$. Then, by the formula of inverse transformation, the joint distribution of $\mathsf{\delta }$ is given by where $f_0\left(x\right)$ is the derivative of $F_0\left(x\right)$ with respect to x. We call (24) the transformed Dirichlet prior and use it as the prior of $\mathsf{\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

With the prior given above, the inference about $\mathsf{\theta }$ is based on the posterior $p^o\left(\mathsf{\theta }\right|\mathbf{U},\mathbf{Z},\mathbf{Y})$, which has no closed form. Motivated by the key idea in [44], 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 $\Omega$, ${\mathbf{U}}^{*}$ and ${\mathbf{Y}}^{*}$ mentioned before, we further let ${\mathbf{Q}}^{*}$ be the collection of latent quantities involved in the specifications of $\mathsf{\beta }$ and $\mathsf{\psi }$, i.e., ${\mathbf{Q}}^{*}=\{{\mathbf{F}}_{\beta }^{*}$, ${\mathbf{F}}_{\psi }^{*}$, ${\mathsf{\eta }}_{\beta }^2$, ${\mathsf{\eta }}_{\psi }^2$, $w_{\beta }$, $w_{\psi }\}$ under SS and $\{{\mathbf{\lambda }}_{\beta }^2$, ${\mathbf{\lambda }}_{\psi }^2\}$ under BaLsso. Rather than working with the posterior $p^o$ directly, we consider the following joint distribution where $p^o$ can be considered as the marginal of $p^{joint}$. We use Markov chain Monte Carlo(MCMC, [45,46]) 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:

Upon convergence, the posterior is approximated by the empirical distribution of the 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 [47] 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: where $\{{\mathsf{\theta }}^{\left(m\right)}: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{\mathsf{\beta }}$ and $\widehat{\mathsf{\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 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^{\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 probability of $f_{\beta j}=1$ and $f_{\psi j}=1$, which can be approximated by 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 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$.

In this section, we will present some technical details on the full conditionals in the MCMC sampling. For ease of exposition, for any scalar or vector x, we use $p\left(x\right|\cdots )$ to denote the conditional distribution of x given `⋯’. Note that under the scenarios SS and BaLsso, the full conditionals of $\Omega$, ${\mathbf{U}}^{*}$, ${\mathbf{Y}}^{*}$ and $\mathsf{\theta }$ are exactly the same. The following derivations are mainly based on the Bayes theorem.

1. Full conditional of $p(\Omega |\cdots )$

It follows from (8), (2) and (5) that where

Let ${\kappa }_i^{*}={\kappa }_i-u_i^{*}(\gamma +{\mathbf{x}}_i^T{\mathsf{\beta }}_x)$, $z_i^{*}=z_i-\gamma -{\mathbf{x}}_i^T{\mathsf{\psi }}_x$. By some algebra, it can be shown that where Hence, draw of $\Omega$ can be obtained by simulating ${\mathsf{\omega }}_i$ independently from the normal distribution (A1).

2. Full conditional of $p\left({\mathbf{U}}^{*}\right|\cdots )$

Following the similar derivation in [41], it can be shown that given $\mathbf{U}$, $\Omega$ and $\mathsf{\theta }$, ${\mathbf{U}}^{*}$ is the entilted Polya-Gamma distribution given by where ${\eta }_i=\alpha +{\mathsf{\beta }}^T{\mathbf{w}}_i$. Drawing $u_i$ from this distribution can be achieved via rejection sampling, see [41] or [52] for more details on this issue.

3. Full of conditional of $p\left({\mathbf{Y}}^{*}\right|\cdots )$

Note that Hence, given $\Omega$, the full conditional of ${\mathbf{Y}}^{*}$ only depends on $\mathsf{\mu }$, $\Lambda$, $\mathbf{Y}$ and $\Omega$, and is given by This is the truncated normal distribution and its draw can be obtained via inverse distribution sampling method, see for example, [53].