Bayesian Elastic Net based on Empirical Likelihood

Regularization methods have been introduced to a linear model to improve the prediction accuracy and interpretability by introducing a penalty term to the least squares criterion; lasso that implements $l_{1}$ penalty [1], ridge that adds $l_{2}$ penalty [2], elastic net (EN) that introduces a combination of $l_{1}$ and $l_{2}$ penalties [3], and non-convex penalties [4, 5]. It is known that EN shows better performances compared to the lasso while preserving a sparse variable selection property [3, 6]. EN also identifies influential variables better than the lasso and has lower false positives than ridge regression [7].

The penalized regression approaches have an intuitive Bayesian counterpart stemming from introducing penalty terms in the form of priors. For instance, Bayesian lasso (BL) [8] uses a double-exponential prior, and Bayesian elastic net (BEN) [9] uses an informative prior that is a compromise between normal and Laplace priors. The non-convex penalties also have been implemented including spike-and-slab prior [10], horseshoe prior [11], and hyperlasso [12]. Motivated by Li and Lin [9], various BEN studies have been conducted including applications to gene-expression data [13], empirical BEN [14], and regressions [15, 16, 17].

The Bayesian approaches on regularization regression have several advantages. First, the Bayesian methods offer a natural interpretation of parameter uncertainty. For example, the standard errors of our parameter of interest can be estimated as the posterior standard deviations because the Bayesian method provides the entire posterior distribution. Second, the Bayesian methods provide valid standard errors compared to the variances of the non-Bayesian penalized regressions estimated by the sandwich or bootstrap [18]. Third, the penalty parameters in the Bayesian methods can be estimated simultaneously with model parameters. When the non-Bayesian models have multiple penalty parameters, such as EN, multi-stage cross-validations are used to estimate the penalty terms. However, the sequential estimation of penalties may incur the over-shrinkage of parameters that leads to larger bias [9]. Lastly, the Bayesian penalized regression approaches show similar or better performances compared to the corresponding non-Bayesian approaches [19, 18, 9].

Owen [20, 21] introduces empirical likelihood (EL), a nonparametric likelihood method to make inference for statistical functionals. EL does not require an assumption that data follow a certain distribution while maintaining similar properties of a conventional likelihood method such as Wilk’s theorem and Bartlett correctability [21, 22, 23]. The EL approach has been applied to complex inferences, including general estimating equation [24], density estimation [25], area under the receiver operating characteristic curve [26], and data imputation [27]. We refer to Chen and Keilegom [28] for a review of EL methods on regression.

EL has extended to high dimensional data. Hjort et al. [29] and Chen et al. [30] show that the asymptotic normality of the EL ratio holds when the dimension of data grows. Tang and Leng [31] and Leng and Tang [32] propose the penalized EL and show that it has the oracle property [4, 33]. Lahiri and Mukhopadhyay [34] propose the penalized EL for estimating population means when the dimension of data is larger than the number of observations. Chang et al. [35] also suggest implementing two penalty functions to penalized models and Lagrange multipliers.

Also, the EL-based Bayesian methods have been developed. Lazar [36] replaces the likelihood function in Bayesian settings by EL and shows the validity of the resulting posterior distribution. In addition, Grendar and Judge [37] show that Bayesian empirical likelihood (BEL) is consistent under the misspecification of the model. BEL has been studied in various areas; the higher-order asymptotic and coverage properties of the posterior distribution for the population mean [38, 39], the survey sampling [40], small area estimation [41], quantile regression [42], Bayesian computation [43], sampling method [44], and lasso and ridge regressions [45].

In this paper, we suggest a Bayesian approach for EN where the likelihood function is replaced by the profile EL of the linear model. We place a special prior distribution on the parameter of interest, which combines the $l_{1}$ and $l_{2}$ penalties leading to the EN approach. The proposed approach takes advantage of the interpretability and robustness of the results achieved by the Bayesian perspective and EL method, respectively. We implement the HMC sampling for BEL suggested by Chaudhuri et al. [44] and propose the leapfrog step size tuning algorithm based on the bisection method. The proposed algorithm enables more efficient sampling than hand-tuning or grid-searching HMC parameters. In addition, our method extends the BEL method for penalized regressions of Bedoui and Lazar [45] by proposing efficient HMC sampling rather than utilizing the tailored Metropolis-Hasting of Chib and Greenberg [46].

The outline of the remaining sections is as follows. In Section 2, we briefly describe the Bayesian linear model based on EL. In Section 3, we propose BEN based on EL and discuss the HMC sampling implementations, along with variable selection procedures. In Section 4, we compare the proposed method with other penalized regression methods using various simulation studies. The air pollution data application is presented in Section 5. Section 6 contains conclusion and future work.

We first outline how EL is implemented in a linear model. Let $X=[\boldsymbol{x_{1}},\cdots,\boldsymbol{x_{p}}]$ be a collection of independent $p$ covariate vectors of size $n$. We consider the linear model

