Variational Bayesian Inference for Quantile Regression Models with Nonignorable Missing Data

Abstract

Quantile regression models are remarkable structures for conducting regression analyses when the data are subject to missingness. Missing values occur because of various factors like missing completely at random, missing at random, or missing not at random. All these may result from system malfunction during data collection or human error during data preprocessing. Nevertheless, it is important to deal with missing values before analyzing data since ignoring or omitting missing values may result in biased or misinformed analysis. This paper studies quantile regressions from a Bayesian perspective. By proposing a hierarchical model framework, we develop an alternative approach based on deterministic variational Bayes approximations. Logistic and probit models are adopted to specify propensity scores for missing manifests and covariates, respectively. Bayesian variable selection method is proposed to recognize significant covariates. Several simulation studies and real examples illustrate the advantages of the proposed methodology and offer some possible future research directions.

1. Introduction

Koenker and Bassett [1] introduced the concept of the quantile regression (QR) model. Unlike the linear regression (OLS) model, OLS focuses solely on modeling the mean of the response variable. QR allows us to estimate regression coefficients at different quartiles of the response variable. This enables us to capture the varying degrees of correlation between the covariates and the response variable. Quantile regression models are widely used in various fields such as finance, biology, and ecology (Baur and Dirk [2]; Huang et al. [3]; Cade and Noon [4]) due to their robust to outliers in the data, robustness, and wide applicability. Yu and Moye [5] examined Bayesian quantile regression (QR) by redefining QR as an asymmetric Laplace distribution. In continuation of this work, Kozumi and Kobayashi [6] proposed an effective Gibbs sampling method for the QR model by leveraging the decomposition property of the asymmetric Laplace distribution. Alhamzawi and Rahim [7] explored the composite quantile regression model with Bayesian estimation, while Yuan et al. [8] focused on studying Bayesian composite quantile regression for a single indicator model. Additionally, Hu et al. [9] introduced a Bayesian joint quantile regression model. An important aspect of building quantile regression models is the selection of predictor variables. Li et al. [10] approached the regularization problem of quantile regression from a Bayesian perspective, focusing on variable selection. To address the instability of posterior estimates caused by variable selection in Gibbs sampling and convergence issues due to fuzzy priors, Alhamzawi and Yu [11] proposed random search variable selection (ISSVS) within a Bayesian framework. Mallick et al. [12] introduced reciprocal Lasso (rLasso) regularization, which offers advantages over Lasso regularization in terms of estimation, prediction, and variable selection within the Bayesian framework. Considerable progress has been made in Bayesian frameworks for both statistical inference and variable selection in quantile regression. However, it is important to note that all the aforementioned studies are based on fully observable data.

Missing data are inevitable in research and can result from various uncontrollable factors. For instance, some results may be lost due to machine malfunction, while in questionnaires, individuals may be unwilling to disclose their income. The study of missing data originated in the 1970s, and Little and Rubin [13] defined three missing data mechanisms based on the causes of missingness: missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). The first two mechanisms, MCAR and MAR, are unrelated to the missing data themselves. However, the third mechanism, MNAR, is directly associated with the missing data and closely reflects real life. Ignoring non-negligible missing data, especially those reflecting practical scenarios through the MNAR mechanism, can lead to erroneous conclusions. Recognizing the distinctive attributes of such data and the multitude of benefits offered by quantile regression, several researchers have explored quantile regression models as a means to address non-negligible missing data. For instance, Yuan and Yin [14] examined Bayesian quantile regression models for longitudinal data with non-negligible missingness. Zhao et al. [15] conducted a study on several estimation methods using inverse probability weighting (IPW) to estimate parameters in quantile regression models when either covariates or response variables more with non-negligible missing values. Wang and Tang [16] investigated nonlinear dynamic factor analysis models with mixed discrete and non-negligible missing covariates, employing Bayesian quantile regression models for statistical inference. Tang et al. [17] and Mulati and Tang [18] explored nonlinear dynamic factor analysis models with non-negligible missing data using Bayesian non-parametric and semi-parametric approaches, respectively. In all of the aforementioned studies, the Bayesian framework was adopted, and Bayesian inference was facilitated using the Markov Chain Monte Carlo (MCMC) algorithm. However, the MCMC algorithm has certain inherent limitations. It can be challenging to sample from, especially when dealing with large datasets or complex models.

The variational Bayesian algorithm transforms the high-dimensional integration problem in Bayesian inference into an optimization problem, enabling efficient computation of the evidence lower bound and obtaining a variational approximation of the posterior distribution, rendering it well-suited for analyzing extensive datasets. The development of variational inference (VI) has also yielded advantages for high-dimensional data, such as stochastic variational inference (SVI) by Hoffman et al. [19], “black box” variational inference by Ranganath et al. [20], and amortized variational inference by Ganguly et al. [21]. Various studies have utilized VI in different contexts, encompassing Bayesian approximation inference for models with missing covariates, as demonstrated in the work of Faes et al. [22], exploring parameter uncertainty in dynamic factorial models with missing data as shown by Erik [23], and conducting statistical inference on identified gene regulatory networks with missing data through variational Bayesian inference as presented in the research by Liu et al. [24]. However, none of the existing works have applied the variational Bayesian algorithm to infer quantile regression (QR) models with missing data.

The main contributions of this paper include the following: (i) a variational Bayesian approach for parameter estimation in quantile regression models; (ii) a variational Bayesian approach for quantile regression models with missing covariates and response variables, in which the benefits of our suggested method are demonstrated through simulations; and (iii) a variational Bayesian approach for variable selection in quantile regression models with missing data, which does the job very well.

This paper is organized as follows: Section 2 presents the variational Bayesian algorithm, QR, lasso penalty term, and missing data. In Section 3, variational inference is performed for QR in the presence of missing covariates and missing response variables, respectively. In Section 4, the variable selection and parameter estimation are performed for QR with lasso penalty term in the presence of missing covariates and response variables, respectively. The feasibility of variational Bayes in QR with missing data and the superiority of comparing Gibbs’ algorithm are illustrated by four experiments in Section 5. Section 6 analyzes the example data using the variational Bayes algorithm and obtains better results.

2. Model and Notation

2.1. Quantile Regression Model

The quantile regression model was initially introduced by Koenker and Bassett [1], and the model can be represented as follows:

(1)$y_i=X_i^T{\beta }_p+{\epsilon }_i,(i=1,\dots ,n)$

$X_i=\left(x_{i1},\dots ,x_{ik}\right)$ represents a $k\times 1$ vector of covariates, and $y_i$ denotes a response variable. Let p ($0<p<1$) represent the quantile level. The parameter vector ${\beta }_p$ is a k-dimensional vector that is associated with the quantile level p. The random error term is denoted as ${\epsilon }_i$. In this case, we consider a more general distribution for the error term that satisfies $P({\epsilon }_i<0|x_i)=p$.

For the quantile regression model, given a specific $X_i$ and $y_i$, the model can be expressed as $Q_p\left(y_i\right|X_i)=X_i^T{\beta }_p$. The estimation of ${\beta }_p$ can be obtained by minimizing the following formula:

$\sum _1^m{\rho }_p\left(y_i-X_i^T{\beta }_p\right).$

The loss function is defined as ${\rho }_p\left(u\right)=u\left(p-I_{(u<0)}\right)$, where $I(·)$ represents the indicator function. Due to the non-differentiability of the loss function, obtaining estimates for ${\beta }_p$ becomes challenging. Yu and Moye [5] demonstrated that the quantile regression model can be estimated within a Bayesian framework when the error term follows an asymmetric Laplacian distribution (ALD). However, ALD is not a standard distribution, making its posterior density function complex and inevitably increasing the computational burden. According to Kozumi and Kobayashi [6], ALD can be expressed as a mixture of exponential and normal distributions. Specifically,

(2)${\epsilon }_{\mathrm{i}}=\phi {\mathrm{z}}_{\mathrm{i}}+\theta \sqrt{\sigma {\mathrm{z}}_{\mathrm{i}}}{\mathrm{v}}_{\mathrm{i}}.$

Among them, $z_i\sim exp\left(1/\sigma \right)$ and $v_i\sim N(0,1)$ are independent of each other. Here, $\phi =\frac{1-2p}{p(1-p)}$ and $\theta =\sqrt{\frac{2}{p(1-p)}}$. Finally, the quantile regression model can be expressed as a layered model, represented as follows:

