# Bayesian Adaptive Lasso for Regression Models with Nonignorable Missing Responses

Authors: Yuanying Zhao; Xingde Duan
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The main purpose of this article is to develop a Bayesian adaptive lasso procedure for analyzing linear regression models with nonignorable missing responses, in which the missingness mechanism is specified by a logistic regression model. A sampling procedure combining the Gibbs sampler and Metropolis-Hastings algorithm is employed to obtain the Bayesian estimates of the regression coefficients, shrinkage coefficients, missingness mechanism models parameters, and their standard errors. We extend the partial posterior predictivep value for goodness-of-fit statistic to investigate the plausibility of the posited model. Finally, several simulation studies and the air pollution data example are undertaken to demonstrate the newly developed methodologies.


---

## Body

## 1. Introduction

As variable selection plays an increasingly important role in modern data analysis, the least absolute shrinkage and selection operator (lasso) is appealing owing to its sparse representation. For example, see Tibshirani [1], Efron et al. [2], Zou [3], Yuan and Lin [4], Meinshausen and Buhlmann [5], etc. In the Bayesian statistical framework, Park and Casella [6] first proposed Bayesian lasso, which exploits model inference via posterior distributions. Afterwards, Guo et al. [7] studied how to select variables using Bayesian lasso in semiparametric structural equation models. Recently, Hans [8] further addressed model uncertainty and variables selection in Bayesian lasso. In particular, Bayesian adaptive lasso proposed by Leng et al. [9] permitted a treatment for variable selection with flexible penalties.On the other hand, missing data are frequently encountered in various settings, consisting of clinical trials, surveys, and longitudinal studies. Covariates and/or responses may be missing, and statistical analysis to handle the missing data is often based on the missing data mechanism [10], such as missing not at random, also referred to as nonignorable missingness. Variable selection procedures in the presence of missing data have received much attention in the recent literature. For instance, Ibrahim et al. [11] employed EM algorithm to propose a model selection criteria in missing data problems. After that, Garcia et al. [12] investigated the variable selection problem for regression models with missing covariate and/or response data. Recently, in the framework of semiparametric varying-coefficient partially linear models with missing response at random, a variable selection procedure presented by Zhao and Xue [13] can be used to simultaneously select significant variables in nonparametric and parametric model, among others. However, to the best of our knowledge, the study of missing data based on Bayesian adaptive Lasso et al. [9] has not been proposed at present. Therefore, our motivation is to construct a Bayesian method to simultaneously select variables and estimate parameters in linear regression models with nonignorable missing responses under the framework of Bayesian adaptive lasso. In our work, missing data mechanism is specified by a logistic regression model, and a MCMC procedure is proposed to obtain the Bayesian estimates of the regression coefficients, shrinkage coefficients, missingness mechanism models parameters, and their standard errors. Contributions of our work mainly include (1) tracking with the computationally intractable integral problem in posterior analysis, especially when nonignorable missing data raises new challenges for variable selection in our considered models; (2) the sampling-based approach could give rise to reliable statistical inferences and does not depend on asymptotic theory.The structure of the article is arranged as follows.We present a hierarchical Bayesian lasso model in the framework of linear regression models with nonignorable missing responses in Section2. In Section 3, a Bayesian procedure is proposed and developed by combining the Gibbs sampler [14] and the Metropolis-Hastings (M-H) algorithm [15, 16] to obtain statistical inference for simultaneous variable selection and parameter estimation. In particular, the partial posterior predictive p-value for goodness-of-fit statistic is also presented to investigate the plausibility of the posited model in Section 3. Several simulation studies and the air pollution data example are employed to demonstrate the preceding proposed Bayesian methodologies in Section 4.The article is ended up with a discussion in Section5.


## 2. Model and Notations

We take into account the linear regression model problemy=μ1n+xβ+ε, in which y=y1,y2,…,yn′ is an n×1 vector of responses, x is an n×p matrix of covariates, β is a p×1 regression parameter, and ε=ε1,ε2,…,εn′ is an n×1 vector of normal errors with mean 0n and covariance σ2In. Here, 1n (or 0n) is an n×1 vector whose each element is to set 1 (or 0), In represents the n×n identity matrix, and the superscript ′ denotes the transposition of a vector (or matrix). Without loss of generality, we just consider the case that μ is 0 and may be deleted from the above proposed model in this article since y and x can be centered. Similar to Leng et al. [9], the following adaptive lasso estimator is used to obtain statistical inference for simultaneous parameter estimation and variable selection(1)miny−xβ′y−xβ+∑j=1pλjβj.Here, flexible penalty parameters are used for different regression parameters. The following hierarchical Bayesian lasso model is considered in this article:(2)y|x,β,σ2∼Nnxβ,σ2In,β|σ2,τ12,…,τp2∼Np0p,σ2Dτ,Dτ=diagτ12,…,τp2,with the following priors on σ2 and τ2=τ12,…,τp2′:(3)σ2,τ12,…,τp2∼πσ2∏j=1pλj22exp−λj2τj22.Park and Casella [6] proposed to employ the improper prior πσ2∝1/σ2 to model the error variance, where Np.,. denotes the p-dimensional normal distribution.In this article, we focus on the case that the responsesy=y1,y2,…,yn′ are incompletely observed, which may be subject to nonignorable missingness, while covariates x=x1,x2,…,xn′ are completely observed, where xi=xi1,…,xip′. For y=y1,y2,…,yn′, we define a vector r=r1,r2,…,rn′ of missing indicators such that ri=0 if it is observed and ri=1 if yi is missing. Therefore, y can be expressed by y=yo,ym, where yo and ym denote the set of the observed data and the missing values associated with the responses, respectively. Let pri|yi,xi,φ be the conditional probability distribution of ri given yi,xi and unknown parameter φ, and then(4)pri|yi,xi,φ=pri=1|yi,xi,φri1−pri=1|yi,xi,φ1−ri.Here,pri=1|yi,xi,φ can be specified by the following logistic regression model:(5)logitpri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+⋯+φ2pxip≜φ′ωi,in which φ=φ0,φ1,φ21,…,φ2p′ and ωi=1,yi,xi1,…,xip′.Letθ=β,σ2,τ2,λ2 that contains all unknown parameter in models (1)–(3), in which τ2=τ12,…,τp2′ and λ2=λ12,…,λp2′. The main objective of this article is to develop a Bayesian procedure to analyze the above given model under the missing data indicators r and the observed data yo,x. To make Bayesian statistical inferences on the above given model, random observations are needed to be generated from the posterior distributions of θ and φ. For the conditional distribution of θ and φ given the missing data indicators r and the observed data yo,x, it follows from models (1)–(5) that pθ,φ|yo,x,r∝(6)pyo,x,r|θ,φpθ,φ∝∫∏i=1npyi|xi,θpri|yi,xi,φdympθ,φ,where pθ,φ represents the joint prior probability density for θ and φ. Obtaining a closed form of integral (6) is rather difficult due to the missing mechanism model involved.


## 3. Bayesian Analysis of the Models

For developing an effective Bayesian approach to analyze the above given model and reducing computational burden, we augment the indicatorsr and the observed data yo,x with the missing values ym to produce full data set yo,x,ym,r=y,x,r in the posterior analysis on the basis of the data augmentation idea [17]. Clearly, the posterior density pθ,φ,ym|yo,x,r is easier to handle than pθ,φ|yo,x,r.However, generating observations frompθ,φ,ym|yo,x,r is still rather difficult due to the complexity of the considered model.The Gibbs sampler [14], one of the most useful tools in statistical calculation, is employed to simulate a sequence of samples from the joint posterior distribution pθ,φ,ym|yo,x,r and then make Bayesian statistical inferences on the above given model.The key to this algorithm is that random samplesθ,φ,ym are iteratively simulated from the following three conditional distribution functions: pθ|yo,ym,x,φ,r,pφ|yo,ym,x,θ,r, and pym|yo,x,θ,φ,r. Note that pθ|yo,ym,x,φ,r=pθ|y,x and pφ|yo,ym,x,θ,r=pφ|y,x,θ,r in our considered model.


### 3.1. Gibbs Sampling and Conditional Distributions