where the error $\epsilon_{i}$ is independent and identically distributed and follows an unknown distribution with mean zero and variance $\sigma^{2}$. The empirical likelihood function for $\boldsymbol{\theta}$ is defined as

and $\boldsymbol{w}=(w_{1},\cdots,w_{n})^{T}$ is the weights on the data points. The coefficients $\boldsymbol{\theta}$ can be estimated by maximizing equation (1) using the Lagrange multipliers. The profile empirical log-likelihood becomes

where $\boldsymbol{\gamma}=\boldsymbol{\gamma}(\boldsymbol{\theta})$ solves the equation

Here, numerical approaches such as the Newton-Raphson method can be used to find the Lagrange multipliers $\boldsymbol{\gamma}$ [47].

The regularization method can be implemented to EL for linear models in different ways. First, the penalty terms can be introduced in the model through the priors $\pi(\boldsymbol{\theta})\propto p_{n}(\boldsymbol{\theta})$, where $p_{n}$ is a penalty function. Second, Tang and Leng [31] introduce the penalty terms in EL $l(\boldsymbol{\theta})\propto-\sum_{i=1}^{n}\log\left\{1+\boldsymbol{\gamma}^{T}\boldsymbol{x_{i}}(y_{i}-\boldsymbol{x_{i}}^{T}\boldsymbol{\theta})\right\}-n\sum_{j=1}^{p}p_{n}(\theta_{j})$. In our study, the $l_{1}$ and $l_{2}$ penalties are introduced in the form of priors.

EN uses both the $l_{1}$ and $l_{2}$ penalties for $p_{t}(\boldsymbol{\theta})$[3]. The EN estimator $\hat{\boldsymbol{\theta}}_{EN}$ is defined as the minimizer of

where

Here, $\lambda_{1}$ and $\lambda_{2}$ are $l_{1}$ and $l_{2}$ penalty parameters that control the amount of shrinkage. Without loss of generality, we assume that

Penalized linear regression models have their Bayesian counterpart models. The Bayesian penalized models introduce priors for $\boldsymbol{\theta}$ where the hyperparameters are functions of the penalty terms. For example, from the form of the EN penalty term in (3), the EN regression parameters can be presented by the following prior

Li and Lin [9] propose the Bayesian counterpart of EN by placing normal and Laplace priors to $l_{1}$ and $l_{2}$ penalties, respectively. The shrinkage parameters, $\lambda_{1}$ and $\lambda_{2}$, are introduced into the model in the form of hyperparameters.

The BEL approach for linear models replaces the likelihood with EL and places parametric priors. Then, the posterior density of EN under the BEL approach $\pi(\boldsymbol{\theta}|\;\mathit{X},\boldsymbol{y})$ becomes

The hierarchical representation of the full model becomes

where $L_{EL}(\boldsymbol{\theta})$ is the profile empirical likelihood for the linear model in (2), IG is an inverse gamma whose probability density function is $\dfrac{b^{a}}{\Gamma\left(a\right)}x^{-a-1}\exp\left(-\dfrac{b}{x}\right),$ for $x>0$, and $X$ and $\boldsymbol{y}$ are omitted in $L_{EL}\left(\boldsymbol{\theta}\right)$ for simplicity. We condition $\boldsymbol{\theta}$ on $\sigma^{2}$ to guarantee the unimodality of the posterior distributions [8, 9]. One can also choose a noninformative prior on $\sigma^{2}$ of the form of $1/\sigma^{2}$.

The prior of $\boldsymbol{\theta}$ given $\sigma^{2}$ defined in (4) is a multiplication of normal and Laplace densities. Andrews and Mallows [48] show that the Laplace distribution can be represented as a scale mixture of normals with an exponential mixing density

Li and Lin [9] also showed that the prior defined in (4) can be presented as a scale mixture of normals with a truncated gamma mixing density. The hierarchical scheme in (4) becomes as follows

where $\boldsymbol{\tau}$ follows the truncated gamma distribution from 1 to $\infty$ with the shape parameter $1/2$ and the rate parameter $\lambda_{1}^{2}/\lambda_{2}\sigma^{2}$. The full joint posterior density is as follows

To sample from the posterior $\pi\left(\boldsymbol{\theta},\;\boldsymbol{\tau},\,\sigma^{2}|X,\;\boldsymbol{y}\right)$, we draw from the following conditional distributions

1. 1.

   Sample $\boldsymbol{\theta}$ from

   $$L_{EL}\left(\boldsymbol{\theta}\right)\exp\left(-\dfrac{\lambda_{2}}{2\sigma^{2}}\sum_{j=1}^{p}\dfrac{\tau_{j}}{\tau_{j}-1}\theta_{j}^{2}\right).$$

2. 2.

   Sampling $\tau_{j}$ is equivalent to sampling $\tau_{j}-1$ from

   $$\text{GIG}\left(\nu=\dfrac{1}{2},\;\psi=\;\dfrac{\lambda_{1}^{2}}{4\lambda_{2}\sigma^{2}},\;\chi=\dfrac{\lambda_{2}}{\sigma^{2}}\theta_{j}^{2}\right)\text{ for }j=1,\cdots,p.$$