(3)$\left\{\begin{array}{c}{y}_i=X_i^T\beta +\phi z_i+\theta \sqrt{\sigma z_i}{v}_i\hfill \\ z_i\mid \sigma \sim exp(1/\sigma )=\sigma exp\left(-\sigma z_i\right)\hfill \\ v_i\sim N(0,1)=1/\sqrt{2\pi }exp\left(-v_i^2/2\right)\hfill \\ \phi =\frac{1-2p}{p(1-p)},\theta =\sqrt{\frac{2}{p(1-p)}}.\hfill \end{array}\right.$

2.2. Elements of Variational Bayes

Taking into account the latent variable $\mathbb{Z}={\mathbb{Z}}_1,\dots ,{\mathbb{Z}}_m$ and the observed values $X=x_1,\dots ,x_n$, the joint density function can be expressed as follows:

$P(\mathbb{Z},X)=P\left(\mathbb{Z}\right)P(X\mid \mathbb{Z}).$

Variational Bayesian inference no longer relies on sampling to approximate the posterior density of complex problems. Instead, it employs optimization techniques. Within an acceptable range of error, a simpler density function $q\left(\mathbb{Z}\right)\in \mathsf{\Theta }$ (where $\mathsf{\Theta }$ represents the space of possible density functions for the latent variable $\mathbb{Z}$) is used to approximate the posterior density $P(\mathbb{Z}\mid X)$. The difference between these two density functions is measured using the Kullback–Leibler divergence.

(4)$D_{KL}\left(q\right|p)=\sum _{\mathbb{Z}}q\left(\mathbb{Z}\right)log\frac{q\left(\mathbb{Z}\right)}{p(\mathbb{Z}\mid X)}=\sum _{\mathbb{Z}}q\left(\mathbb{Z}\right)log\frac{q\left(\mathbb{Z}\right)}{p(\mathbb{Z},X)}+logp\left(\mathbb{Z}\right).$

Alternatively,

(5)$logp\left(\mathbb{Z}\right)=D_{KL}\left(q\right|p)-\sum _{\mathbb{Z}}q\left(\mathbb{Z}\right)log\frac{q\left(\mathbb{Z}\right)}{p(\mathbb{Z}\mid X)}=D_{KL}\left(q\right|p)+L\left(q\right).$

The term $L\left(q\right)$ is referred to as the evidence lower bound. To minimize the KL divergence, we maximize the lower bound.

(6)$L\left(q\right)=\sum _{\mathbb{Z}}q\left(\mathbb{Z}\right)logp(\mathbb{Z},X)-\sum _{\mathbb{Z}}q\left(\mathbb{Z}\right)logq\left(\mathbb{Z}\right)={\mathbb{E}}_q\left[log\frac{p(\mathbb{Z},X)}{q\left(\mathbb{Z}\right)}\right].$

This is equivalent to solving the following optimization problem:

(7)$q^{*}\left(\mathbb{Z}\right)=argmin{D}_{KL}\left(q\right|p)$

By solving the optimization problem, the best approximation function $q^{*}\left({\mathbb{Z}}_i\right)$ can be obtained as follows:

(8)$q^{*}\left(\mathbb{Z}i\right)\propto exp{\mathbb{E}}_{-{\mathbb{Z}}_i}logp(\mathbb{Z},X)(i=1,\dots ,m)$

or

(9)$q^{*}\left(\mathbb{Z}i\right)\propto exp{\mathbb{E}}_{-{\mathbb{Z}}_i}logp({\mathbb{Z}}_i\mid {\mathsf{\Theta }}_r)(i=1,\dots ,m)$

where ${\mathsf{\Theta }}_r=X,{\mathbb{Z}}_1,\dots ,{\mathbb{Z}}_{i-1},{\mathbb{Z}}_{i+1},\dots ,{\mathbb{Z}}_m$. For convenience, in the following paragraphs, $q^{*}(·)$ is denoted as $q(·)$.

Additionally, the complexity of the density of the latent variables determines the complexity of the optimization algorithm. Therefore, we consider the mean-field variational family:

$q\left(\mathbb{Z}\right)=\prod _{i=1}^{m}q\left({\mathbb{Z}}_i\right).$

2.3. Bayesian Lasso Regularized Quantile Regression

Li and Zhu [25] proposed lasso regularized quantile regression. Specifically, they introduced a $L_1$ norm penalty term to control the complexity of the model, aiming to shrink some regression coefficients to zero. The objective function can be expressed as follows: Li and Zhu [25] proposed lasso regularized quantile regression. Specifically, they introduced an $L_1$ norm penalty term to control the complexity of the model, aiming to shrink some regression coefficients to zero. The objective function can be expressed as follows:

(10)$\underset{\beta }{min}\sum _{i=1}^n{\rho }_p\left(y_i-X_i^T\beta \right)+\lambda \sum _{i=1}^k\left|{\beta }_i\right|.$

Here, $\lambda >0$ is the regularization parameter used to adjust the balance between model complexity and the data. The term ${\rho }_p(·)$ represents the loss function proposed by Koenker and Bassett [1]. It is defined as:

${\rho }_p\left(y_i-X_i^T\beta \right)=\left\{\begin{array}{cc}p·\left(y_i-X_i^T\beta \right)\hfill & \mathrm{if}{y}_i-X_i^T\beta >0\hfill \\ -(1-p)·\left(y_i-X_i^T\beta \right)\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$

Bayesian lasso regularization for quantile regression was introduced by Li et al. [10], which improves upon the standard lasso regularization by incorporating a Bayesian framework. The Bayesian lasso regularization tends to capture outliers in the data and improve robustness. By introducing suitable priors on $\beta$, the solution to (10) is equivalent to the Bayesian maximum a posteriori (MAP) estimation.

In this approach, Laplace priors are used for the parameters $\beta$:

(11)$P(\beta \mid \tau ,\lambda )={(\tau \lambda /2)}^{p}exp\left\{-\tau \lambda \sum _{i=1}^k\left|{\beta }_i\right|\right\}.$

By using the equation:

(12)$\frac{a}{2}{e}^{-aH}={\int }_0^{\infty }\frac{1}{\sqrt{2\pi s}}exp\left(-\frac{z^2}{2s}\right)\frac{a^2}{2}exp\left(-\frac{a^{2}s}{2}\right)ds.$

where $\eta =p\lambda$, the Laplace prior on $\beta$ can be written as:

(13)$P(\beta \mid p,\lambda )=\prod _{i=1}^k{\int }_0^{\infty }\frac{1}{\sqrt{2\pi s_i}}exp\left(-\frac{{\beta }_i^2}{2s_i}\right)\frac{{\eta }^2}{2}exp\left(-\frac{{\eta }^2{s}_i}{2}\right)d{s}_i$

where $s=(s_1,\dots ,s_k)$ and $\beta =({\beta }_1,\dots ,{\beta }_k)$.

Considering the error term follows an asymmetric Laplace distribution (ALD), the posterior distribution of $\beta$ is given by:

(14)$P(\beta \mid y,X,p,\lambda )\propto exp\left\{-p\sum _{i=1}^n{\rho }_p\left(y_i-X_i^p\beta \right)-p\lambda \sum _{i=1}^k\left|{\beta }_i\right|\right\}$

To complete the model, gamma priors are placed on $\lambda$ and ${\eta }^2$, resulting in the following hierarchical model:

(15)$\left\{\begin{array}{c}{y}_i=X_i^{\tau }\beta +\phi z_i+\theta \sqrt{\sigma z_i}{v}_i\hfill \\ z_i\mid \sigma \sim exp(1/\sigma )=\sigma exp\left(-\sigma z_i\right)\hfill \\ v_i\sim N(0,1)=1/\sqrt{2\pi }exp\left(-v_i^2/2\right)\hfill \\ {\beta }_j\mid s_j\sim N\left(0,s_j\right)=\frac{1}{\sqrt{2\pi s_j}}exp\left(-{\beta }_j^2/2s_j\right)\hfill \\ s_j\mid {\eta }^2\sim exp\left({\eta }^2/2\right)=\frac{{\eta }^2}{2}exp\left(-\frac{{\eta }^2}{2}{s}_j\right)\hfill \\ {\lambda }_i{\eta }^2-{\lambda }^{i-1}exp(-b\lambda ){\left({\eta }^2\right)}^{c-1}exp\left(-d{\eta }^2\right).\hfill \end{array}\right.$

2.4. Missing Data

In the context of missing data, a comprehensive study was conducted by Little and Rubin [13], who provided insights into missing data mechanisms and patterns. Missing data mechanisms can be classified into three types: missing completely at random (MCAR), missing at random (MAR), and missing not at random (MNAR). In the following, we describe each mechanism and introduce the relevant parameters and their priors for modeling the missing data.

• (1). Missing Completely at Random (MCAR): These types of missing data are not related to either the observed data or the missing data. When $x_{ik}$ or $y_i$ is missing, the probability of missingness, denoted as $P(R_i\mid x_{ik},y_i,\varphi )$, follows a certain link function $H(·)$ parameterized by $\varphi$, i.e., $P(R_i\mid x_{ik},y_i,\varphi )=H(R_i\mid \varphi )$.

• (2). Missing at Random (MAR): These types of missing data are only related to the observed data. For example, when the covariate $x_i$ is missing, the probability of missingness $P(R_i\mid x_{ik},y_i,\varphi )$ is conditional on the observed response $y_i$ and follows a link function $H(·)$ parameterized by $\varphi$, i.e., $P(R_i\mid x_{ik},y_i,\varphi )=H(R_i\mid y_i,\varphi )$.

• (3). Missing Not at Random (MNAR): These types of missing data are related to the missing data itself. For example, when the covariate $x_{jk}$ is missing, the probability of missingness $P(R_i\mid x_{ik},y_i,\varphi )$ depends on the observed covariates $X_{ik}$ and follows a link function $H(·)$ parameterized by $\varphi$, i.e., $P(R_i\mid x_{ik},y_i,\varphi )=H(R_i\mid X_{ik},\varphi )$.

In the vast proportion of previous studies, the emphasis has been on random missing mechanisms. However, missing data frequently exhibit a connection to the data that are missing, suggesting nonignorable missingness. The study of such nonignorable missing data holds significant value as disregarding this data can potentially result in erroneous conclusions during analysis.

Furthermore, the missing data quantile regression models discussed in Section 3 and Section 4 are hierarchical models that can be represented using probabilistic-directed acyclic graphs (DAGs). In these graphs, the nodes correspond to the parameters in the model, and the arrows illustrate the dependencies between the parameters (Bishop [26]; Wasserman [27]). Specifically, in the DAG presented in this paper, the unshaded nodes correspond to the observed data.

Now, let us delve into Equation (3), where $y_i$ and $X_i$ denote the observed values of response variables and covariates for the ith observation. Assuming that both $y_i$ and $X_i$ are prone to nonignorable missingness, but not simultaneously missing, we can establish this model as a quantile regression model that accommodates nonignorable missing data. To effectively model the missing data mechanism, we introduce pertinent parameters and assign their prior distributions. For the random variable ${\delta }_i$, the optimal density function for it $q^{*}\left({\delta }_i\right)$.

$\begin{array}{cc}& \beta \sim N(0,{\sigma }_{\beta }^{2}I),x_{ik}\sim N\left({\mu }_x,{\sigma }_x^2\right),{\varphi }_0,{\varphi }_1\stackrel{\mathrm{ind}.}{\sim }N\left(0,{\sigma }_{\varphi }^2\right),\sigma \sim IG\left(A_{\sigma },B_{\sigma }\right)\hfill \\ & {\mu }_x\sim N\left(0,{\sigma }_{{\mu }_x}^2\right),{\sigma }_x^2\sim IG\left(A_x,B_x\right),z_i\sim exp\left(\frac{1}{\sigma }\right),v_i\sim N(0,1)\hfill \\ & \sigma \sim IG(a^{*},b),{\mu }_q\left({\delta }_i\right)=E_q\left({\delta }_i\right),{\mu }_q\left(1/{\delta }_i\right)=E_q\left(1/{\delta }_i\right)\hfill \\ & {\mu }_q\left(log\left({\delta }_i\right)\right)=E_q\left(log\left({\delta }_i\right)\right),{\sigma }_{q\left({\delta }_i\right)}^2={Var}_q\left({\delta }_i\right)\hfill \end{array}$

3. Variational Bayesian Inference

3.1. Missing Covariate Variables

In Equation (3), we define the matrix:

$\mathit{X}=\left[\begin{array}{cc}1& x_1\\ ⋮& ⋮\\ 1& x_n\end{array}\right]=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_n\end{array}\right],\mathit{Y}=\left[\begin{array}{cc}1& y_1\\ ⋮& ⋮\\ 1& y_n\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ ⋮\\ Y_n\end{array}\right].$

When the covariate x is missing, we introduce an indicator function $R={(R_1,\dots ,R_n)}^T$ to detect if $x_i$ is missing, defined as:

$R_i=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}{x}_i\mathrm{is}\mathrm{observed},\hfill \\ 0,\hfill & \mathrm{if}{x}_i\mathrm{is}\mathrm{missing}.\hfill \end{array}\right.$

Here, $x_{\mathrm{mis}}$ represents the $n_{\mathrm{mis}}$-dimensional subvector consisting of the missing $x_i$, and $n_{\mathrm{obs}}+n_{\mathrm{mis}}=n$. However, $x={(x_{\mathrm{obs}},x_{\mathrm{mis}})}^T$ may not correspond to the original data $x={(x_1,\dots ,x_n)}^T$ in terms of subscripts. Considering the nature of the indicator function, $R_i$ follows a $0-1$ distribution with probability $H\left(X_i\varphi \right)$, where $X_i$ denotes the ith row of matrix $\mathit{X}$. We adopt the probit regression model as the link function $H(·)$, thus:

$R_i\mid \varphi ,X_i\sim \mathsf{\Phi }\left(X_i^T\varphi \right).$

To implement the approach proposed by Albert and Chib [28], we introduce n independent latent variables $a_1,a_2,\dots ,a_n$ in the missing data mechanism, where $a_i$ follows a normal distribution with mean $X_i^T\varphi$ and variance 1. For $i=1,\dots ,n$, we define $R_i=1$ if $a_i>0$, and $R_i=0$ if $a_i<0$. The interactions between the regression parameters and the parameters of the missing data mechanism are shown in Figure 1.

3.2. Variational Bayesian Inference of QR Model Parameters

Based on the above parameters, we consider the mean field variational family,

$q\left(\beta ,z,\sigma ,x_{mis},{\mu }_x,{\sigma }_x^2,\varphi ,a\right)=q\left(\beta \right)q\left(z\right)q\left(\sigma \right)q\left(x_{mis}\right)q\left({\mu }_x\right)q\left({\sigma }_x^2\right)q\left(\varphi \right)q\left(a\right)$

In $z=\left(z_1,\dots ,z_n\right),a=\left(a_1,\dots ,a_n\right),x_{mis}=\left(x_{mis,1},\dots ,x_{mis,n_{mis}}\right)$, $z_1,\dots ,z_n$ are independent of each other, $a_1,\dots ,a_n$ are independent of each other, $x_{mis,1},\dots ,x_{mis,n_{mis}}$ are independent of each other, so there was:

$\begin{array}{c}q\left(z\right)=q\left(z_1\right)q\left(z_2\right)\cdots q\left(z_n\right)\\ q\left(a\right)=q\left(a_1\right)q\left(a_2\right)\cdots q\left(a_n\right)\\ q\left(x_{mis}\right)=q\left(x_{mis,1}\right)q\left(x_{mis,2}\right)\cdots q\left(x_{mis,n_{mis}}\right)\end{array}$

The optimal density function for each parameter is derived as follows. Please refer to Appendix A for a detailed derivation of the optimal density function for each parameter, as well as the calculations of $E_q\left(X_i{X}_i^T\right)$ and $E_q\left(X\right)$, where $l=1,\dots ,n_{\mathrm{mis}}$.

(16)$\begin{array}{cc}\hfill q\left(\beta \right)& \sim N\left({\mu }_{q\left(\beta \right)},{\mathsf{\Sigma }}_{q\left(\beta \right)}\right)\hfill \end{array}$

(17)$\begin{array}{cc}\hfill q\left(z_i\right)& \sim GIG\left(\frac{1}{2},a_{z_i},b_{z_i}\right)\hfill \end{array}$

(18)$\begin{array}{cc}\hfill q\left(\sigma \right)& \sim IG\left(a_{\sigma },b_{\sigma }\right)\hfill \end{array}$

(19)$\begin{array}{cc}\hfill q\left(x_{mis,l}\right)& \sim N\left({\mu }_{q\left(x_{mis,l}\right)},{\sigma }_{q\left(x_{mis,l}\right)}^2\right)\hfill \end{array}$

(20)$\begin{array}{cc}\hfill q\left({\mu }_x\right)& \sim N\left({\mu }_{q\left({\mu }_x\right)},{\sigma }_{q\left({\mu }_x\right)}^2\right)\hfill \end{array}$

(21)$\begin{array}{cc}\hfill q\left({\sigma }_x^2\right)& \sim IG\left(A_q,B_q\right)\hfill \end{array}$

(22)$\begin{array}{cc}\hfill q\left(\varphi \right)& \sim N\left({\mu }_{q\left(\varphi \right)},{\mathsf{\Sigma }}_{q\left(\varphi \right)}\right)\hfill \end{array}$

For the generalized inverse Gaussian distribution $q\left(z_i\right)\sim GIG\left(\frac{1}{2},a_{z_i},b_{z_i}\right)$, where $K_p(·)$ is the Bessel function with order p, with

(23)$\begin{array}{cc}\hfill {\mu }_{q\left(z_i\right)}& =\frac{\sqrt{b_{q\left(z_i\right)}}{K}_{3/2}\sqrt{a_{q\left(z_i\right)}{b}_{q\left(z_i\right)}}}{\sqrt{a_{q\left(z_i\right)}}{K}_{1/2}\sqrt{a_{q\left(z_i\right)}{b}_{q\left(z_i\right)}}}\hfill \end{array}$

(24)$\begin{array}{cc}\hfill {\mu }_{q\left(1/z_i\right)}& =\frac{\sqrt{a_{q\left(z_i\right)}}{K}_{3/2}\sqrt{a_{q\left(z_i\right)}{b}_{q\left(z_i\right)}}}{\sqrt{b_{q\left(z_i\right)}}{K}_{1/2}\sqrt{a_{q\left(z_i\right)}{b}_{q\left(z_i\right)}}}-\frac{1}{b_{q\left(z_i\right)}}\hfill \end{array}$

(25)$\begin{array}{cc}\hfill {\mu }_{q\left(log\left(z_i\right)\right)}& =log\frac{\sqrt{b_{q\left(z_i\right)}}}{\sqrt{a_{q\left(z_i\right)}}}+\frac{\partial }{\partial p}log{K}_p\left(\sqrt{a_{q\left(z_i\right)}{b}_{q\left(z_i\right)}}\right)\hfill \end{array}$

For the inverse gamma distribution $q\left(\sigma \right)\sim IG\left(a_{\sigma },b_{\sigma }\right)$, where $\psi (·)$ denotes the inverse gamma function, we have

(26)$\begin{array}{cc}\hfill {\mu }_{q\left(\sigma \right)}& =\frac{b_{\sigma }}{a_{\sigma }-1}\hfill \\ \hfill {\sigma }_{q\left(1/\sigma \right)}^2& =\frac{a_{\sigma }}{b_{\sigma }}\hfill \\ \hfill {\mu }_{q(log(\sigma \left)\right)}& =log\left(b_{q\left(\sigma \right)}\right)-\psi \left(b_{q\left(\sigma \right)}\right)\hfill \end{array}$

In addition, the required expectations for the additional variable a are as follows:

(27)$\begin{array}{c}\hfill {\mu }_{q\left(a\right)}=E_{q\left(x_{\mathrm{mis}}\right)}\left(X\right){\mu }_{q\left(\varphi \right)}+\frac{{\left(2\pi \right)}^{-2}(2R-1)\odot exp\{-\frac{1}{2}{\left(E_{q\left(x_{\mathrm{mis}}\right)}\left(X\right){\mu }_{q\left(\varphi \right)}\right)}^2\}}{\mathsf{\Phi }((2R-1)\odot \left(E_{q\left(x_{\mathrm{mis}}\right)}\left(X\right){\mu }_{q\left(\varphi \right)})\right)}\end{array}$

In the model, the evidence lower bound by

(28)$\begin{array}{cc}L\left(p\right)\hfill & =\frac{1}{2}log\frac{\left|{\mathsf{\Sigma }}_{q\left(\beta \right)}\right|}{{\sigma }_{\beta }^2}-\frac{1}{2{\sigma }_{\beta }^2}\left\{{∥{\mu }_{q\left(\beta \right)}∥}^2+tr\left({\mathsf{\Sigma }}_{q\left(\beta \right)}\right)\right\}+n{\mu }_{q(log\sigma )}-\frac{1}{2}{\displaystyle \sum _{i=1}^n}\left(a_{z_i}{\mu }_{q\left(z_i\right)}+b_{z_i}{\mu }_{q\left(1/z_i\right)}\right)\hfill \\ \hfill & +{\displaystyle \sum _{i=1}^n}log\frac{{\left(a_{z_i}/b_{z_i}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{z_i}{b}_{z_i}}\right)}-\frac{1}{2}{\displaystyle \sum _{i=1}^n}{\mu }_{q\left(log\left(z_i\right)\right)}+\frac{1}{2}log\frac{{\sigma }_{q\left({\mu }_x\right)}^2}{{\sigma }_{{\mu }_x}^2}-\frac{{\mu }_{q\left({\mu }_x\right)}^2+{\sigma }_{q\left({\mu }_x\right)}^2}{2{\sigma }_{{\mu }_x}^2}\hfill \\ \hfill & +\frac{1}{2}{∥E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\varphi \right)}∥}^2-\frac{1}{2}tr\left\{E_{q\left(x_{mi}\right)}\left(X^{T}X\right)\left({\mu }_{q\left(\varphi \right)}{\mu }_{q\left(\varphi \right)}{}^T+{\mathsf{\Sigma }}_{q\left(\varphi \right)}\right)\right\}\hfill \\ \hfill & +R^{T}log\left\{\mathsf{\Phi }\left(E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\varphi \right)}\right)\right\}+{(1-R)}^{T}log\left\{1-\mathsf{\Phi }\left(E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\varphi \right)}\right)\right\}\hfill \\ \hfill & +\frac{1}{2}log\frac{\left|{\sum }_{q\left(\varphi \right)}\right|}{{\sigma }_{\varphi }^2}-\frac{1}{2{\sigma }_{\varphi }^2}\left\{{∥{\mu }_{q\left(\varphi \right)}∥}^2+tr\left({\mathsf{\Sigma }}_{q\left(\varphi \right)}\right)\right\}+\frac{1}{2}{1}^T{\mu }_{q\left(x_{mis}\right)}-A_{q\left({\sigma }_x^2\right)}log\left(B_{q\left({\sigma }_x^2\right)}\right)\hfill \\ \hfill & +\left(a_{q\left(\sigma \right)}-a\right){\mu }_{q(log(\sigma \left)\right)}+\frac{\left(b_{q\left(\sigma \right)}-b\right)a_{q\left(\sigma \right)}}{b_{q\left(\sigma \right)}}+log\mathsf{\Gamma }\left(a_{q\left(\sigma \right)}\right)-a_{q\left(\sigma \right)}log{b}_{q\left(\sigma \right)}\hfill \\ \hfill & -\frac{{\sum }_{i=1}^n{\mu }_{q\left(log{z}_i\right)}}{2}-\frac{n}{2}{\mu }_{q(log\sigma )}-\frac{1}{2{\theta }^2}{\mu }_{q(1/\sigma )}{1}^T{\mu }_{q(1/z)}{E}_{q\left(x_{mis}\right)q\left(\beta \right)}{∥y-X^T\beta ∥}^2\hfill \\ \hfill & +\frac{\phi {\mu }_{q(1/\sigma )}}{{\theta }^2}{\displaystyle \sum _{i=1}^n}{y}_i-\frac{\phi {\mu }_{q(1/\sigma )}}{{\theta }^2}{1}^T{E}_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\beta \right)}-\frac{{\phi }^2{\mu }_{q(1/\sigma )}}{2{\theta }^2}{1}^T{\mu }_{q\left(z\right)}\hfill \end{array}$