Based on the above discussion, the Gibbs sampler [14] is implemented as follows:Step 1. Specifying the initial values for parameters θ,φ and ym, let θ0=β0,σ20,τ20,λ20, φ0=φ00,φ10,φ210,…,φ2p0′ and ym=ym0.Step 2. Based on θt=βt,σ2t,τ2t,λ2t, and yt=yo,ymt, generating sample θt+1=βt+1,σ2t+1,τ2t+1,λ2t+1 as follows:Generating sampleβt+1 from the following distribution:(7)pβ|yt,x,σ2t,τ2t∼Npμt,σ2tAt−1,whereAt=x′x+Dτ−1 and μt=At−1x′yt in which Dτ=diagτ12t,…,τp2t.Generating sampleσ2t+1 from the following distribution:(8)p1σ2|yt,x,βt+1,τ2t∼Γsh,sct,whereΓa1,a2 represents Gamma distribution with shape parameter a1 and scale parameter a2, in which sh=n+p−1/2, sct=yt−xβt+1′yt−xβt+1/2+βt+1′Dτ−1βt+1/2.Generating sampleτ2t+1 from the following distribution:(9)p1τj2|βt+1,σ2t+1,λj2t∼IGμj∗,λj∗.forj=1,2,…,p,whereIGa3,a4 represents inverse Gaussian distribution with parameters a3 and a4, in which μj∗=λjt+1σt+1/βjt+1 and λj∗=λj2t.Samplingλ2t+1:There are two methods to update parameterλ2t+1: empirical Bayesian method and hierarchical model method; see more details in Leng et al. [9].Step 3. Based on yt=yo,ymt,Generating sampleφt+1=φ0t+1,φ1t+1,φ21t+1,…,φ2pt+1′ from the following distribution:(10)pφ|yt,x,r∝pr|yt,x,φpφ∝exp∑i=1nriφ′ωi−1/2φ−φ0′V−1φ−φ0∏i=1n1+expφ′ωi,whereωi=1,yit,xi′′, pφ∼Np+2φ0,V is the prior distribution of parameter φ.Step 4. Based on θt+1, φt+1 and yt=yo,ymt, generating sample ymt+1 from the following distribution:(11)pyim|xi,θt+1,φt+1∝exp−12σ2t+1yim−xi′βt+12+φt+1Tωi−log1+expφt+1′ωi,whereωi=1,yit,xi′′.Step 5. Repeating steps 2−4.Followed by Geman and Geman [14], the sequence of observation θt,φt,ymt:t=1,2,…,T can be treated as a sample simulated from the joint posterior function pθ,φ,ym|yo,x,r. In practice, the estimated potential scale reduction (EPSR) values [18] can be employed to monitor the convergence of the algorithm.It is easy to find that conditional distributions from (7)–(9) are some standard and familiar distributions. Therefore, generating samples from these conditional distributions are direct. However, those corresponding conditional distributions given in (10) and (11) are complicated and nonstandard distributions. As a result, to directly generate samples from those conditional distributions is rather difficult. The well-known M-H algorithm [15, 16] is one of the most effective ways to generate samples from those corresponding conditional distributions with the help of proposal distributions. Similar to Zhao and Tang [19] and Xu et al. [20], the following proposal distributions Np+20,σφ2Ωφ and N10,σy2Ωy are employed for sampling φ and yim, respectively, where Ωφ−1=1/4∑i=1nωiωi′+V−1 and Ωy−1=expφ′ωi/1+expφ′ωi2φ12+1/σ2yi=0. Specifically, the M-H algorithm for sampling observation φ and yim is implemented by the following process: at the t+1 th iteration with the current value φt and yimt, new candidates φ∗ and yi∗ are generated from Np+2φt,σφ2Ωφ and N1yimt,σy2Ωy, respectively. They are accepted with probabilities min1,pφ∗|y,x,r/pφt|y,x,r and min1+pyi∗|xi,θ,φ/pyimt|xi,θ,φ, respectively. σφ2 and σy2 are chosen such that the average acceptance rate is between 0.25 and 0.45 [21].


### 3.2. Bayesian Estimates and Goodness-of-Fit Statistic

The joint Bayesian estimates of parametersθ,φ and ym are expressed by(12)θ^=1T∑t=1Tθt,φ^=1T∑t=1Tφt,y^m=1T∑t=1Tymt,where θt,φt,ymt:t=1,2,…,T are the samples that are obtained from the joint posterior distribution pθ,φ,ym|yo,x,r via the previous Gibbs sampler procedure. Following Geyer [22], these joint Bayesian estimates and sample covariance matrices are consistent with their corresponding posterior expectations and posterior covariance matrices. Therefore, we can obtain the standard errors by the diagonal elements of the sample covariance matrices.The partial posterior predictive (PPP)p-value proposed by Bayarri and Berger [23] is a goodness-of-fit statistic to investigate the plausibility of the posited model in Bayesian statistical framework; see more details in Lee and Tang [24], Tang and Zhao [25]. Similar to Lee and Tang [24] and Tang and Zhao [25], we can define the PPP p-value for our considered model as(13)pppB=∫PrDyo|θ,φ,ym,x,r≥Dyoobs|θ,φ,ym,x,rp∗θ,φ,ymdθdφdym,in which yoobs is the observed value of the responses yo and(14)p∗θ,φ,ym∝pyoobs,ym,x,r|θ,φpθ,φpDobs|θ,φ,ym,D·|· in (13) denotes a discrepancy variable, and Dobs represents the observed value of D·|· on the basis of yoobs. The following discrepancy variable can be specified in our problem:(15)Dyo|θ,φ,ym,x,r=∑i=1ndi,in which(16)di=yoi−υoi∗Voi−1yoi−υoi∗,ifri=0;0,ifri=1.where υoi∗ and Voi are conditional expectation and conditional variance of yoi given θ,φ,ym,x and r. It is easily seen that(17)p∗θ,φ,ym∝pyoobs,ym,x,r|θ,φpθ,φ∝∏i=1npyi|xi,θpri|yi,xi,θ,φpθ,φ.To reduce computational burden of the PPPp-value in (13), the Monte Carlo integral is employed to approximate pppB. Let θt,φt,ymt:t=1,2,…,T be the samples generated from the distribution p∗θ,φ,ym via the previous Gibbs sampler procedure, and let yot:t=1,…,T be the samplers simulated from pyo|ym,x,θ,φ,r. Then, the PPP p-value in (13) is approximated by(18)ppp^B=1T∑t=1TIDyot|θt,φt,ymt,x,r≥Dyoobs|θt,φt,ymt,x,r,where I· is an indicator function. According to [24, 25], if ppp^B is far away from 0.5, then we reject the posited model.


## 3.1. Gibbs Sampling and Conditional Distributions

Based on the above discussion, the Gibbs sampler [14] is implemented as follows:Step 1. Specifying the initial values for parameters θ,φ and ym, let θ0=β0,σ20,τ20,λ20, φ0=φ00,φ10,φ210,…,φ2p0′ and ym=ym0.Step 2. Based on θt=βt,σ2t,τ2t,λ2t, and yt=yo,ymt, generating sample θt+1=βt+1,σ2t+1,τ2t+1,λ2t+1 as follows:Generating sampleβt+1 from the following distribution:(7)pβ|yt,x,σ2t,τ2t∼Npμt,σ2tAt−1,whereAt=x′x+Dτ−1 and μt=At−1x′yt in which Dτ=diagτ12t,…,τp2t.Generating sampleσ2t+1 from the following distribution:(8)p1σ2|yt,x,βt+1,τ2t∼Γsh,sct,whereΓa1,a2 represents Gamma distribution with shape parameter a1 and scale parameter a2, in which sh=n+p−1/2, sct=yt−xβt+1′yt−xβt+1/2+βt+1′Dτ−1βt+1/2.Generating sampleτ2t+1 from the following distribution:(9)p1τj2|βt+1,σ2t+1,λj2t∼IGμj∗,λj∗.forj=1,2,…,p,whereIGa3,a4 represents inverse Gaussian distribution with parameters a3 and a4, in which μj∗=λjt+1σt+1/βjt+1 and λj∗=λj2t.Samplingλ2t+1:There are two methods to update parameterλ2t+1: empirical Bayesian method and hierarchical model method; see more details in Leng et al. [9].Step 3. Based on yt=yo,ymt,Generating sampleφt+1=φ0t+1,φ1t+1,φ21t+1,…,φ2pt+1′ from the following distribution:(10)pφ|yt,x,r∝pr|yt,x,φpφ∝exp∑i=1nriφ′ωi−1/2φ−φ0′V−1φ−φ0∏i=1n1+expφ′ωi,whereωi=1,yit,xi′′, pφ∼Np+2φ0,V is the prior distribution of parameter φ.Step 4. Based on θt+1, φt+1 and yt=yo,ymt, generating sample ymt+1 from the following distribution:(11)pyim|xi,θt+1,φt+1∝exp−12σ2t+1yim−xi′βt+12+φt+1Tωi−log1+expφt+1′ωi,whereωi=1,yit,xi′′.Step 5. Repeating steps 2−4.Followed by Geman and Geman [14], the sequence of observation θt,φt,ymt:t=1,2,…,T can be treated as a sample simulated from the joint posterior function pθ,φ,ym|yo,x,r. In practice, the estimated potential scale reduction (EPSR) values [18] can be employed to monitor the convergence of the algorithm.It is easy to find that conditional distributions from (7)–(9) are some standard and familiar distributions. Therefore, generating samples from these conditional distributions are direct. However, those corresponding conditional distributions given in (10) and (11) are complicated and nonstandard distributions. As a result, to directly generate samples from those conditional distributions is rather difficult. The well-known M-H algorithm [15, 16] is one of the most effective ways to generate samples from those corresponding conditional distributions with the help of proposal distributions. Similar to Zhao and Tang [19] and Xu et al. [20], the following proposal distributions Np+20,σφ2Ωφ and N10,σy2Ωy are employed for sampling φ and yim, respectively, where Ωφ−1=1/4∑i=1nωiωi′+V−1 and Ωy−1=expφ′ωi/1+expφ′ωi2φ12+1/σ2yi=0. Specifically, the M-H algorithm for sampling observation φ and yim is implemented by the following process: at the t+1 th iteration with the current value φt and yimt, new candidates φ∗ and yi∗ are generated from Np+2φt,σφ2Ωφ and N1yimt,σy2Ωy, respectively. They are accepted with probabilities min1,pφ∗|y,x,r/pφt|y,x,r and min1+pyi∗|xi,θ,φ/pyimt|xi,θ,φ, respectively. σφ2 and σy2 are chosen such that the average acceptance rate is between 0.25 and 0.45 [21].