3. 3.

   Sample $\sigma^{2}$ from

   $$\text{IG}\left(a+p,\;b+\dfrac{1}{2}\sum_{j=1}^{p}\left[\lambda_{2}\dfrac{\tau_{j}}{\tau_{j}-1}\theta_{j}^{2}+\dfrac{\lambda_{1}^{2}}{4\lambda_{2}}\tau_{j}\right]\right).$$

Here, $\text{GIG}(\nu,\;\psi,\;\chi)$ is the generalized inverse Gaussian distribution with probability density function $\dfrac{\left(\psi/\chi\right)^{\nu/2}}{2K_{\nu}\left(\sqrt{\psi\chi}\right)}x^{\nu-1}\exp\left\{-\dfrac{1}{2}\left(\chi x^{-1}+\psi x\right)\right\}$ for $x>0$, where $K_{\nu}\left(\cdot\right)$ is the modified Bessel function of the third kind with order $\nu$. We use the rejection algorithm proposed by Hörmann and Leydold [49] to generate samples from the GIG distribution.

The posterior estimates of BEN based on EL estimates are consistent. As the sample size increases, the estimates converge to the true value of the parameter being estimated. The consistency of the estimators under the BEL framework has been proved for the quantile regression [42] and the penalized regression [45], and it can be easily extended to the proposed estimates.

We conduct Monte Carlo simulations to compare the performance of the proposed Bayesian elastic net based on the empirical likelihood (BEN-EL) model with four different models: BEN, BL, EN, and lasso applied to the least absolute deviation regression (LADL). LADL is known to be robust to heavy tail errors and outliers [59]. We generate data from the linear model

Note that BEN-EL and LADL assume that the mean and the median of errors are zero, respectively. On the other hand, BEN, BL, and EN assume normal errors.

We use the air pollution to mortality data [62]. The data set includes the mortality rate ($\boldsymbol{y}$) and 15 explanatory variables of weather, socioeconomics, and air pollution in 60 metropolitan areas of the United States from 1959 to 1961. See Table A in Appendix A for a description of the variables.

The air pollution set has a complicated correlation structure. The correlation plot is presented in Figure 2. For example, the nitroc oxide level $nox$ is highly correlated with other air pollution variable hydrocarbon pollution $hc$, but also correlated with the weather variable $prec$ and population variable $popn$.

We randomly split the data into a training set with 30 observations and a test set of 30 observations. The five models, BEN-EL, BEN, BL, EN, and LADL, are conducted using the training set, and the prediction errors are obtained using the test set. For Bayesian methods, the EB approach is used to select the penalty parameters, and the Bayesian credible interval with the probability level $\alpha=0.5$ is used to select the variables. We repeat the procedure 100 times.

Table 5 reports the mean prediction error and standard deviation of the five models. BEN-EL yields the smallest mean prediction error, but its standard deviation is the second largest. The relatively large standard deviation of BEN-EL is due to one large prediction error. BEN-EL has the smallest median, mean, the first quartile (Q1), and the third quartile (Q3) than the other methods. Figure 7 in Appendix A shows boxplots of prediction errors of BEN-EL, BEN, BL, LADL. The results imply that the proposed method outperforms the other methods when data have a complex correlation structure and a small number of observations.

We propose a new Bayesian approach for EN based on EL. Accordingly, the profile EL for the linear model replaces the ordinary likelihood function in the Bayesian setting. We use the HMC algorithm because the resulting posterior EL distribution lacks a closed form, and the implementation of even standard MCMC methods is challenging. A new HMC parameter tuning algorithm based on the bisection method is proposed for BEL. Simulation studies and the air pollution data application show that our approach outperforms the other methods when the sample size is small, data are correlated, and error distributions are misspecified.

BEL methods share the same limitations of EL in the sense that the number of variables cannot increase as the sample size decreases. The reason is that the estimating equations are included in the convex hull to maximize the profile EL ratio. However, $X^{T}X$ is not of full rank when $p>n$ and hence it is non-invertible. To address this shortcoming, several approaches have been proposed, such as penalized EL [34] and double penalization [35]. Thus, it might be worth a formal investigation whether the proposed method could be extended for $p>n$ case for future work.

Another interesting topic is implementing the NUTS algorithm of [51] to the BEL framework. NUTS automatically tunes leapfrog steps and step sizes and has been used as a default sampler in popular MCMC software such as Stan [63] and PyMC3 [64]. However, the proper introduction of NUTS into BEL has not been studied. We observe that HMC chains often diverge when the NUTS algorithm is used for BEN-EL. Also, BEL is not supported by most MCMC libraries because EL cannot be used as a likelihood function. Developing a BEL computational pipeline compatible with Stan or PyMC3 for practitioners will be of further interest.

No potential conflict of interest was reported by the authors.

All data used in simulation studies are generated randomly and the air pollution data are obtained from McDonald et al. [62]. The R code used to generate, import, and analyze the data used in this paper is publicly available at the URL: https://github.com/chulmoon/BEN-EL.