where:

(29)$\begin{array}{cc}\hfill E_{q\left(x_{mis}\right)q\left(\beta \right)}{∥y-X^T\beta ∥}^2& ={\parallel y\parallel }^2-2y^T{E}_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\beta \right)}\hfill \\ \hfill & +tr\left\{E_{q\left(x_{m\dot{m}}\right)}\left(X^{T}X\right)\left({\mathsf{\Sigma }}_{q\left(\beta \right)}+{\mu }_{q\left(\beta \right)}{\mu }_{q\left(\beta \right)}{}^T\right)\right\}\hfill \end{array}$

The specific algorithm is shown in Algorithm 1.

Algorithm 1: VI with Nonignorable Missing Covariates in QR Models.

Input: The prior distribution of each parameter, including ${\sigma }_{\beta }^2,{\mu }_x,{\sigma }_x^2,{\sigma }_{\varphi }^2,{\sigma }_{{\mu }_x}^2,A_x,B_x,a^{*},b$, and a given lower bound on the evidence lower bound $tol={10}^{-5}$

3.3. Missing Response Variables

For Equation (3), let us consider the scenario where the response variable $y=\left(y_1,\dots ,y_n\right)$ contains missing values. Similar to missing covariates, we introduce an indicator function $R=\left(R_1,\dots ,R_n\right)$ to detect missing response values, defined as:

$R_i=\left\{\begin{array}{cc}1,\hfill & \mathrm{if}{y}_i\mathrm{is}\mathrm{observed},\hfill \\ 0,\hfill & \mathrm{if}{y}_i\mathrm{is}\mathrm{missing}.\hfill \end{array}\right.$

Here, $y_{mis}$ represents the $n_{mis}$-dimensional subvector consisting of the missing $y_i$, and we have $n_{obs}+n_{mis}=n$. However, $y=\left(y_{obs},y_{mis}\right)$ does not necessarily correspond to the original data $y=\left(y_1,\dots ,y_n\right)$ in terms of subscripts.

Similar to Section 3.1, we introduce the indicator function $R_i$ and consider the link function as a logistic regression model. Consequently, we have:

$R_i\mid \varphi ,X_i\sim \mathrm{Bernoulli}\left(\frac{1}{1+exp(-X_i^T\varphi )}\right).$

Here, $X_i$ represents the i-th row of the matrix $\mathbf{X}$, and $\varphi$ is the parameter vector.

$P\left(R_i\mid \varphi ,Y_i\right)=1/\left(1+exp\left(-Y_i^T\varphi \right)\right)$

Auxiliary Variables

Polson et al. [29] proposed a Bayesian estimation method for logistic regression models by introducing auxiliary variables $\omega =\left({\omega }_1,\dots ,{\omega }_n\right)$ that follow the Pòlya-Gamma distribution. This approach addresses the challenges associated with complex integration and sampling using the Metropolis–Hastings (MH) algorithm.

In their methodology, they consider the following model:

$y_i^{*}\sim \mathrm{Binom}\left(n_i^{*},\frac{1}{1+e^{-{\psi }_i^{*}}}\right)$

where ${\psi }_i^{*}=x_i^{*T}{\beta }^{*}$ and $i=1,\dots ,N^{*}$. They place a Gaussian prior on ${\beta }^{*}$, denoted as ${\beta }^{*}\sim \mathrm{N}(b^{*},B)$. To incorporate the auxiliary variable $\omega$, they introduce it into the model.

The posterior density, after adding the auxiliary variable $\omega$, can be expressed as follows:

$P({\beta }^{*},\omega |y^{*},x^{*})\propto P\left(y^{*}\right|{\beta }^{*},x^{*})·P\left({\beta }^{*}\right)·P\left(\omega \right)·P\left(x^{*}\right).$

Here, $P\left(y^{*}\right|{\beta }^{*},x^{*})$ represents the likelihood function, $P\left({\beta }^{*}\right)$ is the prior distribution of ${\beta }^{*}$, $P\left(\omega \right)$ is the distribution of the auxiliary variable $\omega$, and $P\left(x^{*}\right)$ denotes the prior distribution of the covariates $x^{*}$.

The specific form of the posterior density and the method of inference depend on the choice of prior distributions, likelihood function, and the implementation of the Pòlya-Gamma auxiliary variables.

$\begin{array}{cc}& \left({\omega }_i\mid {\beta }^{*}\right)\sim PG\left(n_i^{*},x_i^{*T}{\beta }^{*}\right)\hfill \\ & \left(\varphi \mid y^{*},\omega \right)\sim N\left(m_{\omega },V_{\omega }\right)\hfill \end{array}$

where

$\begin{array}{cc}& V_{\omega }={\left(X^{*T}\mathsf{\Omega }{X}^{*}+B^{-1}\right)}^{-1}\hfill \\ & m_{\omega }=V_{\omega }\left(X^{*T}\kappa +B^{-1}{b}^{*}\right).\hfill \end{array}$

We define $\kappa =\left(y_1^{*}-n_1^{*}/2,\dots ,y_{N^{*}}^{*}-n_{N^{*}}^{*}/2\right)$, where $\mathsf{\Omega }$ represents the diagonal matrix with ${\omega }_1$ as the diagonal elements. Here, $i=1,\dots ,N^{*}$.

Variational Bayesian inference is typically applicable to specific types of models, particularly conditional covariance exponential family models. However, logistic regression models with Gaussian priors are an exception because there is no covariance between the logistic likelihood and the Gaussian prior in such models. To address this limitation, Daniele et al. [30] proposed introducing an additional variable $\omega =\left({\omega }_1,\dots ,{\omega }_n\right)$ that follows the Pòlya-Gamma distribution. By incorporating this variable into models with logistic components, it strengthens the covariance between the logistic likelihood and the Gaussian prior. This enables Bayesian framework statistical inference for this class of models using a variational Bayesian algorithm.

In the context of this chapter, we consider the missing data mechanism as a logistic regression model. Consequently, we introduce additional variables that follow the Pòlya-Gamma distribution, denoted as $\omega =\left({\omega }_1,\dots ,{\omega }_n\right)$. The interaction between the parameters is shown in Figure 2.

3.4. Variational Bayesian Inference of QR Model

For the above parameters, we consider the mean field variational family,

$q\left(\beta ,z,\sigma ,y_{mis},\varphi ,\omega \right)=q\left(\beta \right)q\left(z\right)q\left(\sigma \right)q\left(y_{mis}\right)q\left(\varphi \right)q\left(\omega \right).$

Consider the variables $z=\left(z_1,\dots ,z_n\right)$, $a=\left(a_1,\dots ,a_n\right)$, and $y_{\mathrm{mis}}=\left(y_{\mathrm{mis},1},\dots ,y_{\mathrm{mis},n_{\mathrm{mis}}}\right)$. Here, $z_1,\dots ,z_n$ are independent of each other, $a_1,\dots ,a_n$ are independent of each other, and $y_{\mathrm{mis},1},\dots ,y_{\mathrm{mis},n_{\mathrm{mis}}}$ are independent of each other. This follows a similar pattern as discussed in Section 3.1.

$\begin{array}{c}q\left(z\right)=q\left(z_1\right)q\left(z_2\right)\cdots q\left(z_n\right)\\ q\left(a\right)=q\left(a_1\right)q\left(a_2\right)\cdots q\left(a_n\right)\\ q\left(y_{mis}\right)=q\left(y_{mis,1}\right)q\left(y_{mis,2}\right)\cdots q\left(y_{mis,n_{mis}}\right).\end{array}$

In the scenario of missing response variables, the posterior probabilities $P(\beta \mid {\mathsf{\Theta }}_r)$, $P(z\mid {\mathsf{\Theta }}_r)$, and $P(\sigma \mid {\mathsf{\Theta }}_r)$ remain the same as in the presence of missing covariates and are not reiterated here. However, we provide the full conditional probability posterior for $P(\varphi \mid {\theta }_r)$, $P(\omega \mid {\theta }_r)$, and $P\left(y_{\mathrm{mix}}\mid {\theta }_r\right)$. Following the approach of Polson et al. [29], we consider the prior distribution of ${\omega }_i$ as $PG(1,0)$.

$\begin{array}{cc}& P\left({\omega }_i\mid {\mathsf{\Theta }}_r\right)\sim PG\left(R_i,Y_i^T\varphi \right)\hfill \\ & P(\varphi \mid {\mathsf{\Theta }}_r)\sim N\left({\mu }_{\varphi },{\mathsf{\Sigma }}_{\varphi }\right)\hfill \end{array}$

where

$\begin{array}{cc}& {\mathsf{\Sigma }}_{\varphi }={\left(Y^T\mathsf{\Omega }Y+{\left({\sigma }_{\varphi }^{2}I\right)}^{-1}\right)}^{-1}\hfill \\ & {\mu }_{\varphi }={\mathsf{\Sigma }}_{\varphi }\left(Y^{T}K\right)\hfill \end{array}$

$\kappa =\left(R_1-1/2,\dots ,R_n-1/2\right)$, $\mathsf{\Omega }$ denotes the diagonal matrix with ${\omega }_i$ diagonal elements, $i=1,\dots ,n$. The optimal density of each parameter in the missing data mechanism section is as follows:

(30)$q\left({\omega }_i\right)\sim PG\left(1,{\xi }_i\right)$

(31)$q\left(\varphi \right)\sim N\left({\mu }_{q\left({\varphi }^y\right)},{\mathsf{\Sigma }}_{q\left({\varphi }^y\right)}\right)$

where:

(32)${\xi }_i=\sqrt{E_q{\left(Y_i\right)}^T\left({\sigma }_{\varphi }^2\right)I{E}_q\left(Y_i\right)+\left(E_q{\left(Y_i\right)}^T{\mu }_{q\left(\varphi \right)}\right)}$

(33)${\mathsf{\Sigma }}_{q\left(\varphi \right)}={\left(E_q{\left(Y\right)}^T\mathsf{\Omega }{E}_q\left(Y\right)+{\left({\sigma }_{\varphi }^{2}I\right)}^{-1}\right)}^{-1}$

(34)${\mu }_{q\left(\varphi \right)}={\mathsf{\Sigma }}_{q\left(\varphi \right)}\left(E_q{\left(Y\right)}^T\kappa \right)$

The optimal density functions of the remaining parameters are as follows, and the detailed derivation of the optimal density function for each parameter and $E_q\left(Y\right)$ is similar to $E_q\left(X\right)$.

(35)$q\left(\beta \right)\sim N\left({\mu }_{q\left(\beta \right)},{\mathsf{\Sigma }}_{q\left(\beta \right)}\right)$

(36)$q\left(z_i\right)\sim GIG\left(\frac{1}{2},a_{z_i},b_{z_i}\right)$

(37)$q\left(\sigma \right)\sim IG\left(a_{\sigma },b_{\sigma }\right)$

(38)$q\left(y_{mis,l}\right)\sim N\left({\mu }_{q\left(y_{mis,l}\right)},{\sigma }_{q\left(y_{mis,l}\right)}^2\right)$

In the model, the evidence lower bound by

(39)$\begin{array}{cc}\hfill L\left(q\right)& =\frac{1}{2}log\frac{\left|{\mathsf{\Sigma }}_{q\left(\beta \right)}\right|}{{\sigma }_{\beta }^2}-\frac{1}{2{\sigma }_{\beta }^2}\left\{{∥{\mu }_{q\left(\beta \right)}∥}^2+{\mathsf{\Sigma }}_{q\left(\beta \right)}\right\}-n{\mu }_{q(log(\sigma \left)\right)}-\frac{1}{2}\sum _{j=1}^{n_{mis}}log{\sigma }_{q\left(y_{mis,j}\right)}^2\hfill \\ \hfill & -\sum _{i=1}^n\left({\mu }_{q\left(z_i\right)}{\mu }_{q(1/\sigma )}+log\frac{{\left(a_{z_i}/b_{z_i}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{z_i}{b}_{z_i}}\right)}-\frac{1}{2}{\mu }_{q\left(log{z}_i\right)}-\frac{1}{2}\left(a_{z_i}{\mu }_{q\left(z_i\right)}+b_{z_i}{\mu }_{q\left(1/z_i\right)}\right)\right)\hfill \\ \hfill & +\left(a_{q\left(\sigma \right)}-a\right){\mu }_{q(log(\sigma \left)\right)}+\frac{\left(b_{q\left(\sigma \right)}-b\right)a_{q\left(\sigma \right)}}{b_{q\left(\sigma \right)}}-a_{q\left(\sigma \right)}log{b}_{q\left(\sigma \right)}\frac{1}{2}log\left|{\sigma }_{q\left(\varphi \right)}^{2}I\right|-\frac{1}{2}{\mu }_{q\left(\varphi \right)}^T{\left({\mathsf{\Sigma }}_{\varphi }\right)}^{-1}{\mu }_{q\left(\varphi \right)}\hfill \\ \hfill & +\sum _{i=1}^n\left[\left(R_i-\frac{1}{2}\right)E_q\left(Y_i\right){\mu }_{q\left(\varphi \right)}-log\left(\frac{exp\left({\xi }_i\right)}{1+exp\left({\xi }_i\right)}\right)-\frac{1}{2}{\xi }_i\right]-\frac{1}{2}tr\left({\left({\mathsf{\Sigma }}_{\varphi }\right)}^{-1}{\mathsf{\Sigma }}_{q\left(\varphi \right)}\right)\hfill \\ \hfill & -\frac{{\mu }_{q(1/\sigma )}}{2}\sum _{i=1}^n\left[E_q\left(y_i^2\right){\mu }_{q\left(1/{\tilde{z}}_i\right)}/{\theta }^2-\frac{2}{{\theta }^2}{E}_q\left(y_i\right)X_i^T{\mu }_{q\left(\beta \right)}{\mu }_{q\left(1/z_i\right)}-2E_q\left(y_i\right)\phi /{\theta }^2\right]\hfill \\ \hfill & -\frac{{\mu }_{q(1/\sigma )}}{2}\sum _{i=1}^n\left[\frac{\phi }{{\theta }^2}{X}_i^T{\mu }_{q\left(\beta \right)}+\frac{{\mu }_{q\left(1/z_i\right)}}{{\theta }^2}tr\left(X_i^T{X}_i\left({\mathsf{\Sigma }}_{q\left(\beta \right)}+{\mu }_{q\left(\beta \right)}{\mu }_{q\left(\beta \right)}^T\right)\right)+\frac{{\phi }^2}{{\theta }^2}{\mu }_{q\left(z_i\right)}\right]\hfill \\ \hfill & -\frac{1}{2}\sum _{i=1}^n{\mu }_{q\left(log\left(z_i\right)\right.}-\frac{n}{2}{\mu }_{q(log(\sigma \left)\right).}\hfill \end{array}$

The specific algorithm is shown in Algorithm 2.

Algorithm 2: VI with Nonignorable Missing Response in QR Models.

Input: The prior distribution of each parameter, including ${\sigma }_{\beta }^2,{\sigma }_{\varphi }^2,{\mu }_x,{\sigma }_x^2,a^{*},b$, and a given lower bound on the evidence lower bound $tol={10}^{-5}$

4. Bayesian Variable Selection of the QR Model

4.1. Covariates with Nonignorable Missing Variable Choices

For Equation (15), the following matrix is defined:

$X=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1k}\\ ⋮& ⋮& \ddots & ⋮\\ 1& x_{1n}& \cdots & x_{nk}\end{array}\right]=\left[\begin{array}{c}{X}_1\\ ⋮\\ X_n\end{array}\right]Y=\left[\begin{array}{cc}1& y_1\\ ⋮& ⋮\\ 1& y_n\end{array}\right]=\left[\begin{array}{c}{Y}_1\\ ⋮\\ Y_n\end{array}\right].$