## 3.2. Bayesian Estimates and Goodness-of-Fit Statistic

The joint Bayesian estimates of parametersθ,φ and ym are expressed by(12)θ^=1T∑t=1Tθt,φ^=1T∑t=1Tφt,y^m=1T∑t=1Tymt,where θt,φt,ymt:t=1,2,…,T are the samples that are obtained from the joint posterior distribution pθ,φ,ym|yo,x,r via the previous Gibbs sampler procedure. Following Geyer [22], these joint Bayesian estimates and sample covariance matrices are consistent with their corresponding posterior expectations and posterior covariance matrices. Therefore, we can obtain the standard errors by the diagonal elements of the sample covariance matrices.The partial posterior predictive (PPP)p-value proposed by Bayarri and Berger [23] is a goodness-of-fit statistic to investigate the plausibility of the posited model in Bayesian statistical framework; see more details in Lee and Tang [24], Tang and Zhao [25]. Similar to Lee and Tang [24] and Tang and Zhao [25], we can define the PPP p-value for our considered model as(13)pppB=∫PrDyo|θ,φ,ym,x,r≥Dyoobs|θ,φ,ym,x,rp∗θ,φ,ymdθdφdym,in which yoobs is the observed value of the responses yo and(14)p∗θ,φ,ym∝pyoobs,ym,x,r|θ,φpθ,φpDobs|θ,φ,ym,D·|· in (13) denotes a discrepancy variable, and Dobs represents the observed value of D·|· on the basis of yoobs. The following discrepancy variable can be specified in our problem:(15)Dyo|θ,φ,ym,x,r=∑i=1ndi,in which(16)di=yoi−υoi∗Voi−1yoi−υoi∗,ifri=0;0,ifri=1.where υoi∗ and Voi are conditional expectation and conditional variance of yoi given θ,φ,ym,x and r. It is easily seen that(17)p∗θ,φ,ym∝pyoobs,ym,x,r|θ,φpθ,φ∝∏i=1npyi|xi,θpri|yi,xi,θ,φpθ,φ.To reduce computational burden of the PPPp-value in (13), the Monte Carlo integral is employed to approximate pppB. Let θt,φt,ymt:t=1,2,…,T be the samples generated from the distribution p∗θ,φ,ym via the previous Gibbs sampler procedure, and let yot:t=1,…,T be the samplers simulated from pyo|ym,x,θ,φ,r. Then, the PPP p-value in (13) is approximated by(18)ppp^B=1T∑t=1TIDyot|θt,φt,ymt,x,r≥Dyoobs|θt,φt,ymt,x,r,where I· is an indicator function. According to [24, 25], if ppp^B is far away from 0.5, then we reject the posited model.


## 4. Numerical Examples

Several simulation studies and the air pollution data example are employed to demonstrate our proposed methodologies in this section.


### 4.1. Example 1

The first simulation study is employed to illustrate the performance of the Bayesian estimates with respective to the sparse, dense, and very sparse cases, respectively. Specifically, we generate observationsyi from the following model:(19)yi=xi′β+εi=xi1β1+⋯+xipβp+εi.Fori=1,2,…,n, here, the true values for parameter β are given as follows:Case i (moderate case):β=3,1.5,0,0,2,0,0,0T;Case ii (dense case):β=1,1,1,1,1,1,1,1T;Case iii (sparse case):β=5,0,0,0,0,0,0,0T.Fori=1,2,…,n, j,k=1,2,…,8, it is assumed that xij follows standard normal distribution N0,1 marginally and the correlation between xij and xik is 0.5j−k, and εi is independent and identically distributed as N10,σ2, and the true value of σ2 is set to 1 in all simulation studies. Missing data are generated from the logistic regression model(20)M0:logitpri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+⋯+φ28xi8,with true values φ0=−1.0, φ1=φ21=⋯=φ28=0.3. The missing data are generated as follows: (a) simulate the full data set yi:i=1,2,…,n from model (19); (b) identify whether the observation yi is missing or not via the logistic model (20) with the true values for φ0,φ1,φ21,…,φ28. More specifically, a random number κ is obtained from the uniform distribution U0,1; the observation yi is missing if κ≤prri=1|yi,xi,φ=1/1+exp−φ0−φ1yi−φ21xi1−⋯−φ28xi8.For each of the above generated data sets under various model assumptions and parameter settings, our proposed methodologies combining the Gibbs sampler together with the M-H algorithm are employed to calculate the Bayesian estimates of unknown parameters on the basis of 100 replications forn=200, and the average proportion of missing data simulated are roughly 0.375. In each replication, we collected 10 000 observations after 10 000 iterations in producing Bayesian results. In M-H algorithm, we take σy2=4.0 and σφ2=1.2 obtaining average acceptance rates 29.6% and 27.5%, respectively. The corresponding results of Bayesian estimates are reported in Tables 1–3, in which ‘True’ represents the true value of parameter; “Mean,” “SD,” and “Median” denote the average, standard deviation, and median values of the Bayesian estimates on the basis of 100 data sets, respectively; “5%” and ‘95%’ are the 5% and 95% quantiles of the Bayesian estimates on the basis of 100 replications; ‘RMSE’ is the root mean square error between the Bayesian estimates on the basis of 100 data sets and its true value. Besides, to summarize the variable selection results, in Table 4, we calculate the proportion that parameter β was identified to be zero in 100 replications in terms of the criterion that a parameter is identified to be 0 if its 95% confidence interval contains zero. Inspections of Tables 1–4 indicate that (1) the Bayesian estimates of unknown parameters for β, σ and φ are quite close to their true values under moderate case, dense case, and sparse case, which illustrates that our proposed Bayesian approach is quite effective; (2) the differences between ‘SD’ and ‘RMSE’ are very small, which shows that the Bayesian estimates are rather reliable and robust to our considered cases. (3) Our proposed methods could identify the correct models in most cases, because the proportion values in Table 4 corresponding to the important covariates are all zero, and the proportion values corresponding to unimportant covariates are more than 90%. (4) It is interesting to observe that the hierarchical model approach tends to obtain better results than empirical Bayesian approach in terms of variable selection, but the differences are minor, especially in the very dense case (e.g., case ii). To sum up, the above findings indicate that our preceding proposed Bayesian procedures could well recover the true information in models.Table 1
Bayesian estimates for casei in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9972.9912.7873.1470.1040.104β21.51.4761.4731.2821.6440.1130.112β30.00.0100.001−0.1120.1300.0670.067β40.00.001−0.000−0.0900.0970.0620.062β52.01.9921.9811.8362.1800.1110.111β60.00.0070.000−0.1110.0930.0690.069β70.00.001−0.000−0.1060.0100.0630.062β80.0−0.0070.000−0.1340.0450.0500.050σ1.00.9990.9980.8871.0890.0640.063φ0−1.0−1.091−1.102−1.348−0.8260.1690.191φ10.30.3680.3530.1130.6310.1610.174φ210.30.2650.244−0.5200.9580.4560.456φ220.30.3370.327−0.3290.9220.3610.361φ230.30.3620.377−0.1220.7680.2740.280φ240.30.3290.294−0.1930.8110.3020.302φ250.30.3170.269−0.2560.8330.3380.338φ260.30.3620.373−0.2520.9010.3050.310φ270.30.3400.316−0.1670.8030.2940.295φ280.30.3450.342−0.0980.7490.2720.274(ii) Hierarchical model scenarioβ13.03.0052.9872.8293.1490.1110.110β21.51.5001.5111.3191.6460.1130.113β30.00.0070.001−0.0590.0660.0470.048β40.0−0.0040.001−0.1190.0670.0690.069β52.02.0031.9961.8342.1570.0100.099β60.0−0.005−0.002−0.0750.0580.0410.041β70.0−0.003−0.004−0.0770.0610.0400.039β80.0−0.003−0.003−0.0490.0390.0320.032σ1.01.0191.0200.9081.1080.0630.066φ0−1.0−1.091−1.102−1.358−0.8530.1610.184φ10.30.3570.3570.0530.6100.1590.168φ210.30.2880.257−0.6300.9460.4580.456φ220.30.2810.266−0.2240.9170.3470.346φ230.30.3840.360−0.1840.8510.3060.316φ240.30.3600.337−0.2120.7800.2960.301φ250.30.2990.348−0.3630.7770.3320.330φ260.30.3280.308−0.2050.7820.2860.286φ270.30.3840.417−0.1260.7990.2840.295φ280.30.2990.311−0.0870.6440.2240.223Table 2
Bayesian estimates for case ii in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ11.01.0141.0230.8151.1630.0990.100β21.00.9940.9970.7831.1850.1290.128β31.00.9830.9740.7881.1570.1090.110β41.01.0151.0200.8191.1800.1150.116β51.01.0041.0130.8361.1560.1050.105β61.00.9920.9860.7881.1640.1140.113β71.00.9910.9790.8111.1750.1160.115β81.00.9890.9830.8291.1250.0960.096σ1.01.0101.0060.9151.1310.0650.065φ0−1.0−1.092−1.083−1.395−0.8410.1720.195φ10.30.3510.345−0.0450.6920.2270.231φ210.30.3170.304−0.1630.7650.3170.316φ220.30.3310.278−0.1850.9190.3340.334φ230.30.3710.365−0.1730.8920.3380.344φ240.30.3100.286−0.3750.9960.3880.387φ250.30.3320.302−0.2200.8310.3380.338φ260.30.3160.309−0.2500.8690.3620.361φ270.30.3630.329−0.3100.9040.3940.397φ280.30.3180.286−0.2090.7470.3270.326(ii) Hierarchical model scenarioβ11.00.9991.0000.8161.1630.1120.112β21.00.9790.9630.7981.1920.1240.125β31.01.0170.9990.8421.2240.1250.125β41.01.0010.9820.7861.2240.1290.128β51.00.9871.0040.7671.1430.1220.122β61.00.9900.9730.7611.2370.1390.139β71.00.9930.9910.8141.1650.1140.114β81.00.9961.0020.8071.1650.1110.110σ1.01.0191.0210.9151.1290.0670.069φ0−1.0−1.087−1.094−1.356−0.8430.1670.188φ10.30.3750.3430.0510.7120.2090.221φ210.30.3610.389−0.1680.7390.2950.300φ220.30.3250.340−04820.8490.3870.386φ230.30.3250.360−0.2870.7880.3320.332φ240.30.3230.334−0.3040.9030.3860.385φ250.30.2920.279−0.3420.8470.3600.359φ260.30.2950.282−0.2770.9030.3490.347φ270.30.3400.356−0.2560.8540.3520.353φ280.30.3110.302−0.2400.7800.3060.304Table 3
Bayesian estimates for case iii in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ15.05.0034.9924.8265.1850.1110.111β20.0−0.011−0.000−0.1530.0460.0680.068β30.00.003−0.001−0.0600.0530.0540.054β40.00.012−0.001−0.1930.0850.0760.076β50.00.0080.001−0.0780.0910.0650.066β60.00.0070.001−0.0830.0820.1110.111β70.00.0040.001−0.1070.0800.0530.052β80.00.0000.000−0.1190.0780.0640.063σ1.00.9910.9840.8851.0800.0580.058φ0−1.0−1.065−1.071−1.310−0.7980.1630.175φ10.30.3590.3610.0980.5290.1230.136φ210.30.2730.243−0.4420.9580.4430.442φ220.30.3230.288−0.0730.7460.2850.285φ230.30.3670.3240.0060.7450.2440.252φ240.30.3180.307−0.1710.7100.2620.262φ250.30.3460.296−0.1090.8220.2960.298φ260.30.3220.339−0.1870.7830.2970.296φ270.30.3260.300−0.0230.6680.2280.229φ280.30.3420.341−0.0380.7220.2340.237(ii) Hierarchical model scenarioβ15.05.0095.0064.8755.1630.0910.091β20.0−0.0000.001−0.0580.0370.0380.038β30.0−0.001−0.002−0.0460.0350.0310.031β40.00.004−0.000−0.0420.0480.0290.029β50.00.0020.002−0.0500.0630.0360.036β60.0−0.0010.000−0.0580.0370.0380.038β70.00.0060.005−0.0430.0550.0300.030β80.00.0010.003−0.0440.0430.0290.029σ1.01.0131.0090.9271.1060.0600.061φ0−1.0−1.089−1.085−1.348−0.8610.1500.174φ10.30.3520.3480.1490.5150.1180.129φ210.30.2890.323−0.4831.0310.4690.467φ220.30.3480.338−0.2120.7890.2810.283φ230.30.3690.362−0.0620.7350.2460.254φ240.30.2830.256−0.1760.7310.2840.283φ250.30.3530.343−0.0850.7690.2600.264φ260.30.3580.326−0.1290.7700.2620.267φ270.30.3640.352−0.1470.7590.2800.286φ280.30.2920.270−0.1170.7390.2610.260Table 4
The variable selection proportions results forβ in the first simulation study.
Caseβ1β2β3β4β5β6β7β8(i) Empirical Bayesian scenarioi0.000.000.960.970.000.950.940.98ii0.000.000.000.000.000.000.000.00iii0.000.960.980.930.940.960.970.94(ii) Hierarchical model scenarioi0.000.000.990.980.001.001.001.00ii0.000.000.000.000.000.000.000.00iii0.001.000.991.001.001.001.001.00


### 4.2. Example 2

In order to investigate the sensitivity to proportion of missing data, we conduct the second simulation study in this part. Simulated data sets are generated by using the same setup as specified in the first simulation study except that true value ofφ is set as φ0=−2.5, φ1=φ21=⋯=φ28=0.3. We fit these simulated data sets under cases i–iii via the preceding developed MCMC algorithm and only present results in Table 5 on the basis of case i under 100 replications to save space, where the average proportion of missing data generated in the second simulation study is about 0.208. Besides, we further calculate sums of RMSE under case i based on these two simulation studies, and the results are 3.365 and 3.237 in this simulation study, compared with 3.681 and 3.522 in the first simulation study. Moreover, the median of mean absolute deviations (MMAD) is calculated for 100 replications, i.e., median 1/n∑i=1nxi′β^−xi′βtrue, in which β^ and βtrue represent the Bayesian estimate and true value of parameter β.Table 5
Bayesian estimates for casei in the second simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9882.9972.8463.0960.0830.083β21.51.4821.4731.3281.6380.0890.091β30.0−0.005−0.001−0.0980.1180.0510.051β40.00.0060.000−0.1080.0950.0630.063β52.01.9921.9921.8582.1240.0920.092β60.0−0.008−0.001−0.1580.0560.0640.064β70.0−0.001−0.000−0.0950.0790.0560.055β80.00.0020.000−0.0980.1180.0510.051σ1.00.9890.9860.9141.0790.0580.058φ0−2.5−2.564−2.567−2.689−2.4520.0720.096φ10.30.3390.3450.1040.5250.1350.140φ210.30.2820.261−0.3950.8500.4130.411φ220.30.2740.258−0.2610.8010.3550.354φ230.30.3000.287−0.2050.7100.2840.282φ240.30.3310.329−0.0680.7310.2660.267φ250.30.2830.279−0.2690.7620.3270.326φ260.30.2930.294−0.2550.6980.3040.302φ270.30.2910.311−0.1960.7390.2820.281φ280.30.3870.337−0.0150.8530.2860.298(ii) Hierarchical model scenarioβ13.02.9842.9832.8223.1290.1010.102β21.51.5061.5101.3451.6540.1010.101β30.00.0080.007−0.0670.0890.0470.047β40.0−0.009−0.006−0.0820.0410.0430.043β52.01.9911.9761.8482.1390.1010.101β60.0−0.0010.003−0.0800.0480.0420.042β70.00.001−0.000−0.0750.0770.0460.046β80.00.0070.004−0.0540.0640.0410.042σ1.00.9980.9950.9211.0860.0550.054φ0−2.5−2.573−2.576−2.681−2.4700.0700.101φ10.30.3270.3150.1330.5010.1190.121φ210.30.2870.276−0.3400.7560.3340.333φ220.30.2680.282−0.2600.7100.3200.320φ230.30.3420.314−0.1260.8340.2930.294φ240.30.3070.321−0.2070.7970.3080.307φ250.30.2500.221−0.2510.7350.2810.284φ260.30.3400.354−0.1250.7120.2830.285φ270.30.3230.286−0.1870.8470.3360.335φ280.30.2970.311−0.1850.7070.2810.279The related results of MMAD for the first (second) simulation are 0.143 7 (0.129 3) and 0.134 8 (0.122 1) based on Empirical Bayesian method and Hierarchical model method, respectively. Meanwhile, the box plots for mean absolute deviations (MAD) are shown in Figure1, where “1” and “2” denote case i in the first and second simulation studies, respectively. All these results indicate that cases with the smaller proportion of missing data in this simulation study generally perform better than cases with the larger proportion of missing data in the first simulation study.Figure 1
Box plots summarizing the MAD for cases under (a) empirical Bayesian scenario (left panel) and (b) hierarchical model scenario (right panel) in the first and second simulation studies.
(a)(b)