In the context of our analysis, the variable $x_i$ represents the covariate for the i-th observation, while $y_i$ denotes the response variable corresponding to the same observation. When any of the covariates in X, specifically the j-th dimension covariate where j ranges from 1 to k, is missing, we follow a similar approach as described in Section 3.1. We introduce the indicative function $R={\left(R_1,\dots ,R_n\right)}^T$ to identify the missing covariates and consider the missing data mechanism as a probit regression.

Furthermore, to incorporate the variables’ choices, we introduce the parameters s and ${\eta }^2$. The prior distribution for s follows an exponential distribution, specifically $s\sim exp\left({\eta }^2/2\right)$. On the other hand, the parameter ${\eta }^2$ is assigned a prior distribution with a gamma distribution, denoted as ${\eta }^2\sim Ga(c,d)$. The parameter interactions are shown in Figure 3.

The optimal densities of the parameters $z,\sigma ,{\mu }_x,{\sigma }_x^2,a,\varphi$ are similar to those of Algorithm 1 and will not be repeated. In the following we give the optimal density functions for the parameters $\beta$,s,${\eta }^2$,$x_{j,mis,l}$, where $l=1,\cdots ,n_{mis}$.

• (1). The full conditional distribution of the parameter $\beta$ is as follows:

  $\begin{array}{cc}\hfill P\left({\beta }_j\mid {\mathsf{\Theta }}_r\right)& \propto P\left(y\mid X,{\beta }_j,\sigma ,z\right)P\left({\beta }_j\mid s_j\right)\hfill \\ & \propto exp\left\{-\frac{1}{2}\left[\left(\sum _{i=1}^n\frac{x_{j,i}^2}{{\theta }^2{z}_i\sigma }+\frac{1}{s_k}\right){\beta }_j^2-2\sum _{i=1}^n\frac{y_{i,j}{x}_{j,i}}{{\theta }^2{z}_i\sigma }{\beta }_j\right]\right\}\hfill \end{array}$

  where $x_k={\left(x_{k1},\dots ,x_{kn}\right)}^T,y_{j,i}=y_i-\phi z_i-{\sum }_{i=1,i\ne j}^k{x}_{i,-j}{\beta }_{-j}$. Therefore, the full conditional distribution of the parameter $\beta$ is $N\left({\tilde{\mu }}_{\beta },{\tilde{\sigma }}_{\beta }^2\right)$, and ${\tilde{\mu }}_{\beta }={\sigma }_{\beta }^2{\sum }_{i=1}^n\frac{{\tilde{y}}_{i,j}{x}_{j,i}}{{\theta }^2{z}_i\sigma },{\tilde{\sigma }}_{\beta }^2=1/\left({\sum }_{i=1}^n\frac{x_{j,i}^2}{{\theta }^2{z}_i\sigma }+\frac{1}{s_j}\right)$. From Equation (7), we know that the parameter ${\beta }_j$ optimal density function is normally distributed as follows:

  (40)$\begin{array}{c}\hfill q\left({\beta }_j\right)\sim N\left({\mu }_{q\left({\beta }_j\right)},{\sigma }_{q\left({\beta }_j\right)}^2\right)\end{array}$

  where

  $\begin{array}{cc}& {\sigma }_{q\left({\beta }_j\right)}^2=1/\left(\sum _{i=1}^n\frac{{\mu }_{q(1/\sigma )}{\mu }_{q(1/z_i)}{E}_{q\left(x_{mis}\right)}{\left(X_j\right)}^2}{{\theta }^2}+{\mu }_{q(1/z_j)}\right)\hfill \\ & {\mu }_{q\left({\beta }_j\right)}={\sigma }_{q\left({\beta }_j\right)}^2\left[\frac{{\mu }_{q(1/\sigma )}{E}_{q\left(x_{mis}\right)}\left(X_j\right)}{{\theta }^2}\left(y\odot {\mu }_{q(1/z)}-\phi \mathrm{l}-\left(E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\beta \right)}-E_{q\left(x_{mis}\right)}\left(X_j\right){\mu }_{q\left({\beta }_j\right)}\right)\odot {\mu }_{q(1/z)}\right)\right].\hfill \end{array}$

• (2). The full conditional distribution of the parameter s is as follows:

  $\begin{array}{cc}\hfill P\left(s_j\mid {\mathsf{\Theta }}_r\right)& \propto P\left({\beta }_j\mid s_j\right)P\left(s_j\mid {\eta }^2\right)\hfill \\ & \propto \frac{1}{\sqrt{s_j}}exp\left\{-\frac{1}{2}\left[{\eta }^2{s}_j+{\beta }_j^2{s}_j^{-1}\right]\right\}.\hfill \end{array}$

  Therefore, the full conditional distribution of the parameter $s_j$ is $GIG\left(1/2,{\eta }^2,{\beta }_k^2\right)$. From Equation (7), we know that the parameter $s_j$ optimal density function is a generalized inverse Gaussian distribution

  (41)$\begin{array}{c}\hfill q\left(s_j\right)\sim GIG\left(\frac{1}{2},a_{s_j},b_{s_j}\right)\end{array}$

  where $a_{s_j}={\mu }_{q\left({\eta }^2\right)},b_{s_j}={\mu }_{q\left({\beta }_j\right)}^2+{\left({\sigma }_{q\left({\beta }_j\right)}^2\right)}^2$.

• (3). The full conditional distribution of the parameter ${\eta }^2$ is as follows:

  $\begin{array}{cc}\hfill P\left({\eta }^2\mid {\mathsf{\Theta }}_r\right)& \propto P\left(s\mid {\eta }^2\right)P\left({\eta }^2\right)\hfill \\ & \propto {\left({\eta }^2\right)}^{k+c-1}exp\left\{-{\eta }^2\left(\sum _{i=1}^k\frac{s_i}{2}+d\right)\right\}.\hfill \end{array}$

  The optimal density function is known to be the gamma distribution

  (42)$\begin{array}{c}\hfill q\left({\eta }^2\right)\sim Ga\left(c_q,d_q\right)\end{array}$

  where $c_q=k+c,d_q=\frac{1}{2}{\sum }_{j=1}^k{\sigma }_{q\left({\beta }_j\right)}^2+d$

• (4). The full conditional distribution of $x_{l,mis}$ is as follows:

  $\begin{array}{cc}& P\left(x_{j,mis,l}\mid {\mathsf{\Theta }}_r\right)\propto P(y\mid x,\beta ,z,\sigma )P\left(x_{j,mis,l}\mid {\mu }_x,{\sigma }_x^2\right)P(a\mid \varphi ,x)\hfill \\ & \propto \frac{exp\left\{-\frac{1}{2}\left[x_{j,mis,l}^2\left(\frac{{\beta }_j^2}{{\theta }^2\sigma z_{mis,l}}+\frac{1}{{\sigma }_x^2}+{\varphi }_1^2\right)\right]\right\}}{exp\left\{-\frac{1}{2}\left[-2x_{j,mis,l}\left(-\frac{{\beta }_j{X}_{l,-j}{\beta }_{-j}+{\beta }_j\phi z_{mis,l}-y_{mis,l}{\beta }_j}{{\theta }^2\sigma z_{mis,l}}+\frac{{\mu }_x}{{\sigma }_x^2}+a_{mis,l}{\varphi }_1+{\varphi }_0{\varphi }_1\right)\right]\right\}}.\hfill \\ \end{array}$

The above equation, $X=\left(x_1,\dots ,x_k\right),x_j={\left(x_{j,obs},x_{j,mis}\right)}^T,X_{l,-j}=\left(x_{l,1},\dots ,x_{l,j-1},\right.$$\left.x_{l,j+1},\dots ,x_{l,k}\right)$, ${\beta }_{-j}=\left({\beta }_1,\dots ,{\beta }_{j-1},{\beta }_{j+1},\dots ,{\beta }_k\right)$. The $x_{j,mis,l}$ optimal density function can be found as

(43)$\begin{array}{c}\hfill q\left(x_{j,mis,l}\right)\sim N({\mu }_{q\left(x_{j,mis,l}\right)},{\sigma }_{q\left(x_{j,mis,l}\right)}^2)\end{array}$

where:

$\begin{array}{cc}\hfill {\sigma }_{q\left(x_{j,mis,l}\right)}^2& =1/\left({\mu }_{q\left(1/{\sigma }_x^2\right)}+\frac{{\mu }_{q(1/\sigma )}{\mu }_{q\left(1/z_{mis,l}\right)}}{{\theta }^2}\left({\mu }_{q\left({\beta }_j\right)}^2+{\left({\sigma }_{q\left({\beta }_j\right)}^2\right)}^2\right)+{\mu }_{q\left({\varphi }_1\right)}^2+{\left({\mathsf{\Sigma }}_{q\left(\varphi \right)}\right)}_{22}^2\right)\hfill \\ \hfill {\mu }_{q\left(x_{j,mis,}\right)}& =-\frac{{\sigma }_{q\left(x_{j,mis,l}\right)}^2{\mu }_{q(1/\sigma )}{\mu }_{q\left(1/z_{mis,l}\right)}}{{\theta }^2}\left({\mu }_{q\left({\beta }_j\right)}{X}_{l,-j}{\mu }_{q\left({\beta }_{-j}\right)}+{\mu }_{q\left({\beta }_j\right)}\phi {\mu }_{q\left(1/z_{mis,l}\right)}-y_{mis,l}{\mu }_{q\left({\beta }_j\right)}\right)\hfill \\ & +{\sigma }_{q\left(x_{j,mis,l}\right)}^2\left({\mu }_{q\left({\mu }_x\right)}/{\mu }_{q\left({\sigma }_x^2\right)}+{\mu }_{q\left(a_{mis,j}\right)}{\mu }_{q\left({\varphi }_1\right)}+{\mu }_{q\left({\varphi }_0\right)}{\mu }_{q\left({\varphi }_1\right)}\right).\hfill \\ \end{array}$

In the model, the evidence lower bound by

(44)$\begin{array}{cc}\hfill l\left(p\right)& =\frac{1}{2}log\frac{\left|{\mathsf{\Sigma }}_{q\left(\widehat{\beta }\right)}\right|}{{\sigma }_{\beta }^2}-\frac{1}{2{\sigma }_{\beta }^2}\left\{{∥{\mu }_{q\left(\widehat{\beta }\right)}∥}^2+tr\left({\mathsf{\Sigma }}_{q\left(\widehat{\beta }\right)}\right)\right\}+n{\mu }_{q(log\widehat{\sigma })}-\frac{1}{2}\sum _{i=1}^n\left(a_{{\widehat{z}}_i}{\mu }_{q\left({\widehat{z}}_i\right)}+b_{z_i}{\mu }_{q\left(1/{\widehat{z}}_i\right)}\right)\hfill \\ \hfill & +\sum _{i=1}^{n}log\frac{{\left(a_{{\widehat{z}}_i}/b_{{\widehat{z}}_i}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{{\tilde{z}}_i}{b}_{{\widehat{z}}_i}}\right)}-\frac{1}{2}\sum _{i=1}^n{\mu }_{q\left(log\left({\tilde{z}}_i\right)\right)}+\frac{1}{2}log\frac{{\sigma }_{q\left({\mu }_x\right)}^2}{{\sigma }_{{\mu }_x}^2}-\frac{{\mu }_{q\left({\mu }_x\right)}^2+{\sigma }_{q\left({\mu }_x\right)}^2}{2{\sigma }_{{\mu }_k}^2}\hfill \\ \hfill & +\frac{1}{2}{∥E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\widehat{\varphi }\right)}∥}^2-\frac{1}{2}tr\left\{E_{q\left(x_{mis}\right)}\left(X^{T}X\right)\left({\mu }_{q\left(\widehat{\varphi }\right)}{\mu }_{q\left(\widehat{\varphi }\right)}{}^T+{\mathsf{\Sigma }}_{q\left(\widehat{\varphi }\right)}\right)\right\}\hfill \\ \hfill & +R^{T}log\left\{\mathsf{\Phi }\left(E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\widehat{\varphi }\right)}\right)\right\}+{(1-R)}^{T}log\left\{1-\mathsf{\Phi }\left(E_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\widehat{\varphi }\right)}\right)\right\}\hfill \\ \hfill & +\frac{1}{2}log\frac{\left|{\sum }_{q\left(\widehat{\varphi }\right)}\right|}{{\sigma }_{\varphi }^2}-\frac{1}{2{\sigma }_{\varphi }^2}\left\{{∥{\mu }_{q\left(\widehat{\varphi }\right)}∥}^2+tr\left({\mathsf{\Sigma }}_{q\left(\widehat{\varphi }\right)}\right)\right\}+\frac{1}{2}{1}^T{\mu }_{q\left(x_{m\dot{\mathrm{z}}}\right)}-A_{q\left({\sigma }_x^2\right)}log\left(B_{q\left({\sigma }_x^2\right)}\right)\hfill \\ \hfill & +\left(a_{q\left(\sigma \right)}-a\right){\mu }_{q(log(\sigma \left)\right)}+\frac{\left(b_{q\left(\sigma \right)}-b\right)a_{q\left(\sigma \right)}}{b_{q\left(\widehat{\sigma }\right)}}+log\mathsf{\Gamma }\left(a_{q\left(\sigma \right)}\right)-a_{q\left(\sigma \right)}log{b}_{q\left(\sigma \right)}\hfill \\ \hfill & -\frac{1}{2}\sum _{i=1}^n{\mu }_{q\left(log{z}_i\right)}-\frac{n}{2}{\mu }_{q(log\widehat{\sigma })}-\frac{1}{2{\theta }^2}{\mu }_{q(1/\widehat{\sigma })}{1}^T{\mu }_{q(1/\widehat{z})}{E}_{q\left(x_{mis}\right)q\left(\widehat{\beta }\right)}{∥y-X^T\beta ∥}^2\hfill \\ \hfill & +\frac{\phi {\mu }_{q(1/\sigma )}}{{\theta }^2}\sum _{i=1}^n{y}_i-\frac{\phi {\mu }_{q(1/\sigma )}}{{\theta }^2}{1}^T{E}_{q\left(x_{mis}\right)}\left(X\right){\mu }_{q\left(\beta \right)}-\frac{{\phi }^2{\mu }_{q(1/\sigma )}}{2{\theta }^2}{1}^T{\mu }_{q\left(z\right)}\hfill \\ \hfill & -\sum _{j=1}^k\left(-\frac{1}{2}{\mu }_{q\left(s_j\right)}{\mu }_{q\left({\eta }^2\right)}-log\frac{{\left(a_{s_j}/b_{s_j}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{s_j}{b}_{s_j}}\right)}+\frac{1}{2}{\mu }_{q\left(s_j\right)}\left(log{s}_j\right)+\frac{1}{2}\left(a_{s_j}{\mu }_{q\left(s_j\right)}+b_{s_j}{\mu }_{q\left(1/s_j\right)}\right)\right)\hfill \\ \hfill & +\frac{k}{2}{\mu }_{q(log\widehat{\sigma })}-c_{q}log\left(d_q\right)+\left(c_q-c\right){\mu }_{q\left({\eta }^2\right)}+log\left(\frac{\tau \left(c_q\right)}{\tau \left(c\right)}\right)+\left(d_q-d\right){\mu }_{q\left({\eta }^2\right)}.\hfill \end{array}$

The specific algorithm is shown in Algorithm 3.

4.2. Response with Nonignorable Missing Variable Choices

Only when the response variable has a nonignorable presence of missing values, similar to Section 3.3, we consider the missing data mechanism as a logistic regression model with parameters s and $\eta$. The prior distribution for these parameters follows the same approach as described in Section 4.1. The interactions between the parameters shown in Figure 4.

Algorithm 3: VI with Nonignorable Missing Covariates in Bayesian lasso regularized QR Models.

Input: The prior distribution of each parameter, including ${\sigma }_{\beta }^2,{\sigma }_{\varphi }^2,{\mu }_x,{\sigma }_x^2,{\sigma }_{{\mu }_x}^2,a^{*},b,c,d$, and a given lower bound on the evidence lower bound $tol={10}^{-5}$

For the above parameters, we consider the family of mean field variables.

$q\left(\beta ,s,{\eta }^2,z,\sigma ,y_{mis},\varphi ,\omega \right)=q\left(\beta \right)q\left(s\right)q\left({\eta }^2\right)q\left(z\right)q\left(\sigma \right)q\left(y_{mis}\right)q\left(\varphi \right)q\left(\omega \right)$

Since $s=\left(s_1,\dots ,s_k\right),s_j$ are independent of each other, so $q\left(s\right)=q\left(s_1\right)\cdots q\left(s_k\right)$. In addition, the optimal density functions for the parameters $s,\eta$ are given in Section 4.1, the optimal density functions for the parameters $z,\sigma ,y_{mis},\varphi ,\omega$ are given in Section 3.3, and only the optimal density functions for the parameter $\beta$ are given here as follows.

From Equation (7), we know that the parameter ${\beta }_j$ optimal density function is normally distributed

(45)$\begin{array}{c}\hfill q\left({\beta }_j\right)\sim N\left({\mu }_{q\left({\beta }_j\right)},{\sigma }_{q\left({\beta }_j\right)}^2\right)\end{array}$

where

$\begin{array}{cc}& {\sigma }_{q\left({\beta }_j\right)}^2=1/\left(\sum _{i=1}^n\frac{{\mu }_{q(1/\sigma )}{\mu }_{q(1/z_i)}{X}_j^2}{{\theta }^2}+{\mu }_{q(1/z_j)}\right)\hfill \\ & {\mu }_{q\left({\beta }_j\right)}={\sigma }_{q\left({\beta }_j\right)}^2\left[\frac{{\mu }_{q(1/\sigma )}{X}_j}{{\theta }^2}\left(E_{q\left(y_{mis}\right)}\left(y\right)\odot {\mu }_{q(1/z)}-\phi \mathrm{l}-\left(X{\mu }_{q\left(\beta \right)}-X_j{\mu }_{q\left({\beta }_j\right)}\right)\odot {\mu }_{q(1/z)}\right)\right].\hfill \end{array}$

In the model, the evidence lower bound by the evidence lower bound in this model is as follows, and the specific algorithm is described in Algorithm 4.

Algorithm 4: VI with Nonignorable Missing Response in Bayesian Lasso-Regularized QR Models.

Input: The prior distribution of each parameter, including ${\sigma }_{\beta }^2,{\sigma }_{\varphi }^2,{\mu }_x,{\sigma }_x^2,a^{*},b,c,d$, and a given lower bound on the evidence lower bound $tol={10}^{-5}$