### 4.3. Example 3

To address the influence of the missingness mechanism model to the Bayesian adaptive lasso, we conduct the third simulation in this example. The full data set is created by employing the same setting specified under case i in the first simulation study. But, we generate the missing data via the following 3 types of missingness mechanisms.Type i: logistic regression model that is different fromM0 specified in (20):(21)logitprri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+φ22xi2+φ23xi3,with φ0=−1.0 and φ1=φ21=φ22=φ23=0.3.Type ii: missing at random (MAR) missingness mechanism.Type iii: nonignorable missingness mechanism model specified inM0.For missing data simulated via type i and type ii, we obtain the Bayesian estimates under incorrect logistic regression modelM0. For missing data simulated via type iii, we calculate the Bayesian estimates under incorrect MAR missingness mechanism. Summaries of these results are reported in Tables 6–9. Inspecting these results, the following observations are indicated: (1) results produced using the model M0 are quite robust and accurate even if the true missing data are MAR, or the true missingness mechanisms are much simpler models; (2) performances obtained on the basis of the incorrect MAR assumption are poor. These observations agreed with those obtained by Lee and Tang [24].Table 6
Bayesian estimates for typei in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.03.0052.9842.8423.1920.1110.111β21.51.5021.5091.3051.6390.1110.111β30.00.0040.000−0.1000.1200.0770.077β40.0−0.004−0.001−0.1070.0720.0610.061β52.01.9781.9881.8152.1290.1010.103β60.00.0140.001−0.0390.0910.0540.056β70.00.0020.000−0.1000.1120.0570.057β80.00.0050.001−0.0620.0970.0580.058σ1.00.9940.9890.8871.0810.0590.059φ0−1.0−1.087−1.089−1.383−0.8190.1740.193φ10.30.3660.3840.0870.5930.1590.171φ210.30.3360.324−0.4100.9550.4410.440φ220.30.2580.292−0.4480.7590.3720.371φ230.30.3570.3290.0380.6670.2110.217φ240.00.027−0.007−0.3800.4370.2600.260φ250.0−0.066−0.098−0.6770.4780.3370.342φ260.00.013−0.015−0.4560.3110.2490.248φ270.0−0.018−0.056−0.4870.4100.2820.281φ280.0−0.015−0.004−0.3930.3480.2470.246(ii) Hierarchical model scenarioβ13.03.0103.0012.8133.1850.1200.120β21.51.4681.4571.2821.6330.1100.114β30.00.0280.018−0.1480.1650.1030.106β40.0−0.0040.002−0.1910.1500.1090.108β52.01.9911.9691.8122.2180.1260.126β60.0−0.003−0.001−0.1780.1420.1020.102β70.0−0.007−0.008−0.1610.1400.1040.103β80.00.002−0.007−0.1570.1650.0940.093σ1.01.0191.0200.9081.1080.0630.066φ0−1.0−1.052−1.055−1.290−0.8490.1430.151φ10.30.3190.3000.0310.5810.1750.175φ210.30.3790.351−0.4941.0810.5030.507φ220.30.3710.364−0.2290.9070.3380.344φ230.30.3430.314−0.0820.7090.2370.240φ240.00.0210.012−0.3810.4260.2550.255φ250.00.0200.006−0.6010.5660.3720.371φ260.00.0320.038−0.4970.3990.2720.271φ270.0−0.054−0.071−0.5120.3510.2700.275φ280.00.0740.065−0.2850.4130.2170.229Table 7
Bayesian estimates for type ii in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9702.9762.7453.1180.1100.113β21.51.4881.4771.3241.6390.1090.110β30.0−0.003−0.001−0.1290.0840.0670.066β40.00.008−0.000−0.0390.0860.0560.057β52.01.9681.9741.7932.1000.0940.099β60.00.0040.000−0.0830.0740.1100.113β70.00.003−0.000−0.0810.0870.0500.050β80.00.0050.000−0.0840.0950.0490.049σ1.01.0010.9960.8941.1060.0640.064φ0—−0.622−0.613−0.888−0.3780.160—φ1—−0.052−0.057−0.3170.1680.151—φ21—0.1990.194−0.6020.8940.476—φ22—0.1090.103−0.4120.5980.314—φ23—−0.013−0.022−0.3830.3110.223—φ24—−0.012−0.001−0.3900.3040.217—φ25—0.1070.098−0.4930.6460.335—φ26—0.0150.026−0.3330.2650.120—φ27—0.006−0.010−0.3810.3370.215—φ28—0.0240.013−0.2860.3640.204—(ii) Hierarchical model scenarioβ13.02.9742.9692.8163.1070.0950.098β21.51.4791.4861.2611.6560.1230.124β30.00.0040.000−0.0660.0700.0470.047β40.00.0070.001−0.0660.0990.0560.056β52.01.9731.9681.8062.1180.1010.104β60.00.0080.003−0.0680.0980.0540.055β70.00.0090.004−0.0510.0470.0490.050β80.0−0.003−0.003−0.0690.0400.0490.049σ1.01.0141.0110.9031.1150.0660.068φ0—−0.606−0.603−0.842−0.3740.152—φ1—−0.054−0.078−0.3610.1680.152—φ21—0.1680.135−0.5520.8630.459—φ22—0.0970.075−0.4290.6010.306—φ23—0.0120.037−0.3420.2960.207—φ24—−0.006−0.029−0.3160.2320.176—φ25—0.1350.157−0.5220.6970.344—φ26—−0.037−0.049−0.3750.2990.200—φ27—0.0500.050−0.3480.3930.234—φ28—−0.0330.029−0.3970.2670.202—Table 8
Bayesian estimates for type iii in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9002.9820.0013.6060.7370.740β21.51.3201.4370.0001.7990.6440.666β30.00.1030.003−0.0850.6170.2870.304β40.00.0550.000−0.0820.1810.2400.245β52.01.7451.9640.0012.1280.6400.686β60.00.0610.000−0.0670.4640.2100.218β70.00.007−0.000−0.0130.0190.0420.043β80.0−0.003−0.000−0.0520.0210.0350.035σ1.01.2701.0330.8892.3050.5320.594(ii) Hierarchical model scenarioβ13.03.1953.1852.9743.3950.1400.240β21.51.6291.6061.4141.8810.1430.192β30.00.0500.034−0.1270.3000.1310.139β40.00.0680.048−0.1640.2990.1260.142β52.02.1392.1351.9002.3540.1370.195β60.00.0530.040−0.1570.2700.1230.133β70.00.0520.048−0.1470.2450.1280.137β80.00.0630.058−0.0970.2170.0990.117σ1.01.2501.2421.1451.3640.0750.261


### 4.4. A Real Example

In order to further demonstrate our proposed Bayesian methodologies, we consider a real example from the air pollution data in this subsection, which is available at the websitehttp://lib.stat.cmu.edu/datasets/NO2.data. The air pollution data includes n=500 randomly chosen responses and p=7 explanatory variables. For i=1,2,…,n, the response variable yi is hourly values of the logarithm of the concentration of NO2, measured at Alnabru in Oslo, Norway, between October 2001 and August 2003; the explanatory variable xi1 denotes the logarithm of the number of cars per hour, xi2 represents temperature two meters above ground, xi3 is wind speed, xi4 denotes the temperature difference between 25 and 2 m above ground, xi5 represents wind direction, and xi6 and xi7 are hours of day and day number from October 1st,2001, respectively. Each explanatory variable is standardized to have mean 0 and variance 1 during exploratory data analysis. The air pollution data has been analyzed and fitted by Zhu and Zhu [26], and they observed that all explanatory variables except xi5 and xi7 make a significant contribution to yi. Model equation (19) and the following model(22)logitprri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+⋯+φ27xi7,where φ0=−1.0, φ1=φ21=⋯=φ27=0.5 are used to reanalyze the air pollution data in the article, where the logistic regression model (22) may represent that we create artificial missing data by deliberated deleting some of responses according to the given missing mechanism, and the missing data’s proportion generated in this way is 0.326. Then, we employ the preceding propose Bayesian procedure to make statistical inferences on the assumption of the above model settings. Here, 10 000 samples after 10 000 burn-in iterations are employed to calculate the Bayesian estimates of unknown parameters. In M-H algorithm, we take σy2=2.2 and σφ2=1.0 giving average acceptance rates 29.0% and 26.8%, respectively. Also, we conduct a real data analysis under the assumption for the fully observed data in the framework of Bayesian adaptive Lasso.Table9 reports the corresponding results for regression coefficients by two different methods, in which ‘EB’ and ‘HM’ denote the Bayesian estimates on the basis of Empirical Bayesian method and Hierarchical model method, respectively. From Table 9, we clearly see that (1) the estimate results for the regression coefficients in all two methods are quite close; (2) almost all standard errors obtained under the fully observed data are smaller than those obtained under the missing data, which indicates that the procedures under the fully observed data outperform the procedures under the missing data, but the differences for performances are minor; (3) all these procedures could identify that non-zero significant regression coefficients are β1,β2,β3,β4, and β6, which is consistent with the results given in Zhu and Zhu [26]. In the end, PPP p-value for goodness-of-fit with respect to model (19) with the nonignorable missing mechanism defined in model (20) is 0.446 and 0.454 based on Empirical Bayesian method and Hierarchical model method, respectively, which illustrates that our considered models as well as the proposed goodness-of-fit statistic are feasible.Table 9
Bayesian estimates and standard errors (in the parentheses) of parameterβ in the air pollution data analysis.
Missing dataFully observed dataEBHMEBHMβ10.547 (0.038)0.552 (0.039)0.541 (0.031)0.540 (0.031)β2−0.156 (0.035)−0.163 (0.036)−0.146 (0.028)−0.145 (0.029)β3−0.273 (0.035)−0.268 (0.036)−0.225 (0.025)−0.225 (0.025)β40.179 (0.038)0.200 (0.041)0.172 (0.028)0.171 (0.028)β50.034 (0.028)0.046 (0.035)0.034 (0.025)0.033 (0.028)β6−0.121 (0.038)−0.121 (0.040)−0.123 (0.030)−0.122 (0.031)β70.000 (0.002)0.053 (0.035)0.039 (0.025)0.037 (0.027)


## 4.1. Example 1

The first simulation study is employed to illustrate the performance of the Bayesian estimates with respective to the sparse, dense, and very sparse cases, respectively. Specifically, we generate observationsyi from the following model:(19)yi=xi′β+εi=xi1β1+⋯+xipβp+εi.Fori=1,2,…,n, here, the true values for parameter β are given as follows:Case i (moderate case):β=3,1.5,0,0,2,0,0,0T;Case ii (dense case):β=1,1,1,1,1,1,1,1T;Case iii (sparse case):β=5,0,0,0,0,0,0,0T.Fori=1,2,…,n, j,k=1,2,…,8, it is assumed that xij follows standard normal distribution N0,1 marginally and the correlation between xij and xik is 0.5j−k, and εi is independent and identically distributed as N10,σ2, and the true value of σ2 is set to 1 in all simulation studies. Missing data are generated from the logistic regression model(20)M0:logitpri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+⋯+φ28xi8,with true values φ0=−1.0, φ1=φ21=⋯=φ28=0.3. The missing data are generated as follows: (a) simulate the full data set yi:i=1,2,…,n from model (19); (b) identify whether the observation yi is missing or not via the logistic model (20) with the true values for φ0,φ1,φ21,…,φ28. More specifically, a random number κ is obtained from the uniform distribution U0,1; the observation yi is missing if κ≤prri=1|yi,xi,φ=1/1+exp−φ0−φ1yi−φ21xi1−⋯−φ28xi8.For each of the above generated data sets under various model assumptions and parameter settings, our proposed methodologies combining the Gibbs sampler together with the M-H algorithm are employed to calculate the Bayesian estimates of unknown parameters on the basis of 100 replications forn=200, and the average proportion of missing data simulated are roughly 0.375. In each replication, we collected 10 000 observations after 10 000 iterations in producing Bayesian results. In M-H algorithm, we take σy2=4.0 and σφ2=1.2 obtaining average acceptance rates 29.6% and 27.5%, respectively. The corresponding results of Bayesian estimates are reported in Tables 1–3, in which ‘True’ represents the true value of parameter; “Mean,” “SD,” and “Median” denote the average, standard deviation, and median values of the Bayesian estimates on the basis of 100 data sets, respectively; “5%” and ‘95%’ are the 5% and 95% quantiles of the Bayesian estimates on the basis of 100 replications; ‘RMSE’ is the root mean square error between the Bayesian estimates on the basis of 100 data sets and its true value. Besides, to summarize the variable selection results, in Table 4, we calculate the proportion that parameter β was identified to be zero in 100 replications in terms of the criterion that a parameter is identified to be 0 if its 95% confidence interval contains zero. Inspections of Tables 1–4 indicate that (1) the Bayesian estimates of unknown parameters for β, σ and φ are quite close to their true values under moderate case, dense case, and sparse case, which illustrates that our proposed Bayesian approach is quite effective; (2) the differences between ‘SD’ and ‘RMSE’ are very small, which shows that the Bayesian estimates are rather reliable and robust to our considered cases. (3) Our proposed methods could identify the correct models in most cases, because the proportion values in Table 4 corresponding to the important covariates are all zero, and the proportion values corresponding to unimportant covariates are more than 90%. (4) It is interesting to observe that the hierarchical model approach tends to obtain better results than empirical Bayesian approach in terms of variable selection, but the differences are minor, especially in the very dense case (e.g., case ii). To sum up, the above findings indicate that our preceding proposed Bayesian procedures could well recover the true information in models.Table 1
Bayesian estimates for casei in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9972.9912.7873.1470.1040.104β21.51.4761.4731.2821.6440.1130.112β30.00.0100.001−0.1120.1300.0670.067β40.00.001−0.000−0.0900.0970.0620.062β52.01.9921.9811.8362.1800.1110.111β60.00.0070.000−0.1110.0930.0690.069β70.00.001−0.000−0.1060.0100.0630.062β80.0−0.0070.000−0.1340.0450.0500.050σ1.00.9990.9980.8871.0890.0640.063φ0−1.0−1.091−1.102−1.348−0.8260.1690.191φ10.30.3680.3530.1130.6310.1610.174φ210.30.2650.244−0.5200.9580.4560.456φ220.30.3370.327−0.3290.9220.3610.361φ230.30.3620.377−0.1220.7680.2740.280φ240.30.3290.294−0.1930.8110.3020.302φ250.30.3170.269−0.2560.8330.3380.338φ260.30.3620.373−0.2520.9010.3050.310φ270.30.3400.316−0.1670.8030.2940.295φ280.30.3450.342−0.0980.7490.2720.274(ii) Hierarchical model scenarioβ13.03.0052.9872.8293.1490.1110.110β21.51.5001.5111.3191.6460.1130.113β30.00.0070.001−0.0590.0660.0470.048β40.0−0.0040.001−0.1190.0670.0690.069β52.02.0031.9961.8342.1570.0100.099β60.0−0.005−0.002−0.0750.0580.0410.041β70.0−0.003−0.004−0.0770.0610.0400.039β80.0−0.003−0.003−0.0490.0390.0320.032σ1.01.0191.0200.9081.1080.0630.066φ0−1.0−1.091−1.102−1.358−0.8530.1610.184φ10.30.3570.3570.0530.6100.1590.168φ210.30.2880.257−0.6300.9460.4580.456φ220.30.2810.266−0.2240.9170.3470.346φ230.30.3840.360−0.1840.8510.3060.316φ240.30.3600.337−0.2120.7800.2960.301φ250.30.2990.348−0.3630.7770.3320.330φ260.30.3280.308−0.2050.7820.2860.286φ270.30.3840.417−0.1260.7990.2840.295φ280.30.2990.311−0.0870.6440.2240.223Table 2
Bayesian estimates for case ii in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ11.01.0141.0230.8151.1630.0990.100β21.00.9940.9970.7831.1850.1290.128β31.00.9830.9740.7881.1570.1090.110β41.01.0151.0200.8191.1800.1150.116β51.01.0041.0130.8361.1560.1050.105β61.00.9920.9860.7881.1640.1140.113β71.00.9910.9790.8111.1750.1160.115β81.00.9890.9830.8291.1250.0960.096σ1.01.0101.0060.9151.1310.0650.065φ0−1.0−1.092−1.083−1.395−0.8410.1720.195φ10.30.3510.345−0.0450.6920.2270.231φ210.30.3170.304−0.1630.7650.3170.316φ220.30.3310.278−0.1850.9190.3340.334φ230.30.3710.365−0.1730.8920.3380.344φ240.30.3100.286−0.3750.9960.3880.387φ250.30.3320.302−0.2200.8310.3380.338φ260.30.3160.309−0.2500.8690.3620.361φ270.30.3630.329−0.3100.9040.3940.397φ280.30.3180.286−0.2090.7470.3270.326(ii) Hierarchical model scenarioβ11.00.9991.0000.8161.1630.1120.112β21.00.9790.9630.7981.1920.1240.125β31.01.0170.9990.8421.2240.1250.125β41.01.0010.9820.7861.2240.1290.128β51.00.9871.0040.7671.1430.1220.122β61.00.9900.9730.7611.2370.1390.139β71.00.9930.9910.8141.1650.1140.114β81.00.9961.0020.8071.1650.1110.110σ1.01.0191.0210.9151.1290.0670.069φ0−1.0−1.087−1.094−1.356−0.8430.1670.188φ10.30.3750.3430.0510.7120.2090.221φ210.30.3610.389−0.1680.7390.2950.300φ220.30.3250.340−04820.8490.3870.386φ230.30.3250.360−0.2870.7880.3320.332φ240.30.3230.334−0.3040.9030.3860.385φ250.30.2920.279−0.3420.8470.3600.359φ260.30.2950.282−0.2770.9030.3490.347φ270.30.3400.356−0.2560.8540.3520.353φ280.30.3110.302−0.2400.7800.3060.304Table 3
Bayesian estimates for case iii in the first simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ15.05.0034.9924.8265.1850.1110.111β20.0−0.011−0.000−0.1530.0460.0680.068β30.00.003−0.001−0.0600.0530.0540.054β40.00.012−0.001−0.1930.0850.0760.076β50.00.0080.001−0.0780.0910.0650.066β60.00.0070.001−0.0830.0820.1110.111β70.00.0040.001−0.1070.0800.0530.052β80.00.0000.000−0.1190.0780.0640.063σ1.00.9910.9840.8851.0800.0580.058φ0−1.0−1.065−1.071−1.310−0.7980.1630.175φ10.30.3590.3610.0980.5290.1230.136φ210.30.2730.243−0.4420.9580.4430.442φ220.30.3230.288−0.0730.7460.2850.285φ230.30.3670.3240.0060.7450.2440.252φ240.30.3180.307−0.1710.7100.2620.262φ250.30.3460.296−0.1090.8220.2960.298φ260.30.3220.339−0.1870.7830.2970.296φ270.30.3260.300−0.0230.6680.2280.229φ280.30.3420.341−0.0380.7220.2340.237(ii) Hierarchical model scenarioβ15.05.0095.0064.8755.1630.0910.091β20.0−0.0000.001−0.0580.0370.0380.038β30.0−0.001−0.002−0.0460.0350.0310.031β40.00.004−0.000−0.0420.0480.0290.029β50.00.0020.002−0.0500.0630.0360.036β60.0−0.0010.000−0.0580.0370.0380.038β70.00.0060.005−0.0430.0550.0300.030β80.00.0010.003−0.0440.0430.0290.029σ1.01.0131.0090.9271.1060.0600.061φ0−1.0−1.089−1.085−1.348−0.8610.1500.174φ10.30.3520.3480.1490.5150.1180.129φ210.30.2890.323−0.4831.0310.4690.467φ220.30.3480.338−0.2120.7890.2810.283φ230.30.3690.362−0.0620.7350.2460.254φ240.30.2830.256−0.1760.7310.2840.283φ250.30.3530.343−0.0850.7690.2600.264φ260.30.3580.326−0.1290.7700.2620.267φ270.30.3640.352−0.1470.7590.2800.286φ280.30.2920.270−0.1170.7390.2610.260Table 4
The variable selection proportions results forβ in the first simulation study.
Caseβ1β2β3β4β5β6β7β8(i) Empirical Bayesian scenarioi0.000.000.960.970.000.950.940.98ii0.000.000.000.000.000.000.000.00iii0.000.960.980.930.940.960.970.94(ii) Hierarchical model scenarioi0.000.000.990.980.001.001.001.00ii0.000.000.000.000.000.000.000.00iii0.001.000.991.001.001.001.001.00