(46)$\begin{array}{cc}\hfill & L\left(q\right)=\frac{1}{2}log\frac{\left|{\sum }_{q\left(\beta \right)}\right|}{{\sigma }_{\beta }^2}-\frac{1}{2{\sigma }_{\beta }^2}\left\{{∥{\mu }_{q\left(\beta \right)}∥}^2+tr\left({\mathsf{\Sigma }}_{q\left(\beta \right)}\right)\right\}-n{\mu }_{q(og(\sigma \left)\right)}-\frac{1}{2}\sum _{j=1}^{{\eta }_{m\epsilon }}log{\sigma }_{q\left(y_{m\in \dot{j},}\right)}^2\hfill \\ \hfill & -\sum _{i=1}^n\left({\mu }_{q\left(z_i\right)}{\mu }_{q(1/\sigma )}+log\frac{{\left(a_{z_i}/b_{z_i}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{z_i}{b}_{z_i}}\right)}-\frac{1}{2}{\mu }_{q\left(log{z}_i\right)}-\frac{1}{2}\left(a_{z_i}{\mu }_{q\left(z_i\right)}+b_{z_i}{\mu }_{q\left(1/z_i\right)}\right)\right)\hfill \\ \hfill & +\left(a_{q\left(\sigma \right)}-a\right){\mu }_{q(log(\sigma \left)\right)}+\frac{\left(b_{q\left(\sigma \right)}-b\right)a_{q\left(\sigma \right)}}{b_{q\left(\sigma \right)}}-a_{q\left(\sigma \right)}log{b}_{q\left(\sigma \right)}\frac{1}{2}log\left|{\sigma }_{q\left(\varphi \right)}^{2}I\right|-\frac{1}{2}{\mu }_{q\left(\varphi \right)}{}^T{\left(\sum _{\varphi }\right)}^{-1}{\mu }_{q\left(\mathsf{\Phi }\right)}\hfill \\ \hfill & +\sum _{i=1}^n\left[\left(R_i-\frac{1}{2}\right)E_q\left(Y_i\right){\mu }_{q\left(\varphi \right)}-log\left(\frac{exp\left({\xi }_i\right)}{1+exp\left({\xi }_i\right)}\right)-\frac{1}{2}{\xi }_i\right]-\frac{1}{2}tr\left({\left({\mathsf{\Sigma }}_{\varphi }\right)}^{-1}{\mathsf{\Sigma }}_{q\left(\varphi \right)}^2\right)\hfill \\ \hfill & -\frac{{\mu }_{q(1/\sigma )}}{2}\sum _{i=1}^n\left[E_q\left(y_i^2\right){\mu }_{q\left(1/z_i\right)}/{\theta }^2-\frac{2}{{\theta }^2}{E}_q\left(y_i\right)X_i^T{\mu }_{q\left(\beta \right)}{\mu }_{q\left(1/z_i\right)}-2E_q\left(y_i\right)\phi /{\theta }^2\right]+\left(d_q-d\right){\mu }_{q\left({\eta }^2\right)}\hfill \\ \hfill & -\frac{{\mu }_{q(1/\sigma )}}{2}\sum _{i=1}^n\left[\frac{\phi }{{\theta }^2}{X}_i^T{\mu }_{q\left(\beta \right)}+\frac{{\mu }_{q\left(1/z_i\right)}}{{\theta }^2}tr\left(X_i^T{X}_i\left({\mathsf{\Sigma }}_{q\left(\beta \right)}+{\mu }_{q\left(\beta \right)}{\mu }_{q\left(\beta \right)}^T\right)\right)+\frac{{\phi }^2}{{\theta }^2}{\mu }_{q\left(\left({}_i\right)\right.}\right]\hfill \\ \hfill & -\frac{1}{2}\sum _{i=1}^n{\mu }_{q\left(log\left(z_i\right)\right.}-\frac{n}{2}{\mu }_{q(log(\sigma \left)\right)}+\frac{k}{2}{\mu }_{q\left(\sigma \right)}(log\sigma )-c_{q}log\left(d_q\right)+\left(c_q-c\right){\mu }_{q\left({\eta }^2\right)}+log\left(\frac{\tau \left(c_q\right)}{\tau \left(c\right)}\right)\hfill \\ \hfill & -\sum _{i=1}^k\left(-\frac{1}{2}{\mu }_{q\left(s_i\right)}{\mu }_{q\left({\eta }^2\right)}-log\frac{{\left(a_{s_i}/b_{s_i}\right)}^{1/4}}{2{\mathrm{K}}_{1/2}\left(\sqrt{a_{s_i}{b}_{s_i}}\right)}+\frac{1}{2}{\mu }_{q\left(s_i\right)}\left(log{s}_i\right)+\frac{1}{2}\left(a_{s_i}{\mu }_{q\left(s_i\right)}+b_{s_i}{\mu }_{q\left(1/s_i\right)}\right)\right)\hfill \end{array}$

5. Simulation Studies

In this section, we conducted four experiments to validate the effectiveness of variational inference for parameter estimation in a quantile regression model with non-negligible missing data. The computer involved used the Windows operating system with an Intel® Core i5-8400 six-core processor. The experimental data were generated based on the following quantile regression models:

$\begin{array}{c}{y}_i=X_i^T\beta +{\epsilon }_i,i=1,\dots ,n\\ {\epsilon }_i\sim ALD(0,\sigma ,p)\end{array}$

In the above equation, $X_i$$=\left(x_{i1},\dots ,x_{ik}\right)$ and $\beta =\left({\beta }_0,\dots ,{\beta }_k\right)$, where $x_{i1}=1$. Furthermore, we also performed Gibbs sampling simulations for different scenarios in Simulations 1 and II. We compared the results in terms of CPU time, the deviation of parameter estimates from the true values (BIAS, $\left|\frac{1}{100}{\sum }_{i=1}^{100}\widehat{\beta }-\beta \right|$), the root mean square error (RMSE, $\sqrt{\frac{1}{100}{\sum }_{i=1}^{100}{(\widehat{\beta }-\beta )}^2}$), and the standard deviation (SD). For Simulations I and II, we introduced approximately 10% missing data as M1 and approximately 20% missing data as M2 under different error terms.

• Simulation I:

In this simulation, we take $n=200$ or 30, $k=2,x_{i2}$ randomly generated from the normal distribution $N(0,0.6)$ and consider the case of missing covariates $x_{i2}$. When $n=200$, we take the parameters $\mathit{\beta }=\left({\beta }_0,{\beta }_1\right)$, and the truth values of $\sigma$ are ${\beta }_0={\beta }_1=1,\sigma =1$.

The prior distributions for the parameters $\beta$ and $\sigma$ are set as $\mathit{\beta }\sim N(0,I)$ and $\sigma \sim IG(1,0.0009)$. We consider different distributions for ${\epsilon }_i$ and examine various deletion cases at different quartiles within each distribution. The details of these cases are as follows: (C1): ${\epsilon }_i$ follows a standard normal distribution $N(0,0.7)$. We set the true value of $\mathit{\varphi }$ as $(2.5,-2.7)$, and the missing rates of $x_{i2}$ are approximately $10.9\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(0.75,-1.05)$, the missing rate of $x_{i2}$ is about $23\%$. (C2): ${\epsilon }_i$ follows a mixed normal distribution $0.9N(0,0.7)+0.1N(0,1)$. The true value of $\mathit{\varphi }$ is set as $(2.4,-2.8)$, and the missing rates of $x_{i2}$ are approximately $11\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(0.75,-1)$, the missing rate of $x_{i2}$ is about $22\%$. (C3): ${\epsilon }_i$ follows a chi-squared distribution $0.3{\chi }_5^2$. We set the true value of $\mathit{\varphi }$ as $(2.55,-2.75)$, and the missing rates of $x_{i2}$ are approximately $9\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(0.75,-1)$, the missing rate of $x_{i2}$ is about $23\%$. When $n=30$, We examined the M1 situation, maintaining the parameters and error distributions as previously detailed.

• Simulation II:

In this simulation, we consider a scenario where $n=200$ or 30, $k=2$. The covariate $x_{i2}$ is randomly generated from a normal distribution $N(0,0.6)$. When n = 200, we focus on the case of missing response variable $y_i$. The true values of the parameters $\mathit{\beta }=\left({\beta }_0,{\beta }_1\right)$ and $\sigma$ are set as follows: $n{\beta }_0={\beta }_1=1,\sigma =1.$ The prior distributions for the parameters $\beta$ and $\sigma$ are set as $\mathit{\beta }\sim N(0,0.8I)$ and $\sigma \sim IG(1,0.05)$. We consider different distributions for ${\epsilon }_i$ and examine various deletion cases at different quartiles within each distribution. The details of these cases are as follows: (C1): ${\epsilon }_i$ follows a standard normal distribution $N(0,0.3)$. We set the true value of $\mathit{\varphi }$ as $(1.5,1.5)$, and the missing rates of $x_{i2}$ are approximately $11\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(0.82,0.85)$, the missing rates of $x_{i2}$ are approximately $19\%$. (C2): ${\epsilon }_i$ follows a mixed normal distribution $0.9N(0,0.1)+0.1N(0,0.3)$. We set the true value of $\mathit{\varphi }$ as $(2.1,0.29)$, and the missing rates of $x_{i2}$ are approximately $12\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(1.5,0.3)$, the missing rates of $x_{i2}$ are approximately $23\%$. (C3): ${\epsilon }_i$ follows a chi-square distribution $0.1{\chi }_3^2$. We set the true value of $\mathit{\varphi }$ as $(2,0.45)$, and the missing rates of $x_{i2}$ are approximately $9\%$. Alternatively, when we set the true value of $\mathit{\varphi }$ as $(2,-0.15)$, the missing rates of $x_{i2}$ are approximately $19\%$. When $n=30$, we examined the M1 situation, maintaining the parameters and error distributions as previously detailed.

By comparing the results with and without the inclusion of the scale parameter $\sigma$ in Simulation I and Simulation II, when $n=200$, we observed the following findings for different missing cases and error distributions. The parameter estimation results obtained using the VI algorithm for Simulation I and Simulation II are presented in Table 1 and Table 2, respectively. Meanwhile, the parameter estimation results obtained using Gibbs are displayed in Table A1 and Table A2 (Appendix B). The Bayesian estimates obtained by the proposed method we present perform satisfactorily in that (i) they have smaller BIAS, RMS, and SD, and the SD and RMS values are close to each other; and (ii) the comparison between Table 1 and Table A1, and Table 2 and Table A2 reveals that different error distributions, various missing cases, and the inclusion of the parameter $\sigma$ have a significant impact on CPU time. When considering covariates with nonignorable missingness and including the parameter $\sigma$, VI performs over ten times faster than Gibbs. Without including $\sigma$, VI is nearly 50 times faster than Gibbs. Similarly, in the presence of nonignorable missing response variables and including sigma, VI is nearly 50 times faster than Gibbs. Without $\sigma$, VI speeds up by almost 80 times compared to Gibbs. This means our proposed variational Bayesian method is faster than MCMC while ensuring parameter estimation accuracy. When $n=30$, although VI is almost a hundred times faster than Gibbs, its accuracy is somewhat lower. This difference is even more noticeable when the response variable is absent, for specific displays in Table 3.

Based on the good performance of variational inference observed in Simulation I and Simulation II, we employed Bayesian lasso for variable selection in the models with missing covariates and response variables in Simulations III and IV, respectively. Furthermore, the variational Bayesian inference achieved convergence in both Simulations I and II. For convenience, we provide the convergence results of C1 at different quartiles of M1 when the covariates have missing values and the response variables have missing values, as depicted in Figure 5, see Figure 6 for run time.

• Simulation III:

In this simulation, we utilize a Bayesian lasso quantile regression model to perform variable selection and parameter estimation in the presence of missing covariates. The data generation process follows the equation $y_i={\beta }_0+X_i^T\beta +{\epsilon }_i,i=1,\dots ,n$, where we omit the scale parameter $\sigma$. The true values are set as ${\beta }_0=1$, $\mathit{\beta }={\left({\beta }_1,\dots ,{\beta }_7\right)}^T$, $X_i=\left(x_{i1},\dots ,x_{i7}\right)$, with $n=200$. The covariates $x_{i1},\dots ,x_{i7}$ are independently sampled from a normal distribution $N\left(0,0.6\right)$, and we take $\beta =(3,1.5,0,0,2,0,0)$ with a prior distribution of $N(0,1.5I)$. The error term ${\epsilon }_i$ is considered under three different distributions, similar to Simulations I and II, when $x_{i1}$ is missing. (C1), the error term ${\epsilon }_i$ follows a standard normal distribution $N\left(0,0.7\right)$. By setting the true values as $\mathit{\varphi }=(0.70,-1.05)$, the missing rates of $x_{i1}$ are approximately $28\%$. (C2), the error term ${\epsilon }_i$ follows a mixed normal distribution $0.9N\left(0,0.7\right)+0.1N\left(0,1\right)$. Taking the true values as $\mathit{\varphi }=(0.75,-1.00)$, the missing rates of $x_{i1}$ are approximately $25\%$. (C3), the error term ${\epsilon }_i$ follows a chi-square distribution $0.3{\chi }_5^2$. With the true values $\mathit{\varphi }=(0.75,-1.00)$, the missing rates of $x_{i1}$ are approximately $25\%$.

• Simulation IV:

In this simulation, we employ a Bayesian lasso-quantile regression model to handle variable selection and parameter estimation in the presence of missing response variables. The data generation process is similar to Simulations III and does not include the scale parameter $\sigma$. We set $\mathit{\beta }=(1,-1,1,0,0,0,0)$ and assume a prior distribution of $N(0,I)$. The error term ${\epsilon }_i$ is considered under three different distributions, following the pattern of Simulation I and Simulation II when $y_i$ is missing. (C1), the error term ${\epsilon }_i$ follows a standard normal distribution $N\left(0,0.3\right)$. With the true values set as $\mathit{\varphi }=(0.82,0.80)$, the missing rates of $y_i$ are approximately $23\%$. (C2), the error term ${\epsilon }_i$ follows a mixed normal distribution $0.9N\left(0,0.1\right)+0.1N\left(0,0.9\right)$. Taking the true values as $\mathit{\varphi }=(1.08,0.60)$, the missing rates of $y_i$ are approximately $16\%$. (C3), the error term ${\epsilon }_i$ follows a chi-square distribution $0.1{\chi }_3^2$. With the true values $\mathit{\varphi }=(1.00,0.45)$, the missing rates of $y_i$ re about $18\%$. In Simulation III and Simulation IV, we repeated 100 times and evaluated the performance of variable selection and parameter estimation using the following metrics: (1) the $L_2$ parametrization between the parameter estimates and the true value $\left(L_2={(\widehat{\beta }-\beta )}^T(\widehat{\beta }-\beta )\right)$; (2) mean square error (MSE) $\left(MSE={(\widehat{\beta }-\beta )}^T{X}^{T}X(\widehat{\beta }-\beta )/n\right)$; (3) the number of parameters that have a true value of zero and are correctly identified as zero is recorded as “C”; (4) the number of parameters whose true value is not zero but is incorrectly identified as zero is recorded as “IC”.

We investigated three interquartile levels (0.25, 0.5, 0.75), as shown in Table 4 and Table 5, yielding the following results:

• (1). The “IC” value is extremely close to zero, and “C” closely approximates the number of zeros in the real parameter.

• (2). The $L_2$ values for the three interquartile levels are below 0.03, indicating a negligible difference from zero.

• (3). The “MSE” value at all three quantile levels is less than 0.09.

These findings demonstrate the strong performance of our proposed VI for variable selection.

6. A Real Example

In this section, we demonstrate the application of Algorithm 3 using data obtained from the American Press Institute, specifically the 1995 American University News Report dataset. The dataset can be accessed at http://lib.stat.cmu.edu/datasets/colleges (accessed on 21 April 2023). Our focus is on examining the relationship between graduation rates (y) and various factors.

We select the following factors for analysis: the public/private indicator as $x_1$, the student-to-faculty ratio as $x_2$, the natural logarithm of the number of applicants accepted as $x_3$, the average Combined SAT score as $x_4$, the average ACT score as $x_5$, the natural logarithm of the number of applications received as $x_6$, the natural logarithm of the number of new students enrolled as $x_7$, room and board costs as $x_8$, and the natural logarithm of instructional expenditure per student as $x_9$.

To ensure a complete set of observations, we assume $x_1$, $x_2$, and y to be fully observed. However, we remove one erroneous case from the response variable y, resulting in a total of 1203 cases. The final missing rates for the remaining variables are as follows: $x_3$ (0.41%), $x_4$ (39.23%), $x_5$ (45.22%), $x_6$ (0.49%), $x_7$ (0.16%), $x_8$ (4.90%), and $x_9$ (2.07%).

Furthermore, to facilitate statistical inferences on the data, we apply two preprocessing steps. First, we take the natural logarithm of the variables $x_3$, $x_6$, $x_7$, $x_8$, and $x_9$. Second, we standardize all the variables to eliminate any scale differences among them. These preprocessing steps ensure that the data are appropriately transformed and normalized for subsequent analysis and inference.

We consider the following quantile regression model.

$Q\left(y_i\mid x_i,{\beta }_z\right)={\beta }_{z0}+x_i^T{\beta }_{\tau }$

where $i=1,\dots ,1203,x_i={\left(x_{i1},\dots ,x_{i9}\right)}^T,{\beta }_{\tau }={\left({\beta }_{\tau 1},\dots ,{\beta }_{\tau 9}\right)}^T,{\beta }_{\tau 0}$ denotes the intercept term. For $x_3,x_4,x_5,x_6,x_7,x_8,x_9$ with missing $x_3,x_4,x_5,x_6,x_7,x_8,x_9$, they are assumed to follow a normal distribution, and the expectation and variance are the mean and variance of the observable data, respectively. The missing data mechanism is set as in the algorithm. We consider the above missingness as non-negligible missingness. We calculate the algorithm’s success based on $\tau =0.5$ and the estimates of the last known parameters: (1) $x_2,x_3,x_6,x_7$ have a small effect on the response variable, and the 0.5 quantile regression can be written as $y=-0.0334+0.3947x_1+0.3185x_4+0.0649x_5+0.148x_8+0.034x_9$ (2) with missing $x_3,x_6,x_7$ missing mechanisms for MAR, $x_4,x_5,x_8,x_9$ for MNAR. Lower bound convergence in Figure 7, Bayesian estimates (EST) and 95% confidence intervals (CI) of the parameters in Table 6.

7. Discussion

In existing research on quantile regression models with significant missing data, most studies have utilized the Markov Chain Monte Carlo (MCMC) algorithm for Bayesian inference, despite some drawbacks such as sampling difficulties and extensive computational time. In this paper, we propose employing the variational Bayesian algorithm for statistical inference of the aforementioned model while incorporating variable selection. Specifically, in the presence of missing covariates, we consider a probit regression model for the missing data mechanism, and for the missing response variable, we consider a logistic regression model. Additionally, we employ lasso regularization for variable selection in both cases. Simulation and real-world experiments demonstrate that when covariates or response variables have nonignorable missing values in quantile regression models, the variational Bayesian method not only ensures inference precision but also consumes less time compared to MCMC.

Based on our experience, we often encounter situations where both covariates and response variables have missing values. The aforementioned developed variational Bayesian approach faces a well-known ill-posed problem in practical applications. To address this issue, we explore variational Bayesian parameter estimation for quantile regression models in the presence of simultaneous missing covariates or response variables.

Our proposed method has some potential drawbacks. For instance, as the number of dimensions increases, the convergence speed of variational Bayes may slow down. Furthermore, the calculation of expectations may lack an analytical solution, and assumptions about the correctness of the missing data mechanism must be considered. To solve these issues, we can employ commonly used machine learning techniques like deep learning and neural networks to address missing data mechanisms. Furthermore, in the future, the stochastic variational Bayesian method can enhance the algorithm’s efficiency, while the Bayesian dimensionality reduction approach can tackle issues related to high-dimensional quantile regression problems.

Author Contributions

Conceptualization, M.T.; methodology, X.L. and M.T.; software, X.L.; data curation, X.L.; writing—original draft, X.L.; supervision, M.T. and X.H.; funding acquisition, X.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Not applicable.

Acknowledgments

We are grateful to the editor, associated editor, and four referees for their valuable suggestions and comments that greatly improved the article, grateful to both StatLib and the original contributor of the material.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. DAG for the missing covariates for Quantile regression model.

Figure 2. DAG for the missing response for Quantile regression model.

Figure 3. DAG for the missing covariates for Bayesian lasso regularized quantile regression model.

Figure 4. DAG for the missing response for Bayesian lasso regularized quantile regression model.

Figure 5. Variable fractional Bayesian convergence effect.

Figure 6. Comparison of CPU time of variational Bayesian over Gibbs when covariates with missing.

Figure 7. Lower bound convergence of evidence in the case analysis.

Table 1

Parameter estimation and CPU time cases in Simulation I using VI, $n=200$.

	Method	Missing	p	BIAS		RMSE		SD		CPU Time
				${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$
without	C1	M1	0.25	0.006	0.078	0.062	0.138	0.061	0.114	124.4
			0.5	0.049	0.040	0.079	0.109	0.062	0.101	111
			0.75	0.104	0.0001	0.126	0.120	0.071	0.121	133.2
		M2	0.25	0.014	0.082	0.069	0.152	0.068	0.128	185
			0.5	0.071	0.022	0.095	0.117	0.063	0.116	176.6
			0.75	0.129	0.001	0.147	0.134	0.072	0.136	201.8
	C2	M1	0.25	0.066	0.076	0.088	0.132	0.059	0.109	138.4
			0.5	0.056	0.039	0.076	0.101	0.050	0.095	139.8
			0.75	0.036	0.030	0.079	0.105	0.072	0.101	137.4
		M2	0.25	0.068	0.077	0.088	0.127	0.055	0.101	178.2
			0.5	0.063	0.050	0.084	0.123	0.053	0.109	165
			0.75	0.053	0.039	0.083	0.120	0.064	0.113	236.6
	C3	M1	0.25	0.028	0.040	0.058	0.101	0.053	0.096	124.8
			0.5	0.081	0.015	0.142	0.111	0.073	0.126	139.6
			0.75	0.121	0.010	0.136	0.131	0.117	0.208	136.6
		M2	0.25	0.045	0.019	0.121	0.130	0.062	0.103	172.2
			0.5	0.110	0.018	0.124	0.109	0.072	0.153	158.8
			0.75	0.137	0.018	0.105	0.162	0.108	0.109	268.4
with	C1	M1	0.25	0.025	0.009	0.062	0.074	0.058	0.074	387.83
			0.5	0.048	0.019	0.066	0.083	0.045	0.081	342.17
			0.75	0.035	0.011	0.054	0.088	0.040	0.088	359.41
		M2	0.25	0.045	0.011	0.062	0.085	0.043	0.085	1068.91
			0.5	0.088	0.025	0.097	0.082	0.041	0.079	963.02
			0.75	0.125	0.015	0.133	0.093	0.046	0.093	955.78
	C2	M1	0.25	0.098	0.009	0.105	0.074	0.038	0.074	384.02
			0.5	0.059	0.004	0.070	0.063	0.039	0.064	335.35
			0.75	0.039	0.004	0.057	0.072	0.041	0.073	385.50
		M2	0.25	0.098	0.016	0.105	0.075	0.040	0.073	1019.75
			0.5	0.078	0.008	0.088	0.076	0.039	0.076	885.35
			0.75	0.069	0.013	0.084	0.083	0.048	0.083	925.44
	C3	M1	0.25	0.026	0.004	0.035	0.046	0.023	0.046	320.22
			0.5	0.072	0.004	0.085	0.075	0.046	0.075	318.81
			0.75	0.151	0.039	0.162	0.109	0.059	0.102	414.64
		M2	0.25	0.045	0.006	0.054	0.050	0.029	0.050	903.49
			0.5	0.117	0.037	0.125	0.086	0.042	0.078	865.73
			0.75	0.175	0.036	0.186	0.130	0.063	0.125	1109.56

Table 2

Parameter estimation and CPU time cases in Simulation II using VI, $n=200$.

	Method	Missing	p	BIAS		RMSE		SD		CPU Time
				${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$
with	C1	M1	0.25	0.002	0.051	0.042	0.095	0.042	0.081	153.4
			0.5	0.063	0.019	0.074	0.072	0.038	0.070	132.2
			0.75	0.128	0.058	0.132	0.098	0.030	0.079	152.6
		M2	0.25	0.011	0.025	0.035	0.086	0.033	0.082	211.2
			0.5	0.076	0.025	0.085	0.088	0.038	0.084	190.8
			0.75	0.151	0.019	0.156	0.079	0.042	0.077	200
	C2	M1	0.25	0.015	0.003	0.026	0.043	0.021	0.043	115.8
			0.5	0.003	0.001	0.020	0.038	0.020	0.039	103.6
			0.75	0.077	0.009	0.079	0.057	0.020	0.057	148.4
		M2	0.25	0.038	0.022	0.045	0.056	0.023	0.051	168.4
			0.5	0.048	0.032	0.040	0.030	0.019	0.030	157.2
			0.75	0.117	0.030	0.012	0.064	0.029	0.055	173
	C3	M1	0.25	0.015	0.003	0.026	0.043	0.021	0.043	127.8
			0.5	0.003	0.001	0.020	0.038	0.020	0.039	84
			0.75	0.077	0.009	0.079	0.057	0.020	0.057	134.6
		M2	0.25	0.043	0.022	0.053	0.055	0.030	0.051	150
			0.5	0.048	0.032	0.054	0.057	0.026	0.048	134
			0.75	0.117	0.032	0.121	0.064	0.029	0.055	214
without	C1	M1	0.25	0.024	0.060	0.038	0.102	0.031	0.087	69.4
			0.5	0.066	0.055	0.078	0.078	0.031	0.055	58.8
			0.75	0.057	0.061	0.065	0.088	0.030	0.064	96
		M2	0.25	0.049	0.090	0.128	0.116	0.064	0.076	88.4
			0.5	0.037	0.087	0.040	0.100	0.030	0.051	77.6
			0.75	0.049	0.018	0.057	0.185	0.028	0.046	113.5
	C2	M1	0.25	0.027	0.012	0.028	0.045	0.008	0.046	65
			0.5	0.019	0.019	0.025	0.031	0.015	0.025	57.8
			0.75	0.002	0.013	0.017	0.037	0.017	0.035	81.6
		M2	0.25	0.102	0.051	0.104	0.070	0.022	0.048	105.17
			0.5	0.078	0.050	0.081	0069	0.021	0.047	97.80
			0.75	0.057	0.044	0.061	0.070	0.023	0.054	120.80
	C3	M1	0.25	0.010	0.036	0.023	0.072	0.022	0.066	83.8
			0.5	0.009	0.015	0.024	0.040	0.023	0.038	86.4
			0.75	0.076	0.027	0.080	0.059	0.025	0.052	103.2
		M2	0.25	0.033	0.011	0.041	0.070	0.026	0.074	108
			0.5	0.052	0.042	0.057	0.066	0.024	0.051	112.2
			0.75	0.122	0.055	0.012	0.083	0.026	0.062	150.6

Table 3

Parameter estimation and CPU time cases in Simulation I and Simulation II, $n=30$.

	Method	Missing	p	BIAS		RMSE		SD		CPU Time
				${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$
Simulation I	C1	VI	0.25	0.199	0.412	0.244	0.481	0.143	0.262	18.4
			0.5	0.032	0.250	0.148	0.344	0.149	0.232	22.2
			0.75	0.033	0.264	0.111	0.346	0.108	0.226	25.6
		Gibbs	0.25	0.033	0.065	0.123	0.267	0.128	0.264	3353.2
			0.5	0.022	0.070	0.104	0.261	0.108	0.253	3347.8
			0.75	0.035	0.141	0.132	0.268	0.105	0.223	4275
	C2	VI	0.25	0.086	0.333	0.145	0.116	0.113	0.256	15.8
			0.5	0.065	0.281	0.11	0.365	0.090	0.239	17.6
			0.75	0.007	0.299	0.099	0.357	0.093	0.257	18.4
		Gibbs	0.25	0.038	0.062	0.125	0.266	0.123	0.251	3668.4
			0.5	0.018	0.132	0.100	0.230	0.119	0.230	3357.2
			0.75	0.117	0.112	0.192	0.264	0.169	0.235	3873
	C3	VI	0.25	0.145	0.433	0.186	0.493	0.121	0.243	16.8
			0.5	0.056	0.271	0.120	0.358	0.090	0.229	18
			0.75	0.037	0.239	0.149	0.327	0.130	0.237	17.6
		Gibbs	0.25	0.063	0.082	0.133	0.225	0.110	0.211	4150
			0.5	0.158	0.172	0.234	0.357	0.176	0.318	3834
			0.75	0.127	0.192	0.361	0.464	0.359	0.425	3614
Simulation II	C1	VI	0.25	0.124	0.350	0.178	0.442	0.121	0.277	14.7
			0.5	0.046	0.275	0.138	0.418	0.131	0.315	16.7
			0.75	0.005	0.021	0.145	0.348	0.140	0.284	17
		Gibbs	0.25	0.030	0.040	0.088	0.116	0.074	0.126	1450.4
			0.5	0.007	0.017	0.050	0.090	0.060	0.091	1473.6
			0.75	0.007	0.030	0.077	0.145	0.078	0.146	1539.5
	C2	VI	0.25	0.167	0.312	0.198	0.425	0.118	0.286	15.2
			0.5	0.009	0.269	0.155	0.341	0.155	0.225	15.8
			0.75	0.042	0.273	0.187	0.372	0.187	0.275	16.2
		Gibbs	0.25	0.062	0.006	0.084	0.050	0.082	0.048	1566.17
			0.5	0.028	0.005	0.004	0.059	0.031	0.057	1600.80
			0.75	0.003	0.024	0.021	0.060	0.023	0.054	1666.80
	C3	VI	0.25	0.090	0.246	0.123	0.312	0.092	0.266	15.4
			0.5	0.059	0.135	0.193	0.290	0.193	0.268	16.3
			0.75	0.009	0.272	0.240	0.499	0.225	0.422	15.9
		Gibbs	0.25	0.033	0.001	0.041	0.060	0.046	0.064	1600
			0.5	0.012	0.008	0.057	0.086	0.044	0.081	1600.2
			0.75	0.012	0.025	0.062	0.003	0.066	0.082	1500.6

Table 4

Results of variable selection for Simulation III.

	0.25				0.5				0.75
	L2	MSE	C	IC	L2	MSE	C	IC	L2	MSE	C	IC
C1	0.0276	0.0072	4.85	0	0.0056	0.0022	4.94	0	0.0251	0.0080	4.94	0
C2	0.0252	0.0080	4.89	0	0.0063	0.0062	4.93	0	0.0184	0.0074	4.88	0
C3	0.0134	0.0048	4.86	0	0.0048	0.0019	4.75	0	0.0124	0.0045	4.9	0

Table 5

Results of variable selection for Simulation IV.

	0.25				0.5				0.75
	L2	MSE	C	IC	L2	MSE	C	IC	L2	MSE	C	IC
C1	0.0369	0.0118	4.93	0	0.0035	0.0011	4.93	0	0.0170	0.0053	4.95	0
C2	0.0071	0.0021	4.95	0	0.0036	0.0012	4.89	0	0.0063	0.0020	4.93	0
C3	0.0174	0.0053	4.97	0	0.0044	0.0016	4.91	0	0.0054	0.0018	4.95	0

Table 6

Bayesian estimates (EST) and 95% confidence intervals (CI) of the parameters in real example.

Par	Est (CI)	Par	Est (CI)
${\beta }_0$	−0.0334 (−0.0369,−0.0300)	${\varphi }_{4,0}$	0.2552 (0.2538,0.2567)
${\beta }_1$	0.3947 (0.3885,0.4009)	${\varphi }_{4,1}$	−0.1527 (−0.1622,0.1432)
${\beta }_2$	−0.0500 (−0.1102,0.0105)	${\varphi }_{5,0}$	0.1066 (0.1066,0.1067)
${\beta }_3$	0.0121 (−0.1023,0.1265)	${\varphi }_{5,1}$	−0.0284 (−0.0307,−0.0261)
${\beta }_4$	0.3185 (0.3114,0.3256)	${\varphi }_{6,0}$	1.5119 (1.5117,1.5121)
${\beta }_5$	0.0649 (0.0616,0.0682)	${\varphi }_{6,1}$	0.0036 (0.0028,0.0044)
${\beta }_6$	0.0335 (−0.0428,0.1138)	${\varphi }_{7,0}$	1.5435 (1.5344,1.5346)
${\beta }_7$	−0.0600 (−0.1299,0.0098)	${\varphi }_{7,1}$	0.0002 (0,0.0003)
${\beta }_8$	0.1480 (0.1387,0.1573)	${\varphi }_{8,0}$	1.2650 (1.2648,1.2652)
${\beta }_9$	0.0340 (0.0208,0.0472)	${\varphi }_{8,1}$	−0.0416 (−0.0442,−0.0399)
${\varphi }_{3,0}$	1.5177 (1.5174,1.5180)	${\varphi }_{9,0}$	1.4132 (1.4130,1.4133)
${\varphi }_{3,1}$	0.0024 (0.0014,0.0034)	${\varphi }_{9,1}$	−0.0042 (−0.0061,−0.0024)

Appendix A

The derivations of Algorithms 1–4 are similar to each other; thus, we only refine the steps of Algorithm 1. The full conditional probability posterior are as follows:(A1)$$\begin{array}{cc}\hfill & p(\beta \mid {\mathsf{\Theta }}_r)\propto exp\left\{-\frac{1}{2}\left[{\beta }^T\left(\frac{1}{\sigma }\sum _{i=1}^n\frac{X_i{X}_i^T}{{\theta }^2{z}_i}+\frac{1}{{\sigma }_{\beta }^2}I\right)\beta -2{\beta }^T\left(\frac{1}{\sigma }\sum _{i=1}^n\frac{\left(y_i-\phi z_i\right)x_i}{{\theta }^2{z}_i}\right)\right]\right\}\hfill \end{array}$$(A2)$$\begin{array}{cc}\hfill & p(z\mid {\mathsf{\Theta }}_r)\propto \prod _{i=1}^n{z}_i^{-\frac{1}{2}}exp\left\{-\frac{1}{2}\left[\frac{{\left(y_i-X_i^T\beta \right)}^2}{{\theta }^2\sigma z_i}+\frac{2{\theta }^2+{\phi }^2}{{\theta }^2\sigma }{z}_i\right]\right\}\hfill \end{array}$$(A3)$$\begin{array}{cc}\hfill & p(\sigma \mid {\mathsf{\Theta }}_r)\propto {\left(\frac{1}{\sigma }\right)}^{\frac{3n}{2}+a+1}exp\left\{-\frac{1}{\sigma }\left(\frac{\left(y-X^T\beta -\phi z_i\right){\sum }_1^{-1}\left(y-X^T\beta -\phi z_i\right)}{2}+\sum _{i=1}^n{z}_i+b\right)\right\}:\hfill \end{array}$$(A4)$$\begin{array}{cc}\hfill & p\left(x_{mis}\right|{\mathsf{\Theta }}_r)\propto \prod _{i=1}^{Nmis}\frac{exp\left[-\frac{1}{2}{x}_i(\frac{{\beta }_i^2}{{\theta }^2\sigma z_i}+\frac{1}{{\sigma }_x^2}+{\varphi }_1^2)x_i\right]}{exp\left[-x_i(\frac{{\mu }_x}{{\sigma }_x^2}-\frac{{\beta }_0{\beta }_1}{{\theta }^2\sigma z_i}-{\varphi }_0{\varphi }_1+a_i{\varphi }_1+\frac{{\beta }_1{y}_i}{{\theta }^2\sigma z_i}-\frac{{\beta }_1\phi }{2{\theta }^2\sigma })\right]}\hfill \end{array}$$(A5)$$\begin{array}{cc}\hfill & p\left({\mu }_x\mid {\mathsf{\Theta }}_r\right)\propto exp\left\{-\frac{1}{2}\left[{\mu }_x\left(\frac{n}{{\sigma }_x^2}+\frac{1}{{\sigma }_{{\mu }_x}^2}\right){\mu }_x-2{\mu }_x\frac{{\sum }_i^{N_{obs}}{x}_{obs,i}}{{\sigma }_{{\mu }_x}^2}\right]\right\}\hfill \end{array}$$(A6)$$\begin{array}{cc}\hfill & p\left({\sigma }_x^2\mid {\mathsf{\Theta }}_r\right)\sim IG\left(A_x+\frac{1}{2}n,B_x+\frac{1}{2}{∥X_j-{\mu }_{x}1∥}^2\right)\hfill \end{array}$$(A7)$$\begin{array}{cc}\hfill & p(\varphi \mid {\mathsf{\Theta }}_r)\propto exp\left\{-\frac{1}{2}\left[{\varphi }^T\left(x^{T}x+\frac{1}{{\sigma }_{\varphi }^2}\right)\varphi -2{\varphi }^{T}xa\right]\right\}\hfill \end{array}$$(A8)$$\begin{array}{cc}\hfill & p(a\mid {\mathsf{\Theta }}_r)={\left(2\pi \right)}^{-\frac{n}{2}}\prod _{i=1}^n\left[{\left\{\frac{I\left(a_i>0\right)}{\mathsf{\Phi }\left(x_i\varphi \right)}\right\}}^{R_i}{\left\{\frac{I\left(a_i<0\right)}{1-\mathsf{\Phi }\left(x_i\varphi \right)}\right\}}^{1-R_i}exp\left\{-\frac{1}{2}{\left(a_i-{\varphi }_0-{\varphi }_1{x}_i\right)}^2\right\}\right]\hfill \end{array}$$

In Step 4 of Algorithm 1, $$\begin{array}{cc}& E_{q\left(x_{mis}\right)q\left(\beta \right)\phi \left(z\right)}\sum _{i=1}^n\left[\frac{{\left(y_i-X_i^T\beta -\phi z_i\right)}^2}{{\theta }^2{z}_i}\right]\hfill \\ & =\frac{1}{{\theta }^2}\sum _{i=1}^n\left(y_i^2{\mu }_{q(1/z_i)}-2y_i{E}_{q\left(x_{\mathrm{mis}}\right)}\left(X_i\right){\mu }_{q\left(\beta \right)}{\mu }_{q(1/z_i)}\right)\hfill \\ & -\frac{1}{{\theta }^2}\sum _{i=1}^n\left(-2y_i\phi +2E_{q\left(x_{mis}\right)}\left(X_i\right){\mu }_{q\left(\beta \right)}\phi +{\mu }_{q(1/z_i)}{\phi }^2+{\mu }_{q(1/z_i)}{E}_q{\left(X_i^T\beta \right)}^2\right)\hfill \end{array}$$ and, $$E_q{\left(X_i^T\beta \right)}^2=tr\left[E_{q\left(x_{mis}\right)}\left(X_i^T{X}_i\right)({\sigma }_{q\left(\beta \right)}^2+{\mu }_{q\left(\beta \right)}^T{\mu }_{q\left(\beta \right)})\right]$$

In Step 7 of Algorithm 1, $$E_q{∥x-{\mu }_{x}1∥}^2={∥x_{obs}-{\mu }_{q\left({\mu }_x\right)}1∥}^2+{∥{\mu }_{q\left(x_{mis}\right)}-{\mu }_{q\left({\mu }_x\right)}1∥}^2+n{\sigma }_{q\left({\mu }_x\right)}^2+\sum _{j=1}^{n_{mis}}{\sigma }_{q\left(x_{mis,j}\right)}^2$$

Expressions for $E_{q\left(x_{mis}\right)}\left(X\right),E_{q\left(x_{mis}\right)}\left(X_i^T{X}_i\right)$, $E_q\left(X\right)=\left(\begin{array}{cc}1\hfill & E_q\left(x\right)\hfill \end{array}\right)$, where $$E_q\left(x\right)=\left(\begin{array}{c}{x}_{obs}\hfill \\ {\mu }_{q\left(x_{\mathrm{mis}}\right)}\hfill \end{array}\right)$$

In other words, by replacing $x_i$ with ${\mu }_{q\left(x_{mis,i}\right)}$ in X with the missing $x_i$, we obtain $E_q\left(X\right)$. Since $X_i=\left(1x_i\right)$ is a vector of $1\times 2$, we have, $$E_q\left(X_i{X}_i^T\right)=\left(\begin{array}{cc}1& a_{12}\\ a_{21}& a_{22}\end{array}\right)$$

It is known that $E_q\left(X_i{X}_i^T\right)$ is a symmetric matrix. When $x_i$ is missing, $a_{12}={\mu }_{q\left(x_{mis,i}\right)},$ $a_{22}={\left({\mu }_{q\left(x_{mis,i}\right)}\right)}^2+{\sigma }_{q\left(x_{mis,i}\right)}^2$. then, $$E_q\left(X{X}^T\right)=\left(\begin{array}{cc}n& 1^T{x}_{obs}+1^T{\mu }_{q\left(x_{mij}\right)}\\ 1^T{x}_{obs}+1^T{\mu }_{q\left(x_{mbs}\right)}& {∥x_{obs}∥}^2+{∥{\mu }_{q\left(x_{mis}\right)}∥}^2+\sum _{j=1}^{n_{mis}}{\sigma }_{q\left(x_{mbs,j}\right)}^2\end{array}\right)$$

Appendix B

Table A1

Parameter estimation and CPU time cases in Simulation I using Gibbs, $n=200$.

	Method		P	BIAS		RMSE		SD		CPU Time
				${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$
with	C1	M1	0.25	0.009	0.034	0.059	0.107	0.059	0.102	7273
			0.5	0.013	0.047	0.055	0.092	0.054	0.079	8157
			0.75	0.039	0.114	0.076	0.162	0.066	0.115	8925
		M2	0.25	0.010	0.078	0.070	0.141	0.070	0.119	16,166
			0.5	0.054	0.046	0.078	0.109	0.056	0.099	16,270
			0.75	0.075	0.017	0.098	0.112	0.063	0.111	15,185
	C2	M1	0.25	0.065	0.045	0.089	0.120	0.062	0.112	8202
			0.5	0.011	0.051	0.056	0.122	0.055	0.111	8962
			0.75	0.029	0.077	0.060	0.134	0.053	0.110	9922
		M2	0.25	0.085	0.095	0.115	0.159	0.076	0.127	16,220
			0.5	0.048	0.145	0.080	0.179	0.063	0.105	15,989
			0.75	0.038	0.193	0.078	0.109	0.067	0.122	15,570
	C3	M1	0.25	0.077	0.034	0.078	0.062	0.029	0.055	8646
			0.5	0.167	0.088	0.165	0.117	0.042	0.085	9535
			0.75	0.248	0.131	0.123	0.102	0.060	0.131	8627
		M2	0.25	0.137	0.037	0.127	0.093	0.040	0.090	16,323
			0.5	0.129	0.078	0.130	0.103	0.042	0.074	16,661
			0.75	0.139	0.105	0.168	0.126	0.044	0.081	16,555
without	C1	M1	0.25	0.038	0.010	0.056	0.073	0.041	0.073	5721
			0.5	0.071	0.005	0.081	0.070	0.038	0.070	6958
			0.75	0.112	0.012	0.120	0.078	0.045	0.078	5717
		M2	0.25	0.038	0.101	0.056	0.073	0.041	0.073	12,503
			0.5	0.071	0.005	0.081	0.070	0.081	0.070	12,324
			0.75	0.112	0.012	0.120	0.073	0.045	0.078	12,901
	C2	M1	0.25	0.101	0.004	0.108	0.070	0.108	0.070	6989
			0.5	0.101	0.008	0.108	0.070	0.038	0.069	6934
			0.75	0.033	0.007	0.049	0.081	0.036	0.081	6892
		M2	0.25	0.110	0.010	0.116	0.073	0.039	0.073	12,901
			0.5	0.067	0.048	0.079	0.088	0.041	0.074	12,484
			0.75	0.083	0.046	0.092	0.100	0.042	0.089	12,560
	C3	M1	0.25	0.029	0.001	0.038	0.055	0.025	0.055	6922
			0.5	0.082	0.007	0.091	0.081	0.040	0.081	6888
			0.75	0.142	0.008	0.153	0.105	0.058	0.105	6896
		M2	0.25	0.047	0.016	0.056	0.052	0.031	0.049	12,359
			0.5	0.113	0.032	0.120	0.089	0.042	0.084	12,426
			0.75	0.175	0.044	0.187	0.118	0.065	0.110	12,462

Table A2

Parameter estimation and CPU time cases in Simulation II using Gibbs, $n=200$.

	Method		P	BIAS		RMSE		SD		CPU Time
				${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$	${\mathit{\beta }}_{\mathbf{0}}$	${\mathit{\beta }}_{\mathbf{1}}$
with	C1	M1	0.25	0.032	0.019	0.042	0.045	0.028	0.041	6159.5
			0.5	0.001	0.044	0.023	0.055	0.023	0.033	6180.5
			0.75	0.047	0.12	0.054	0.13	0.027	0.05	9567.5
		M2	0.25	0.113	0.09	0.128	0.116	0.064	0.076	9144.5
			0.5	0.037	0.087	0.04	0.1	0.03	0.051	10,985
			0.75	0.049	0.018	0.057	0.185	0.028	0.046	10,882
	C2	M1	0.25	0.025	0.027	0.026	0.034	0.01	0.021	7521
			0.5	0.007	0.02	0.012	0.024	0.009	0.013	6396.5
			0.75	0.006	0.032	0.012	0.037	0.01	0.018	7632
		M2	0.25	0.083	0.087	0.084	0.087	0.028	0.05	9013
			0.5	0.019	0.036	0.021	0.04	0.01	0.018	8655.5
			0.75	0.017	0.061	0.021	0.067	0.011	0.028	9447
	C3	M1	0.25	0.085	0.05	0.089	0.064	0.028	0.042	7583
			0.5	0.022	0.041	0.032	0.057	0.024	0.039	6798.5
			0.75	0.022	0.076	0.035	0.089	0.027	0.046	7240
		M3	0.25	0.033	0.013	0.041	0.044	0.025	0.041	11,252.5
			0.5	0.029	0.083	0.035	0.091	0.018	0.036	11,173.5
			0.75	0.047	0.011	0.053	0.117	0.023	0.044	11,873.5
without	C1	M1	0.25	0.034	0.032	0.047	0.06	0.032	0.051	4545
			0.5	0.019	0.013	0.034	0.047	0.028	0.045	4625
			0.75	0.007	0.037	0.031	0.068	0.03	0.057	5358
		M2	0.25	0.079	0.03	0.112	0.099	0.105	0.108	8190
			0.5	0.082	0.012	0.089	0.056	0.033	0.055	8599.5
			0.75	0.062	0.023	0.076	0.052	0.039	0.047	8877.5
	C2	M1	0.25	0.027	0.008	0.031	0.027	0.014	0.025	4754
			0.5	0.023	0.001	0.033	0.029	0.023	0.029	4944.5
			0.75	0.028	0.001	0.04	0.05	0.029	0.05	5175
		M2	0.25	0.04	0.003	0.044	0.031	0.018	0.032	8334.5
			0.5	0.035	0.002	0.04	0.03	0.019	0.03	8072.5
			0.75	0.036	0.004	0.049	0.043	0.031	0.043	8232
	C3	M1	0.25	0.028	0.01	0.031	0.026	0.014	0.024	4822
			0.5	0.023	0.002	0.032	0.029	0.023	0.029	4756
			0.75	0.028	0.001	0.04	0.05	0.029	0.05	4955.5
		M2	0.25	0.04	0.02	0.041	0.029	0.01	0.023	7966
			0.5	0.035	0.002	0.04	0.03	0.019	0.03	7895
			0.75	0.036	0.004	0.048	0.043	0.031	0.043	8127.5

References

1. Koenker, R.; Bassett, G. Regression quantiles. Econom. Econom. Soc.; 1987; 46, pp. 33-50.

2. Baur, D.G.; Dimpfl, T.; Jung, R.C. Stock return autocorrelations revisited: A quantile regression approach. J. Empir. Finance; 2012; 19, pp. 254-265. [DOI: https://dx.doi.org/10.1016/j.jempfin.2011.12.002]

3. Huang, L.; Zhu, W.; Saunders, C.P.; MacLeod, J.N.; Zhou, M.; Stromberg, A.J.; Bathke, A.C. A novel application of quantile regression for identification of biomarkers exemplified by equine cartilage microarray data. BMC Bioinform.; 2008; 9, 300. [DOI: https://dx.doi.org/10.1186/1471-2105-9-300] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18597687]

4. Cade, B.S.; Barry, R.N. A gentle introduction to quantile regression for ecologists. Front. Ecol. Environ.; 2003; 1, pp. 412-420. [DOI: https://dx.doi.org/10.1890/1540-9295(2003)001[0412:AGITQR]2.0.CO;2]

5. Yu, K.; Stander, J. Bayesian analysis of a Tobit quantile regression model. J. Econom.; 2007; 137, pp. 260-276. [DOI: https://dx.doi.org/10.1016/j.jeconom.2005.10.002]

6. Hideo, K.; Kobayashi, G. Gibbs sampling methods for Bayesian quantile regression. J. Stat. Comput. Simul.; 2011; 81, pp. 1565-1578.

7. Alhamzawi, R. Bayesian Analysis of Composite Quantile Regression. Stat. Biosci.; 2016; 8, pp. 358-373. [DOI: https://dx.doi.org/10.1007/s12561-016-9158-8]

8. Yuan, X.; Xiang, X.; Zhang, X. Bayesian composite quantile regression for the single-index model. PLoS ONE; 2023; 18, e0285277. [DOI: https://dx.doi.org/10.1371/journal.pone.0285277]

9. Hu, Y.; Wang, H.J.; He, X.; Guo, J. Bayesian joint-quantile regression. Comput. Stat.; 2020; 36, pp. 2033-2053. [DOI: https://dx.doi.org/10.1007/s00180-020-00998-w]

10. Li, Q.; Lin, N.; Xi, R. Bayesian regularized quantile regression. Bayesian Anal.; 2010; 5, pp. 533-556. [DOI: https://dx.doi.org/10.1214/10-BA521]

11. Alhamzawi, R.; Yu, K. Variable selection in quantile regression via Gibbs sampling. J. Appl. Stat.; 2012; 39, pp. 799-813. [DOI: https://dx.doi.org/10.1080/02664763.2011.620082]

12. Alhamzawi, R.; Mallick, H. Bayesian reciprocal LASSO quantile regression. Commun. Stat. Simul. Comput.; 2020; 51, pp. 6479-6494. [DOI: https://dx.doi.org/10.1080/03610918.2020.1804585]

13. Little, R.J.A.; Rubin, D.B. Statistical Analysis with Missing Data; 3rd ed. John Wiley & Sons: New York, NY, USA, 2019.

14. Ying, Y.; Yin, G. Bayesian Quantile Regression for Longitudinal Studies with Nonignorable Missing Data. Biometrics; 2010; 66, pp. 105-114.

15. Zhao, P.-Y.; Tang, N.-S.; Jiang, D.-P. Efficient inverse probability weighting method for quantile regression with nonignorable missing data. Statistics; 2017; 51, pp. 363-386. [DOI: https://dx.doi.org/10.1080/02331888.2016.1268615]

16. Wang, Z.; Tang, N. Bayesian Quantile Regression with Mixed Discrete and Nonignorable Missing Covariates. Bayesian Anal.; 2020; 15, pp. 579-604. [DOI: https://dx.doi.org/10.1214/19-BA1165]

17. Tang, N.; Chow, S.-M.; Ibrahim, J.G.; Zhu, H. Bayesian Sensitivity Analysis of a Nonlinear Dynamic Factor Analysis Model with Nonparametric Prior and Possible Nonignorable Missingness. Psychometrika; 2017; 82, pp. 875-903. [DOI: https://dx.doi.org/10.1007/s11336-017-9587-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29030749]

18. Tuerde, M.; Tang, N. Bayesian semiparametric approach to quantile nonlinear dynamic factor analysis models with mixed ordered and nonignorable missing data. Statistics; 2022; 56, pp. 1166-1192. [DOI: https://dx.doi.org/10.1080/02331888.2022.2121399]

19. Hoffman, M.D.; Blei, D.M.; Wang, C.; Paisley, J. Stochastic variational inference. J. Mach. Learn. Res.; 2012; 14, pp. 1303-1347.

20. Beal, M.J. Variational Algorithms for Approximate Bayesian Inference; University of London, University College London: London, UK, 2003.

21. Ganguly, A.; Jain, S.; Watchareeruetai, U. Amortized Variational Inference: Towards the Mathematical Foundation and Review. arXiv; 2022; arXiv: 2209.10888

22. Faes, C.; Ormerod, J.T.; Wand, M.P. Variational Bayesian Inference for Parametric and Nonparametric Regression with Missing Data. J. Am. Stat. Assoc.; 2011; 106, pp. 959-971. [DOI: https://dx.doi.org/10.1198/jasa.2011.tm10301]

23. Spaanberg, E. Variational Inference of Dynamic Factor Models with Arbitrary Missing Data. arXiv; 2022; arXiv: 2207.01976

24. Liu, Q.; Li, J.; Dong, M.; Liu, M.; Chai, Y. Identification of gene regulatory networks using variational bayesian inference in the presence of missing data. IEEE/ACM Trans. Comput. Biol. Bioinform.; 2022; 20, pp. 399-409. [DOI: https://dx.doi.org/10.1109/TCBB.2022.3144418] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35061589]

25. Li, Y.; Zhu, J. L1-Norm Quantile Regression. J. Comput. Graphical Stat.; 2008; 17, pp. 163-185. [DOI: https://dx.doi.org/10.1198/106186008X289155]

26. Bishop, C.M.; Nasrabadi, N.M. Pattern Recognition and Machine Learning; Springer: New York, NY, USA, 2006.

27. Wasserman, L. All of Statistics: A Concise Course in Statistical Inference; Springer: New York, NY, USA, 2004.

28. Jim, A.; Chib, S. Bayesian analysis of binary and polychotomous response data. J. Am. Stat. Assoc.; 1993; 88, pp. 669-679.

29. Polson, N.G.; Scott, J.G.; Windle, J. Bayesian Inference for Logistic Models Using Pólya–Gamma Latent Variables. J. Am. Stat. Assoc.; 2012; 108, pp. 1339-1349. [DOI: https://dx.doi.org/10.1080/01621459.2013.829001]

30. Durante, D.; Rigon, T. Conditionally Conjugate Mean-Field Variational Bayes for Logistic Models. Statist. Sci.; 2019; 34, pp. 472-485. [DOI: https://dx.doi.org/10.1214/19-STS712]