## 4.2. Example 2

In order to investigate the sensitivity to proportion of missing data, we conduct the second simulation study in this part. Simulated data sets are generated by using the same setup as specified in the first simulation study except that true value ofφ is set as φ0=−2.5, φ1=φ21=⋯=φ28=0.3. We fit these simulated data sets under cases i–iii via the preceding developed MCMC algorithm and only present results in Table 5 on the basis of case i under 100 replications to save space, where the average proportion of missing data generated in the second simulation study is about 0.208. Besides, we further calculate sums of RMSE under case i based on these two simulation studies, and the results are 3.365 and 3.237 in this simulation study, compared with 3.681 and 3.522 in the first simulation study. Moreover, the median of mean absolute deviations (MMAD) is calculated for 100 replications, i.e., median 1/n∑i=1nxi′β^−xi′βtrue, in which β^ and βtrue represent the Bayesian estimate and true value of parameter β.Table 5
Bayesian estimates for casei in the second simulation study.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9882.9972.8463.0960.0830.083β21.51.4821.4731.3281.6380.0890.091β30.0−0.005−0.001−0.0980.1180.0510.051β40.00.0060.000−0.1080.0950.0630.063β52.01.9921.9921.8582.1240.0920.092β60.0−0.008−0.001−0.1580.0560.0640.064β70.0−0.001−0.000−0.0950.0790.0560.055β80.00.0020.000−0.0980.1180.0510.051σ1.00.9890.9860.9141.0790.0580.058φ0−2.5−2.564−2.567−2.689−2.4520.0720.096φ10.30.3390.3450.1040.5250.1350.140φ210.30.2820.261−0.3950.8500.4130.411φ220.30.2740.258−0.2610.8010.3550.354φ230.30.3000.287−0.2050.7100.2840.282φ240.30.3310.329−0.0680.7310.2660.267φ250.30.2830.279−0.2690.7620.3270.326φ260.30.2930.294−0.2550.6980.3040.302φ270.30.2910.311−0.1960.7390.2820.281φ280.30.3870.337−0.0150.8530.2860.298(ii) Hierarchical model scenarioβ13.02.9842.9832.8223.1290.1010.102β21.51.5061.5101.3451.6540.1010.101β30.00.0080.007−0.0670.0890.0470.047β40.0−0.009−0.006−0.0820.0410.0430.043β52.01.9911.9761.8482.1390.1010.101β60.0−0.0010.003−0.0800.0480.0420.042β70.00.001−0.000−0.0750.0770.0460.046β80.00.0070.004−0.0540.0640.0410.042σ1.00.9980.9950.9211.0860.0550.054φ0−2.5−2.573−2.576−2.681−2.4700.0700.101φ10.30.3270.3150.1330.5010.1190.121φ210.30.2870.276−0.3400.7560.3340.333φ220.30.2680.282−0.2600.7100.3200.320φ230.30.3420.314−0.1260.8340.2930.294φ240.30.3070.321−0.2070.7970.3080.307φ250.30.2500.221−0.2510.7350.2810.284φ260.30.3400.354−0.1250.7120.2830.285φ270.30.3230.286−0.1870.8470.3360.335φ280.30.2970.311−0.1850.7070.2810.279The related results of MMAD for the first (second) simulation are 0.143 7 (0.129 3) and 0.134 8 (0.122 1) based on Empirical Bayesian method and Hierarchical model method, respectively. Meanwhile, the box plots for mean absolute deviations (MAD) are shown in Figure1, where “1” and “2” denote case i in the first and second simulation studies, respectively. All these results indicate that cases with the smaller proportion of missing data in this simulation study generally perform better than cases with the larger proportion of missing data in the first simulation study.Figure 1
Box plots summarizing the MAD for cases under (a) empirical Bayesian scenario (left panel) and (b) hierarchical model scenario (right panel) in the first and second simulation studies.
(a)(b)


## 4.3. Example 3

To address the influence of the missingness mechanism model to the Bayesian adaptive lasso, we conduct the third simulation in this example. The full data set is created by employing the same setting specified under case i in the first simulation study. But, we generate the missing data via the following 3 types of missingness mechanisms.Type i: logistic regression model that is different fromM0 specified in (20):(21)logitprri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+φ22xi2+φ23xi3,with φ0=−1.0 and φ1=φ21=φ22=φ23=0.3.Type ii: missing at random (MAR) missingness mechanism.Type iii: nonignorable missingness mechanism model specified inM0.For missing data simulated via type i and type ii, we obtain the Bayesian estimates under incorrect logistic regression modelM0. For missing data simulated via type iii, we calculate the Bayesian estimates under incorrect MAR missingness mechanism. Summaries of these results are reported in Tables 6–9. Inspecting these results, the following observations are indicated: (1) results produced using the model M0 are quite robust and accurate even if the true missing data are MAR, or the true missingness mechanisms are much simpler models; (2) performances obtained on the basis of the incorrect MAR assumption are poor. These observations agreed with those obtained by Lee and Tang [24].Table 6
Bayesian estimates for typei in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.03.0052.9842.8423.1920.1110.111β21.51.5021.5091.3051.6390.1110.111β30.00.0040.000−0.1000.1200.0770.077β40.0−0.004−0.001−0.1070.0720.0610.061β52.01.9781.9881.8152.1290.1010.103β60.00.0140.001−0.0390.0910.0540.056β70.00.0020.000−0.1000.1120.0570.057β80.00.0050.001−0.0620.0970.0580.058σ1.00.9940.9890.8871.0810.0590.059φ0−1.0−1.087−1.089−1.383−0.8190.1740.193φ10.30.3660.3840.0870.5930.1590.171φ210.30.3360.324−0.4100.9550.4410.440φ220.30.2580.292−0.4480.7590.3720.371φ230.30.3570.3290.0380.6670.2110.217φ240.00.027−0.007−0.3800.4370.2600.260φ250.0−0.066−0.098−0.6770.4780.3370.342φ260.00.013−0.015−0.4560.3110.2490.248φ270.0−0.018−0.056−0.4870.4100.2820.281φ280.0−0.015−0.004−0.3930.3480.2470.246(ii) Hierarchical model scenarioβ13.03.0103.0012.8133.1850.1200.120β21.51.4681.4571.2821.6330.1100.114β30.00.0280.018−0.1480.1650.1030.106β40.0−0.0040.002−0.1910.1500.1090.108β52.01.9911.9691.8122.2180.1260.126β60.0−0.003−0.001−0.1780.1420.1020.102β70.0−0.007−0.008−0.1610.1400.1040.103β80.00.002−0.007−0.1570.1650.0940.093σ1.01.0191.0200.9081.1080.0630.066φ0−1.0−1.052−1.055−1.290−0.8490.1430.151φ10.30.3190.3000.0310.5810.1750.175φ210.30.3790.351−0.4941.0810.5030.507φ220.30.3710.364−0.2290.9070.3380.344φ230.30.3430.314−0.0820.7090.2370.240φ240.00.0210.012−0.3810.4260.2550.255φ250.00.0200.006−0.6010.5660.3720.371φ260.00.0320.038−0.4970.3990.2720.271φ270.0−0.054−0.071−0.5120.3510.2700.275φ280.00.0740.065−0.2850.4130.2170.229Table 7
Bayesian estimates for type ii in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9702.9762.7453.1180.1100.113β21.51.4881.4771.3241.6390.1090.110β30.0−0.003−0.001−0.1290.0840.0670.066β40.00.008−0.000−0.0390.0860.0560.057β52.01.9681.9741.7932.1000.0940.099β60.00.0040.000−0.0830.0740.1100.113β70.00.003−0.000−0.0810.0870.0500.050β80.00.0050.000−0.0840.0950.0490.049σ1.01.0010.9960.8941.1060.0640.064φ0—−0.622−0.613−0.888−0.3780.160—φ1—−0.052−0.057−0.3170.1680.151—φ21—0.1990.194−0.6020.8940.476—φ22—0.1090.103−0.4120.5980.314—φ23—−0.013−0.022−0.3830.3110.223—φ24—−0.012−0.001−0.3900.3040.217—φ25—0.1070.098−0.4930.6460.335—φ26—0.0150.026−0.3330.2650.120—φ27—0.006−0.010−0.3810.3370.215—φ28—0.0240.013−0.2860.3640.204—(ii) Hierarchical model scenarioβ13.02.9742.9692.8163.1070.0950.098β21.51.4791.4861.2611.6560.1230.124β30.00.0040.000−0.0660.0700.0470.047β40.00.0070.001−0.0660.0990.0560.056β52.01.9731.9681.8062.1180.1010.104β60.00.0080.003−0.0680.0980.0540.055β70.00.0090.004−0.0510.0470.0490.050β80.0−0.003−0.003−0.0690.0400.0490.049σ1.01.0141.0110.9031.1150.0660.068φ0—−0.606−0.603−0.842−0.3740.152—φ1—−0.054−0.078−0.3610.1680.152—φ21—0.1680.135−0.5520.8630.459—φ22—0.0970.075−0.4290.6010.306—φ23—0.0120.037−0.3420.2960.207—φ24—−0.006−0.029−0.3160.2320.176—φ25—0.1350.157−0.5220.6970.344—φ26—−0.037−0.049−0.3750.2990.200—φ27—0.0500.050−0.3480.3930.234—φ28—−0.0330.029−0.3970.2670.202—Table 8
Bayesian estimates for type iii in the example 3.
ParameterTrueMeanMedian5%95%SDRMSE(i) Empirical Bayesian scenarioβ13.02.9002.9820.0013.6060.7370.740β21.51.3201.4370.0001.7990.6440.666β30.00.1030.003−0.0850.6170.2870.304β40.00.0550.000−0.0820.1810.2400.245β52.01.7451.9640.0012.1280.6400.686β60.00.0610.000−0.0670.4640.2100.218β70.00.007−0.000−0.0130.0190.0420.043β80.0−0.003−0.000−0.0520.0210.0350.035σ1.01.2701.0330.8892.3050.5320.594(ii) Hierarchical model scenarioβ13.03.1953.1852.9743.3950.1400.240β21.51.6291.6061.4141.8810.1430.192β30.00.0500.034−0.1270.3000.1310.139β40.00.0680.048−0.1640.2990.1260.142β52.02.1392.1351.9002.3540.1370.195β60.00.0530.040−0.1570.2700.1230.133β70.00.0520.048−0.1470.2450.1280.137β80.00.0630.058−0.0970.2170.0990.117σ1.01.2501.2421.1451.3640.0750.261


## 4.4. A Real Example

In order to further demonstrate our proposed Bayesian methodologies, we consider a real example from the air pollution data in this subsection, which is available at the websitehttp://lib.stat.cmu.edu/datasets/NO2.data. The air pollution data includes n=500 randomly chosen responses and p=7 explanatory variables. For i=1,2,…,n, the response variable yi is hourly values of the logarithm of the concentration of NO2, measured at Alnabru in Oslo, Norway, between October 2001 and August 2003; the explanatory variable xi1 denotes the logarithm of the number of cars per hour, xi2 represents temperature two meters above ground, xi3 is wind speed, xi4 denotes the temperature difference between 25 and 2 m above ground, xi5 represents wind direction, and xi6 and xi7 are hours of day and day number from October 1st,2001, respectively. Each explanatory variable is standardized to have mean 0 and variance 1 during exploratory data analysis. The air pollution data has been analyzed and fitted by Zhu and Zhu [26], and they observed that all explanatory variables except xi5 and xi7 make a significant contribution to yi. Model equation (19) and the following model(22)logitprri=1|yi,xi,φ=φ0+φ1yi+φ21xi1+⋯+φ27xi7,where φ0=−1.0, φ1=φ21=⋯=φ27=0.5 are used to reanalyze the air pollution data in the article, where the logistic regression model (22) may represent that we create artificial missing data by deliberated deleting some of responses according to the given missing mechanism, and the missing data’s proportion generated in this way is 0.326. Then, we employ the preceding propose Bayesian procedure to make statistical inferences on the assumption of the above model settings. Here, 10 000 samples after 10 000 burn-in iterations are employed to calculate the Bayesian estimates of unknown parameters. In M-H algorithm, we take σy2=2.2 and σφ2=1.0 giving average acceptance rates 29.0% and 26.8%, respectively. Also, we conduct a real data analysis under the assumption for the fully observed data in the framework of Bayesian adaptive Lasso.Table9 reports the corresponding results for regression coefficients by two different methods, in which ‘EB’ and ‘HM’ denote the Bayesian estimates on the basis of Empirical Bayesian method and Hierarchical model method, respectively. From Table 9, we clearly see that (1) the estimate results for the regression coefficients in all two methods are quite close; (2) almost all standard errors obtained under the fully observed data are smaller than those obtained under the missing data, which indicates that the procedures under the fully observed data outperform the procedures under the missing data, but the differences for performances are minor; (3) all these procedures could identify that non-zero significant regression coefficients are β1,β2,β3,β4, and β6, which is consistent with the results given in Zhu and Zhu [26]. In the end, PPP p-value for goodness-of-fit with respect to model (19) with the nonignorable missing mechanism defined in model (20) is 0.446 and 0.454 based on Empirical Bayesian method and Hierarchical model method, respectively, which illustrates that our considered models as well as the proposed goodness-of-fit statistic are feasible.Table 9
Bayesian estimates and standard errors (in the parentheses) of parameterβ in the air pollution data analysis.
Missing dataFully observed dataEBHMEBHMβ10.547 (0.038)0.552 (0.039)0.541 (0.031)0.540 (0.031)β2−0.156 (0.035)−0.163 (0.036)−0.146 (0.028)−0.145 (0.029)β3−0.273 (0.035)−0.268 (0.036)−0.225 (0.025)−0.225 (0.025)β40.179 (0.038)0.200 (0.041)0.172 (0.028)0.171 (0.028)β50.034 (0.028)0.046 (0.035)0.034 (0.025)0.033 (0.028)β6−0.121 (0.038)−0.121 (0.040)−0.123 (0.030)−0.122 (0.031)β70.000 (0.002)0.053 (0.035)0.039 (0.025)0.037 (0.027)


## 5. Conclusion and Discussion

We develop a Bayesian adaptive lasso approach to analyze linear regression models with nonignorable missing responses in this article, where the nonignorable missingness mechanism is specified by a logistic regression model. With the help of Gibbs sampler [14] and the Metropolis-Hastings algorithm [15, 16], we deal with high dimensional integral with not having a closed form in the posterior analysis that is often computationally intractable and then obtain the Bayesian estimates of unknown parameters and their corresponding standard error estimates, the PPP p value [23] for goodness-of-fit. In short, empirical performances of proposed Bayesian approaches in our simulation studies are feasible and effective, and applications for the newly developed statistical tools to the proposed Bayesian adaptive lasso with nonignorable missing responses are plausible.



---

*Source: 3168735-2022-02-11.xml*