A novel Bayesian approach for variable selection in linear regression models

TL;DR

This paper introduces a new Bayesian method for variable selection in linear regression that incorporates prior beliefs, uses a stable g-prior, and employs an advanced MCMC algorithm, demonstrating competitive performance and consistency.

Contribution

It presents a hierarchical Bayesian framework with a novel MCMC algorithm for variable selection, allowing direct prior specification and ensuring numerical stability.

Findings

Performs at least as well as existing methods on real and simulated data.

Proves consistency of the variable selection approach under certain conditions.

Introduces a stable g-prior with ridge parameter for improved numerical stability.

Abstract

We propose a novel Bayesian approach to the problem of variable selection in multiple linear regression models. In particular, we present a hierarchical setting which allows for direct specification of a-priori beliefs about the number of nonzero regression coefficients as well as a specification of beliefs that given coefficients are nonzero. To guarantee numerical stability, we adopt a $g$-prior with an additional ridge parameter for the unknown regression coefficients. In order to simulate from the joint posterior distribution an intelligent random walk Metropolis-Hastings algorithm which is able to switch between different models is proposed. Testing our algorithm on real and simulated data illustrates that it performs at least on par and often even better than other well-established methods. Finally, we prove that under some nominal assumptions, the presented approach is consistent in terms of model selection.

Tables

Table 1. Table 1: Diabetes data performance comparison: Median of mean squared prediction errors (MMSE) and median of mean absolute prediction deviations (MMAD) based on a 5 5 5 -fold cross-validation.

Method	MMSE	MMAD
New approach	0.4873534	0.5678801
Lasso	0.492067	0.571596
Adaptive lasso	0.4939229	0.5736721
Elastic net	0.4922686	0.5706994
Bayesian lasso	0.4924316	0.5736084
Bayesian adaptive lasso	0.4997307	0.5786672
Bayesian elastic net	0.4895844	0.5727555

Table 2. Table 2: Prostate cancer data performance comparison: Median of mean squared prediction errors (MMSE) and median of mean absolute prediction deviations (MMAD) based on a 5 5 5 -fold cross-validation.

Method	MMSE	MMAD
New approach	0.5068169	0.6089882
Lasso	0.5214503	0.6023388
Adaptive lasso	0.5429631	0.5990569
Elastic net	0.5293173	0.5837592
Bayesian lasso	0.5462163	0.5916963
Bayesian adaptive lasso	0.547966	0.5614767
Bayesian elastic net	0.5267334	0.6138624

Table 3. Table 3: Ozone data performance comparison: Median of mean squared prediction errors (MMSE) and median of mean absolute prediction deviations (MMAD) based on a 5 5 5 -fold cross-validation.

Method	MMSE	MMAD
New approach	0.2471677	0.3837716
Lasso	0.2451734	0.3898133
Adaptive lasso	0.2539498	0.3893223
Elastic net	0.2498921	0.3898747
Bayesian lasso	0.256831	0.3926603
Bayesian adaptive lasso	0.2458003	0.3865291
Bayesian elastic net	0.2409481	0.3853373

Table 4. Table 4: Performance comparison for different specifications of n 𝑛 n : Median of mean squared prediction errors (MMSE) and median of mean absolute prediction deviations (MMAD) based on a 100 100 100 simulated datasets.

Method	MMSE	MMAD	MMSE	MMAD	MMSE	MMAD
	$n=50$	$n=50$	$n=100$	$n=100$	$n=200$	$n=200$
New approach	0.3390823	0.4618555	0.2436617	0.3933677	0.234171	0.3868592
Lasso	0.4128203	0.5142703	0.2972025	0.4385642	0.2612206	0.4064837
Adaptive lasso	0.4146856	0.5174635	0.2796941	0.422432	0.2606342	0.4076061
Elastic net	0.4557527	0.5405664	0.3204951	0.4521248	0.2667704	0.4108516
Bayesian lasso	0.7189461	0.6758523	0.3032926	0.439959	0.2596804	0.4053971
Bayesian adaptive lasso	0.7795937	0.7110935	0.8896619	0.7637107	0.3793791	0.4924849
Bayesian elastic net	0.8148157	0.7146019	0.3324611	0.4647419	0.2668035	0.4122355

Table 5. Table 5: Performance comparison for different specifications of n 𝑛 n : Median of mean squared prediction errors (MMSE) and median of mean absolute prediction deviations (MMAD) based on a 100 100 100 simulated datasets.

Method	MMSE	MMAD	MMSE	MMAD	MMSE	MMAD
	$n=50$	$n=50$	$n=100$	$n=100$	$n=200$	$n=200$
New approach	0.9155098	0.7634292	0.4057127	0.5094444	0.343676	0.4678344
Lasso	0.8979915	0.7481239	0.4926306	0.5699708	0.4065196	0.5115676
Adaptive lasso	0.9159057	0.7556247	0.5108755	0.5731952	0.4093501	0.5127194
Elastic net	0.8918046	0.7524384	0.5090338	0.574375	0.4080551	0.5138157
Bayesian lasso	1.043241	0.8165255	0.5140319	0.572721	0.4343926	0.5303905
Bayesian elastic net	1.147992	0.8520899	0.7743223	0.6907885	0.8825549	0.7553813

Equations

$$\mathbf{y}=\mathbf{X}\mathbf{\beta }+\mathbf{\epsilon }$$

$$p(\mathbf{\beta }|{\sigma }^2)=\prod _{i=1}^p\frac{\lambda }{2\sigma }{e}^{-\lambda \frac{|{\beta }_i|}{\sigma }}$$

$$\mathcal{P}(\mathbf{\beta },\mathbf{\lambda }):=\sum _{i=1}^p{\lambda }_i|{\beta }_i|$$

$$\underset{\mathbf{\beta }}{\min }\sum _{i=1}^n(y_i-{\text{x}}_i^T\mathbf{\beta }{)}^2+\lambda ||\mathbf{\beta }|{|}_2^2$$

$$p(\mathbf{\beta }|{\sigma }^2)=\prod _{i=1}^p\frac{{\lambda }_i}{2\sigma }{e}^{-{\lambda }_i\frac{|{\beta }_i|}{\sigma }}$$

$$\mathcal{P}(\mathbf{\beta },\mathbf{\lambda }):={\lambda }_1||\mathbf{\beta }|{|}_1+{\lambda }_2||\mathbf{\beta }|{|}_2^2$$

$$p(\mathbf{\gamma })=\prod _{i=1}^p{\gamma }_i{\pi }_i+(1-{\gamma }_i)(1-{\pi }_i)$$

$$p(\mathcal{A}=\{{\alpha }_1,...,{\alpha }_k\})\propto (p_{{\alpha }_1}+...+p_{{\alpha }_k})\frac{1}{k}\tilde{p}(k)$$

$$p(\mathcal{A}=\{{\alpha }_1,...,{\alpha }_k\})\propto \tilde{p}(k).$$

$$p({\sigma }^2)\propto ({\sigma }^2{)}^{-(a+1)}\exp \left(-\frac{b}{{\sigma }^2}\right)$$

$${\mathbf{\beta }}_{\mathcal{A}}|g,{\sigma }^2,{\mathbf{X}}_{\mathcal{A}}\sim \mathcal{N}(0,(g^{-1}{\sigma }^{-2}{\mathbf{X}}_{\mathcal{A}}^T{\mathbf{X}}_{\mathcal{A}}+\lambda {\text{I}}_k{)}^{-1})$$

$$\lambda =\{\begin{array}{ll}\max (\frac{1}{k},\frac{1}{300})& \text{if }n\le \zeta \\ 0& \text{else}\end{array}$$

$$\mathbf{y}|\mathbf{X},\mathbf{\beta },{\sigma }^2\sim \mathcal{N}(\mathbf{X}\mathbf{\beta },{\sigma }^2{\text{I}}_n)$$

$$p(\mathbf{\beta }|\mathcal{A},g,{\sigma }^2,{\mathbf{X}}_{\mathcal{A}})\propto \exp \left(-\frac{1}{2}{\mathbf{\beta }}_{\mathcal{A}}^T\left(g^{-1}{\sigma }^{-2}{\mathbf{X}}_{\mathcal{A}}^T{\mathbf{X}}_{\mathcal{A}}+\lambda {\text{I}}_k\right){\mathbf{\beta }}_{\mathcal{A}}\right)h({\mathbf{\beta }}_{{\mathcal{A}}^C})$$

$$p(\mathcal{A}=\{{\alpha }_1,...,{\alpha }_k\})\propto (p_{{\alpha }_1}+...+p_{{\alpha }_k})\frac{1}{k}\tilde{p}(k)$$

$$g\sim \text{IG}\left(\frac{1}{2},\frac{n}{2}\right)$$

$${\sigma }^2\sim \text{IG}\left(a,b\right)$$

$$h({\mathbf{\beta }}_{{\mathcal{A}}^C})=\{\begin{array}{ll}1& \text{if }{\mathbf{\beta }}_{{\mathcal{A}}^C}=0\\ 0& \text{else}\end{array}$$

$$\mathbf{y}={\mathbf{X}}_{\mathcal{A}}{\mathbf{\beta }}_{\mathcal{A}}+\mathbf{\epsilon }\text{ with }\mathbf{\epsilon }\sim \mathcal{N}(0,{\sigma }^2{\text{I}}_n).$$

$$\underset{n\to \infty }{\mathrm{p}\lim }p(M_{\mathcal{A}}|\mathbf{y},\mathbf{X})=1\text{ and }\underset{n\to \infty }{\mathrm{p}\lim }p(M_{{\mathcal{A}}^{\prime }}|\mathbf{y},\mathbf{X})=0\text{ for all }{\mathcal{A}}^{\prime }\ne \mathcal{A}$$

$$\underset{n\to \infty }{\mathrm{p}\lim }\frac{{\mathbf{y}}^T({\text{I}}_n-{\mathbf{X}}_{{\mathcal{A}}^{\prime }}{\left({\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T{\mathbf{X}}_{{\mathcal{A}}^{\prime }}\right)}^{-1}{\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T)\mathbf{y}}{n}={\sigma }^2.$$

$$\underset{n\to \infty }{\lim }\frac{{\beta }_{\mathcal{A}}^T{\mathbf{X}}_{\mathcal{A}}^T({\text{I}}_n-{\mathbf{X}}_{{\mathcal{A}}^{\prime }}{\left({\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T{\mathbf{X}}_{{\mathcal{A}}^{\prime }}\right)}^{-1}{\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T){\mathbf{X}}_{\mathcal{A}}{\beta }_{\mathcal{A}}}{n}=b_{{\mathcal{A}}^{\prime }}\in (0,\infty ),$$

$$\underset{n\to \infty }{\mathrm{p}\lim }\frac{{\mathbf{y}}^T({\text{I}}_n-{\mathbf{X}}_{{\mathcal{A}}^{\prime }}{\left({\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T{\mathbf{X}}_{{\mathcal{A}}^{\prime }}\right)}^{-1}{\mathbf{X}}_{{\mathcal{A}}^{\prime }}^T)\mathbf{y}}{n}={\sigma }^2+b_{{\mathcal{A}}^{\prime }}.$$

$$q({\sigma }_{t+1}^2|{\sigma }_t^2)$$

$$q(g_{t+1}|g_t)$$

$$c_h\sim \text{Bernoulli}(p_h)$$

$$q(\alpha |{\mathcal{A}}_t)={\text{I}}_{\{1,...,p\}}(\alpha )\{\begin{array}{ll}{\tilde{p}}_{\alpha }& \text{if }{k}_t>1,\alpha \notin {\mathcal{A}}_t\\ \frac{\sum _{i\in {\mathcal{A}}_t}{\tilde{p}}_i}{{\tilde{p}}_{\alpha }\sum _{i\in {\mathcal{A}}_t}1/{\tilde{p}}_i}& \text{if }{k}_t>1,\alpha \in {\mathcal{A}}_t\\ 0& \text{if }{k}_t=1,\alpha \in {\mathcal{A}}_t\\ \frac{{\tilde{p}}_{\alpha }}{\sum _{i\in \{1,...,p\}\setminus {\mathcal{A}}_t}{\tilde{p}}_i}& \text{else}\end{array}$$

$$q({\mathcal{A}}_{t+1}|{\mathcal{A}}_t,c_h=1)=q([{\mathcal{A}}_{t+1}\setminus {\mathcal{A}}_t]\cup [{\mathcal{A}}_t\setminus {\mathcal{A}}_{t+1}]|{\mathcal{A}}_t){\text{I}}_{\{1\}}(|k_{t+1}-k_t|)$$

$$q({\mathcal{A}}_{t+1}|{\mathcal{A}}_t,c_h=0)=\{\begin{array}{ll}1& \text{if }{\mathcal{A}}_{t+1}={\mathcal{A}}_t\\ 0& \text{else}\end{array}$$

$$q({\mathcal{A}}_{t+1}|{\mathcal{A}}_t)=\sum _{c_h\in \{0,1\}}q({\mathcal{A}}_{t+1},c_h|{\mathcal{A}}_t)$$

$$=\sum _{c_h\in \{0,1\}}q({\mathcal{A}}_{t+1}|c_h,{\mathcal{A}}_t)q(c_h)$$

$$\mathbf{\beta }|{\sigma }^2,\mathbf{y},\mathbf{X}\sim \mathcal{N}(\tilde{\mathbf{m}},{\sigma }^2\tilde{\mathbf{M}})$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

Nowadays, where computing becomes less expensive, huge amounts of data are routinely collected by companies and organizations of all kinds. The scale of the data being collected increases greatly and thus also the number of measured features (variables/predictors/covariates). Fitting statistical models, which include a large number of features, often results in over-fitting since the resulting model complexity cannot be supported by the available sample size (Buehlmann et al. [2016]). Moreover, it becomes hard to detect the most predictive features and their importance with respect to the model. Thus, variable selection (i.e. the extraction of the relevant features for a given task) plays an ever more important role in many fileds including genetics (Lee et al. [2003]), astronomy (Zheng and Zhang [2007]), and economics (Foster and Stine [2004]).

In this paper we are interested in variable selection in the multiple linear regression model

[TABLE]

where $\mathbf{y}$ is a $n$-dimensional response vector, $\mathbf{X}=(\mathbf{x}_{1},...,\mathbf{x}_{p})\in\mathbb{R}^{n\times p}$ denotes the so-called design matrix which holds the $p$ potential predictors, $\boldsymbol{\beta}=(\beta_{1},...,\beta_{p})^{T}$ is the vector of unknown regression coefficients, and $\boldsymbol{\varepsilon}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{\text{I}}_{n})$ denotes the noise vector with independent and identically distributed components. To avoid the need for an intercept $\beta_{0}$ it is assumed that the response $\mathbf{y}$ is zero centered. Moreover, to allow for an easy evaluation of the influence of single predictors $x_{1},...,x_{p}$ on the model based on the magnitude of the regression coefficients the predictors are treated as zero centered $\sum_{i=1}^{n}x_{ji}=0$ and standardized with $||\mathbf{x}_{j}||\texttwosuperior=n-1$ for $1\leq j\leq p$.

In the linear model $(1)$ one assumes that only some regression coefficients are different from zero. Thus, the problem of variable selection reduces to the identification of the nonzero regression coefficients (Alhamzawi and Taha Mohammad Ali [2018]). Especially shrinkage approaches such as the lasso (Tibshirani [1996]), the adaptive lasso (Zou [2006]), the elastic net (Zou and Hastie [2005]) and their Bayesian analogues (Park and Casella [2008], Alhamzawi et al. [2012], Leng et al. [2014], Huang et al. [2015]) which simultaneously perform variable selection and coefficient estimation have been shown to be effective and are often the methods of choice in linear regression. These methods estimate $\boldsymbol{\beta}$ as minimizer of the objective function $\mathcal{L}(\boldsymbol{\beta})+\mathcal{P}(\boldsymbol{\beta},\boldsymbol{\lambda})$, where $\mathcal{L}$ denotes the quadratic loss function (negative log-likelihood) $\sum_{i=1}^{n}(y_{i}-\textbf{x}_{i}^{T}\boldsymbol{\beta})^{2}$ and $\mathcal{P}$ denotes a method specific penalty function that encourages a sparse solution. The parameter vector $\boldsymbol{\lambda}$ controls the penalization strength and thus the level of sparsity. Commonly, a good choice for $\boldsymbol{\lambda}$ is determined via cross-validation in practical applications. Among above mentioned approaches the most popular one is the least absolute shrinkage and selection operator (lasso), which was proposed by Tibshirani [1996]. In the lasso regression the penalty function is defined as $\mathcal{P}(\boldsymbol{\beta},\lambda):=\lambda||\boldsymbol{\beta}||_{1}=\lambda\sum_{i=1}^{p}|\beta_{i}|$. Thus, the regression coefficients are continuously shrunken to zero and for sufficiently large $\lambda$ some coefficients take exactly the value zero. The lasso can be viewed as a convex and therefore more efficiently solvable reformulation of the best subset selection approach, in which the penalty function is given by $\mathcal{P}(\boldsymbol{\beta},\lambda):=\lambda||\boldsymbol{\beta}||_{0}=\lambda\#(i|\beta_{i}\neq 0)$. Moreover, the lasso estimate for $\boldsymbol{\beta}$ can be interpreted as a Bayesian posterior mode estimate when independent Laplace priors all with zero mean and the same scale parameter $\lambda>0$ are assigned to the regression coefficients:

[TABLE]

In contrast to the frequentist approach the Bayesian one (Park and Casella [2008]) provides credible intervals for the model parameters which can be used for variable selection. A further generalization of the classical lasso is the adaptive lasso (Zou [2006]), which allows for different penalization factors of the regression coefficients:

[TABLE]

For instance, the penalization vector $\boldsymbol{\lambda}=(\lambda_{1},...,\lambda_{p})^{T}$ can be chosen based on the ordinary least squares estimator $\lambda_{j}^{-1}=|\widehat{\boldsymbol{\beta}}_{j}(\text{ols})|^{\gamma}=|((\mathbf{X}^{T}\mathbf{X})^{-1}\mathbf{X}^{T}\mathbf{y})_{j}|^{\gamma}$ with $\gamma>0$, or in more complicated settings (collinear predictors, or $n<p$) similarly based on the minimizer of the ridge-regularized problem given by:

[TABLE]

The adaptive lasso was proposed to address potential weaknesses of the classical lasso. It has been shown (Zou [2006]) that there are situations in which the classical lasso either tends to select inactive predictors, or over-shrinks the regression coefficients of correct predictores. In contrast, the adaptive lasso satisfies the so-called oracle properties introduced by Fan and Li [2001]. A procedure satisfies these properties, if it asymptotically identifies the right subset model and has optimal estimation rate. For more details please refer to Fan and Li [2001]. Recently, a Bayesain version of the adaptive lasso was proposed by Alhamzawi et al. [2012] and also by Leng et al. [2014]. The Bayesian adaptive lasso generalizes the classical (i.e. non adaptive) Bayesian lasso by allowing different scale parameters in the Laplace priors of the regression coefficients:

[TABLE]

The scale parameters $\lambda_{i}$ are then either choosen via marginal maximal likelihood in an emprical Bayesian setting or by assigning appropriate hyperpriors in a hierarchical Bayes approach. Another generalization of the classical lasso is the elastic net proposed by Zou and Hastie [2005]. Here the penalty function is given by:

[TABLE]

The elastic net encourages a grouping of strongly correlated predictors, such that they tend to be in or out the model together. Morover, it works better than the classical lasso in the case where $p\gg n$. There are also some Bayesian versions of the elastic net proposed in literature including the one proposed by Huang et al. [2015]. In particular, we used their R-package EBglmnet to validate our novel method.

Besides above decribed shrinkage approaches, Bayesian methods that are based on the introduction of a random indicator vector $\boldsymbol{\gamma}=(\gamma_{1},...,\gamma_{p})^{T}\in\{0,1\}^{p}$ gain increasing popularity (Liang et al. [2008], Guan and Stephens [2011], Fisher and Mehta [2014], Wang et al. [2015]). If $\gamma_{i}$ equals zero the regression coefficient $\beta_{i}$ of the $i$-th predictor is assumed to equal zero. This is equivalent to the assumption that the $i$-th predictor does not explain the target $y$. A common choice is to assign independent Bernoulli priors to the indicator variables:

[TABLE]

For instance, Guan and Stephens [2011] used prior $(2)$ with the simplification $\pi=\pi_{1}=...=\pi_{p}$ and assigned an uniform Prior $U(a,b)$ to $\operatorname{ln}(\pi)$. In the work of Fisher and Mehta [2014] the flat prior $\boldsymbol{\gamma}\propto 1$ was used which is equivalent to seting $\pi_{i}=0.5$ in $(2)$ for $i=1,...,p$. For given $\boldsymbol{\gamma}$ the vector of non-zero coefficients denoted by $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$ is usally assigned a conventional $g$-prior $\boldsymbol{\beta}_{\boldsymbol{\gamma}}|g,\sigma^{2},\mathbf{X}_{\mathcal{A}}\sim\mathcal{N}(\mathbf{0},g\sigma^{2}(\mathbf{X}_{\boldsymbol{\gamma}}^{T}\mathbf{X}_{\boldsymbol{\gamma}})^{-1})$ (introduced by Zellner [1986]) where $\mathbf{X}_{\boldsymbol{\gamma}}$ denotes the design matrix with all columns deleted that correspond to indicator variables being equal to zero. The reason for the common choice of this prior is that it often leads to a computationally tractable Bayes factor.

Inspired by above described approaches, we propose a setting which is not based on a random indicator vector $\boldsymbol{\gamma}$, but on a random set $\mathcal{A}\subseteq\{1,...,p\}$ that holds the indices of the active predictors, i.e. the predictors with coefficients different from zero. We assign a prior to $\mathcal{A}$ which depends on the cardinality of the set $|\mathcal{A}|$ as well as on the actual elements of $\mathcal{A}$. This enables an easy formulation of useful a priori knowledge. In section 3 we show how the results of other methods can be used in order to define such a useful prior in terms of an empirical Bayes approach. To simulate from the joint posterior of our model we propose a novel random walk Metropolis-Hastings algorithm, in contrast to Guan and Stephens [2011], Fisher and Mehta [2014], Wang et al. [2015] where Gibbs sampling is applied for this task. This is especially useful when it is difficult or even impossible to determine the conditional distributions needed for the Gibbs sampling algorithm in an analytically closed form. In these cases the corresponding integrals have to be approximated and in particular for practitioners it can become a difficult task to choose a proper approximation technique among the various possible ones. For instance, in Liang et al. [2008] the Laplace approximation is used for this task. Wang et al. [2015] assigned on purpose a beta-prime prior with specific parameters to g in order to achive a closed-form expression of the marginal posterior of $\boldsymbol{\gamma}$. Since we don’t need to care about closed-form expressions by using a random walk Metropolis algorithm, our approach is based on the more popular Zellner-Siow prior for $g$ (Zellner and Siow [1980]), which is an inverse gamma distribution with shape parameter $\alpha=1/2$ and scale parameter $\beta=n/2$.

However, the $g$-prior depends on the inverse of the empirical covariance matrix of the selected predictors. This matrix is singular if the number of selected covariates is greater than the number of observations $n$ and, further, may be almost rank deficient given that the predictors are highly correlated. To overcome this problem Wang et al. [2015] replaced the classical inverse with the Moore-Penrose generalized inverse and thus ended up with the so-called $gsg$-prior (see West [2003]). In contrast to them, we adopt a $g$-prior with an additional ridge parameter for the unknown regression coefficients to guarantee nonsingularity of the empirical covariance matrix. This modification of the classical $g$-prior was first proposed by Gupta and Ibrahim [2007] and further investigated by Baragatti and Pommeret [2012].

Finally, in Section 2.2 we state that our approach is consistent in terms of model selection according to the consistency definition given by Fernández et al. [2001]. The proof of this result is deferred to the appendix. Moreover, in Section 3, we evaluate our approach on the basis of real and simulated data and compare the results with the already described shrinkage methods. We show that our approach performs at least on par and often better than the comparative methods.

2 Methodology

In this section we describe the hierarchical Bayesian approach we propose for the task of variable selection in multiple linear regression models. Moreover, we state that the novel model presented is consistent in terms of model selection. Finally, we propose an intelligent random walk Metropolis-Hastings algorithm to simulate from the joint posterior of the model parameters.

2.1 Hierarchical representation of the full model

In order to perform variable selection, a random set $\mathcal{A}\subseteq\{1,...,p\}$ which holds the indices of the active predictors, i.e. the predictors with coefficients different from zero, is used. Moreover, it is assumed that the cardinality $k=|\mathcal{A}|$ is greater than zero. In this paper the focus lies on problems where at least one of the measured features is predictive. Nevertheless, one can easily generalize the proposed approach in such a way that also the null model (all regression coefficients equal zero) is valid. Thus, we propose the following prior for $\mathcal{A}$:

[TABLE]

where

1. (i)

$\{p_{\alpha_{1}},...,p_{\alpha_{k}}\}\subseteq\{p_{1},...,p_{p}\}$ with $\sum_{i=1}^{p}p_{i}=1$ and $p_{i}\geq 0$ for $i=1,...,p$ 2. (ii)

$\tilde{p}:\{1,...,p\}\rightarrow\mathbb{R}_{0}^{+}$.

The mapping $\tilde{p}$ can be arbitrarily chosen and should represent the a-priori belief in the model size, i.e. the number of parameters different from zero. For instance one could define $\tilde{p}$ based on the probability mass function of a zero truncated binomial distribution with size parameter equal to the number of predictors $p$ and the second parameter chosen in such a way that the mean of the distribution equals the a-priori most probable model size. On the other hand the parameters $p_{1},...,p_{p}$ are used to represent the a-priori belief in the importance of the predictors $x_{1},...,x_{p}$ to the model. In terms of empirical Bayes these parameters could be chosen based on the regression coefficients of a ridge regularized model (compare to the adaptive lasso) or based on the correlations (i.e. the linear dependencies) between the target variable and the predictors. Please note that in the case of equal a-priori importance of the predictors $\frac{1}{p}=p_{1}=...=p_{p}$, the prior $(3)$ reduces to

[TABLE]

Thus, in this case the prior only reflects the a-priori belief regarding $k$, which shows the necessity of the factor $\frac{1}{k}$ in $(3)$.

For the variance $\sigma^{2}$ of the error terms $\varepsilon_{1},...,\varepsilon_{n}$ an inverse gamma prior is chosen:

[TABLE]

For $b\rightarrow 0$ and $a\rightarrow 0$ the inverse gamma prior converges to the non-informative and improper Jeffreys prior $p(\sigma^{2})\propto\frac{1}{\sigma^{2}}$ (Jeffreys [1961]) which is invariant under a change of scale. In the experimental studies in Section 3, $a$ and $b$ are both set to $0.001$ such that the chosen prior is a proper approximation of the Jeffreys prior.

As already stated in the introduction, for given $\mathcal{A}$ the vector of non-zero coefficients denoted by $\boldsymbol{\beta}_{\mathcal{A}}$ is commonly assigned a conventional $g$-prior $\boldsymbol{\beta}_{\mathcal{A}}|g,\sigma^{2},\mathbf{X}_{\mathcal{A}}\sim\mathcal{N}(\mathbf{0},g\sigma^{2}(\mathbf{X}_{\mathcal{A}}^{T}\mathbf{X}_{\mathcal{A}})^{-1})$ (introduced by Zellner [1986]), where $\mathbf{X}_{\mathcal{A}}$ denotes the submatrix of $\mathbf{X}$ consisting of all columns corresponding to predictors with index in $\mathcal{A}$. Since the $g$-prior depends on the inverse of the empirical covariance matrix of the selected predictors it is singular if the number of selected covariates $k$ is greater than the number of observations $n$ and further may be almost rank deficient provided that the predictors are highly correlated (see West [2003], Gupta and Ibrahim [2007], Baragatti and Pommeret [2012]). To overcome these two problems, based on the ideas of Gupta and Ibrahim [2007] and Baragatti and Pommeret [2012], we consider a ridge penalized version of the $g$-prior:

[TABLE]

where

[TABLE]

with $\zeta$ denoting an appropriately large constant. It should be explained what appropriately large means. Obviously, for $n>\zeta=p-1$ the first possible problem ($k>n$) cannot appear. The second problem, i.e. an almost rank deficient matrix $\mathbf{X}_{\mathcal{A}}^{T}\mathbf{X}_{\mathcal{A}}$, also should not appear for $n\gg k$, i.e. for sufficiently large $\zeta$. Even if the predictors are highly correlated, a huge sample size will prevent $\mathbf{X}_{\mathcal{A}}^{T}\mathbf{X}_{\mathcal{A}}$ from beeing nearly singular and thus maybe computationally singular (compare Carsey and Harden [2013]). Although it has to be assumed that there is no predictor that is an exact linear combination of the other predictors. However, this assumption can be easily achieved by excluding such non-informative predictors from the beginning. It should be mentioned that setting $\lambda$ to zero is not necessary at all in practical applications. This is just a theoretical consideration which will ease the proof of model selection consistency in A. The definition of $\lambda$ for $n\leq\zeta$ is taken from Baragatti and Pommeret [2012], who suggested to reduce the influence of the ridge parameter when the number of predictors is large up to the threshold $1/300$.

Finally, we assign the popular Zellner-Siow prior (Zellner and Siow [1980]), i.e. an inverse gamma distribution with shape parameter $\alpha=1/2$ and scale parameter $\beta=n/2$ to $g$, resulting in the complete hierarchical representation of the model

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\boldsymbol{\beta}_{\mathcal{A}^{C}}$ denotes the components of $\boldsymbol{\beta}$ with indices not in $\mathcal{A}$,

[TABLE]

and $\lambda$ satisfies equation $(6)$.

2.2 Model consistency

In this section we state that model $(7-12)$ is consistent in terms of model selection. The proof of the result is is given in A. Therefore, throughout the section it is assumed that the sample $\mathbf{y}$ is generated by model $M_{\mathcal{A}}$ with parameters $\boldsymbol{\beta}_{\mathcal{A}}$ and $\sigma^{2}$, i.e.,

[TABLE]

In rough words, model selection consistency means that the true model will be selected provided that enough data is available. A sound mathematical definition was given by Fernández et al. [2001]. According to these authors the posterior probability is said to be consistent if

[TABLE]

where the probability limit is taken with respect to the true sampling distribution $(13)$. Moreover, some preliminary results which are required in our proof are formulated as Lemma A.1 in their work. In our notation this Lemma reads as:

Under the sampling model (also true model) $M_{\mathcal{A}}$ $(13)$,

If $M_{\mathcal{A}}$ is nested within or equal to model $M_{\mathcal{A}^{\prime}}$,

[TABLE] 2. 2.

Under the assumption that for any model $M_{\mathcal{A}^{\prime}}$ that does not nest $M_{\mathcal{A}}$,

[TABLE]

one obtains

[TABLE]

Theorem 1

Assume that $(16)$ holds. Then the proposed model setup $(7-12)$ is consistent in terms of model selection.

2.3 Metropolis-Hastings algorithm

As usual in Bayesian regression models, the joint posterior $p(\boldsymbol{\beta},\mathcal{A},g,\sigma^{2}|\mathbf{y},\mathbf{X})$ is analytically intractable. The most common approach to overcome this problem is to use Markov chain Monte Carlo (MCMC) methods to simulate from this unknown distribution. Such a method generates realizations of a Markov chain which converges in distribution to the (conditional) random vector of interest, i.e. in this paper to $\boldsymbol{\beta},\mathcal{A},g,\sigma^{2}|\mathbf{y},\mathbf{X}$. Thus, after some time of convergence, the so-called burn-in phase, the generated random numbers can be considered as a (dependent) sample from the posterior distribution. To finally obtain an independent sample only simulated numbers with a given distance $d$ are chosen, a so-called thinning is performed. In this paper a special Metropolis Hastings (MH) algorithm, and thus a special MCMC method, is proposed to simulate from $p(\boldsymbol{\beta},\mathcal{A},g,\sigma^{2}|\mathbf{y},\mathbf{X})$. MH algorithms are iterative and in each iteration a random number is drawn from a proposal distribution $q$ which is then accepted or not based on a given criterion. If the proposal $q$ depends on the respectively last accepted random number one also speaks of a random walk MH algorithm. For further details on MCMC methods please refer to Fahrmeir et al. [2013]. In order to define a MH algorithm for a given problem basically the only thing one has to specify is the proposal distribution $q$ used in the algorithm. Therefore, we define the proposal $q(\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},g_{t+1},\sigma_{t+1}^{2}|\mathcal{A}_{t},g_{t},\sigma_{t}^{2})$ as follows:

Inspired by the work of Hastings [1970] the proposal distribution for the variance parameter $\sigma^{2}$ as well as the proposal distribution of the parameter $g$ are defined to be uniform distributions

[TABLE]

where $\varepsilon_{\sigma}>0$ and $\varepsilon_{g}>0$ denote tuning parameters. Tuning parameters are always chosen in such a way, that the acceptance rate of the algorithm is neither very low nor very high which should guarantee a fast convergence of the algorithm. In order to specify the proposal for the random set $\mathcal{A}$ at first a Bernoulli distributed random variable $c_{h}$ is introduced

[TABLE]

where $p_{h}\in[0,1]$ can be considered as a tuning parameter. The event $c_{h}=1$ means that the model size (number of nonzero regression coefficients) changes, i.e. $k_{t+1}\neq k_{t}$. Moreover, a random variable $\alpha$ with support on the index set $\{1,...,p\}$ of the regression coefficients is introduced in order to describe the model transition probabilities in case of changing model size. A realization of this random variable corresponds to the index of a predictor which is going to be added to or removed from the model. Transitions between models, which differ by two or more parameters, are not allowed. The probability mass function of $\alpha$ is defined as

[TABLE]

where the parameters $\tilde{p}_{1},...,\tilde{p}_{p}$ are greater or equal to zero, sum up to one, and can but must not be identical to the parameters $p_{1},...,p_{p}$ already used in the prior $p(\mathcal{A})$. Again, as the parameters $p_{1},...,p_{p}$ the parameters $\tilde{p}_{1},...,\tilde{p}_{p}$ can be used to represent a-priori beliefs about the importance of the predictors $x_{1},...,x_{p}$ and a good choice of them should lead to a fast convergence of the MH algorithm. Especially, in problems with a high number of predictors it is essential to assign $\tilde{p}_{1},...,\tilde{p}_{p}$ with informative values, since otherwise it might take a very large number of iterations until a MCMC method (see Section 2.3) used to simulate from the joint posterior might converge. According to $(21)$ in the case where $k_{t}>1$ a predictor $x_{\alpha}$ is added with probability $p_{\alpha}$ which might correspond to the a-priori importance of this predictor. Furthermore, a predictor $x_{\alpha}$ is removed with a probability proportional to $p_{\alpha}^{-1}$ which might be the inverse a-priori importance of $x_{\alpha}$. In the case where $k_{t}=1$ a predictor $x_{\alpha}$ is added with a probability proportional to $p_{\alpha}$ and it is impossible that the only predictor in the model is removed, this would violate the prior assumption that $k=|\mathcal{A}|>0$ with probability $1$. Using $(20)$ and $(21)$ the proposal distribution $q(\mathcal{A}_{t+1}|\mathcal{A}_{t})$ is finally defined by:

[TABLE]

Note that the probability mass function $q$ in equation $(22)$ is given by equation $(21)$, since $|k_{t+1}-k_{t}|=1$. Besides the easy evaluation of $q(\mathcal{A}_{t+1}|\mathcal{A}_{t})$ for a given input, one can easily simulate from this distribution. This is done by sampling from the distribution of $c_{h}$ followed by sampling from the conditional distribution of $\mathcal{A}_{t+1}|\mathcal{A}_{t},c_{h}$ where the condioning takes places with respect to the sample of $c_{h}$ generated in the previous step. The only thing remaining to define the overall proposal $q(\boldsymbol{\beta}_{t+1},\mathcal{A}_{t+1},g_{t+1},\sigma_{t+1}^{2}|\mathcal{A}_{t},g_{t},\sigma_{t}^{2},\mathbf{X})$ is the conditional proposal $q(\boldsymbol{\beta}_{t+1}|\mathcal{A}_{t+1},g_{t+1},\sigma_{t+1}^{2},\mathbf{X}_{\mathcal{A}_{t+1}})$. It is well known (see Fahrmeir et al. [2013]) that under the observation model $\mathbf{y}|\boldsymbol{\beta},\sigma^{2},\mathbf{X}\sim\mathcal{N}(\mathbf{X}\boldsymbol{\beta},\sigma^{2}\mathbf{\text{I}})$ with prior distributions $\boldsymbol{\beta}|\sigma^{2}~{}\sim\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{M})$ and $\sigma^{2}\sim\text{IG}(a,b)$ the conditional posterior of $\boldsymbol{\beta}$ given $\sigma^{2}$ is given by

[TABLE]

with

[TABLE]

[TABLE]

Thus, a natural choice for $q(\boldsymbol{\beta}_{t+1}|\mathcal{A}_{t+1},g_{t+1},\sigma_{t+1}^{2},\mathbf{X}_{\mathcal{A}_{t+1}})$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Finally, the overall proposal distribution can be written as

[TABLE]

One can easily simulate from the overall proposal by iteratively simulating from the factors in $(27)$ from right to left and conditioning on the values sampled up to the current step of execution.

3 Experimental studies

In this section the (predictive) performance of the proposed hierarchical Bayesian approach is analysed and compared with some other Bayesian and non-Bayesian methods. For this task we have implemented the MH algorithm described in Section 2.3 by using the interface between the programming languages R (R Core Team [2015]) and C. The comparison includes the lasso (Tibshirani [1996]), the adaptive lasso (alasso) (Zou [2006]), the elastic net (elastic) (Zou and Hastie [2005]), the Bayesian lasso (blasso) (Park and Casella [2008]), the Bayesian adaptive lasso (balasso) (Alhamzawi et al. [2012], Leng et al. [2014]) and the Bayesian elastic net (belastic) (Huang et al. [2015]). For our extensive comparisons three real-world datasets and also simulated datasets from two different artificial models are used. The performance comparison on the real-world datasets is carried out on the basis of a $5$-fold cross-validation. Performing a $5$-fold cross-validation on a given dataset means partitioning the dataset in $5$ equal sized subsets, iteratively selecting each of these subsets exactly once as testing data, whilst using the remaining subsets for training. Aggregating some measure of accuracy, respectively computed from each of the $5$ testing datasets, results in the overall performance of the models to evaluate. The performance comparison for the artificial models is accomplished by simulating multiple datasets from these models, then performing a train/test split for each simulated dataset, and finally again aggregating the single accuracies achieved. The accuracy measures chosen in this paper are the mean squared error (MSE) and the mean absolute deviation (MAD)

[TABLE]

where $n_{te}$ denotes the cardinality of the testing dataset and the $\widehat{y}_{i}$ denote the predicted target values, for $i=1,...,n_{te}$. For the aggregation of single accuracies we are using the median since it is more robust to outliers than the mean of the values. The median of MSEs is then denoted by MMSE and the median of MADs by MMAD. For our novel approach it should be stressed that, for obtaining a prediction $\widehat{y}^{*}$ corresponding to a test sample $\mathbf{x}^{*}$, an estimate of the expected value of the posterior predictive distribution $\mathbb{E}(y^{*}|\mathbf{x}^{*},\mathbf{y},\mathbf{X})$ is used. Using Monte Carlo integration this expected value can be estimated as follows:

[TABLE]

where $\boldsymbol{\beta}_{i},...,\boldsymbol{\beta}_{N}$ is a sample of size $N$ from the marginal posterior $p(\boldsymbol{\beta}|\mathbf{y},\mathbf{X})$. For all the non-Bayesian comparsion models the predicted target values $\widehat{y}_{1},...,\widehat{y}_{te}$ are computed by calculating the inner product of the test samples with the estimate $\widehat{\boldsymbol{\beta}}$ obtained by minimizing the corresponding model-specific objective function. For the Bayesian comparsion models the way the predictions are computed depend on the output provided by the R-packages used in this paper. In this work the R-function blasso included in the R-package monomvn (Gramacy [2018]) is applied in order to perform Bayesian lasso regression. Since the function blasso returns a sample from the marginal posterior $p(\boldsymbol{\beta}|\mathbf{y},\mathbf{X})$, as for the model proposed by us, the predictions are computed by estimating the posterior predictive mean. Moreover, we use the R-Function brq included in the package Brq (Alhamzawi [2018]) for applying the Bayesian adaptive lasso regression. It should be mentioned that this R-function is not developed for the classical linear regression setting where the conditional mean of the response variable is modelled, rather it is developed for the more general setting where any conditional quantile is modelled. For purposes of the comparision study in this paper the conditional median, i.e. the conditional $50\%$-quantile, is selected. The R-function produces a Bayes estimate of $\boldsymbol{\beta}$ which is then used for prediction (i.e. inner product computation with the test samples). In order to use the Bayesian elastic net the function EBglmnet included in the package EBglmnet (Huang and Liu [2016]) is used in this work. This function returns the posterior mean of the nonzero regression coefficients, which is then used for prediction (again via inner product computation with the test samples).

It is also important to mention how (inside the $5$-fold cross-validation for real-world datasets, or for the multiple simulated datasets) the hyperparameters of the trained models are determined. For the frequentist models, which are all implemented in the R-package glmnet (Friedman et al. [2010]), the R-function cv.glmnet is used, which is also included in this package. The penalization factor $\lambda$ of the lasso regression is chosen based on a $10$-fold cross-validation on a one dimensional grid with the MAD as accuracy measure. Analogously, the penalization factors $\lambda_{1},\lambda_{2}$ of the elastic net are selected based on a $10$-fold cross-validation on a now two dimensional grid, again with MAD as accuracy measure. For the adaptive lasso the chosen penalization factors are the absolute values of the inverse regression coefficients obtained from the respective ridge-penalized problem, the penalization factor of which is also determined via $10$-fold cross-validation. For the Bayesian version of the lasso and the adaptive lasso the default values of the corresponding functions blasso and brq are assigned to the hyperparameters of the a-priori distributions. For further details please refer to the help pages of these packages. However, it should be mentioned that for the Bayesian lasso the option of additional model selection based on a Reversible Jump MCMC algorithm is set to true for the experimental studies in this paper. Moreover, for the Bayesian elastic net the hyperparameters in the prior distribution are selected via $5$-fold cross-validation. This is done by using the function cv.EBglmnet, which is also available in the package EBglmnet. The hyperparameters of our proposed Bayesian model are chosen differently in the following experimental studies to indicate useful possible choices. These specific choices are described later on in this article.

In order to obtain i.i.d. samples from the model-specific posterior distributions the R-functions used in this paper merely save a fraction of all samples generated by the MCMC methods implemented. For the Bayesian elastic net, the decision which samples to store is already taken by the implementation itself and cannot be user-defined. All other Bayesian models go along with implementations which allow for a user specific determination of this fraction. To do so the following procedure is chosen: In a first step, respectively, the first $10,000$ samples are deleted and, in a second step, all samples except of every $10$-th one are deleted. Note that samples generated in the first step before convergence of the used MCMC algorithms are deleted, while thinning done in the second step should guarantee that the remaining samples can be considered as being independent. For each of the observed Bayesian models $50,000$ (dependent) samples are generated, except for the Bayesian adaptive lasso where $70,000$ ones are produced. This results in $4,000$ i.i.d. samples each, except of $6,000$ samples for the adaptive lasso. The reason for the special treatment is that in our experiments the Bayesian adaptive lasso turned out to require more samples to get stable estimates of the target variables.

3.1 Real-world studies

In this section the performance of the proposed method is compared with the performance of other well-established methods based on three real-world datasets. The comparison is made along the lines of the above considerations (see beginning of Section 3).

3.1.1 The diabetes dataset

In this section the popular diabetes dataset, originally used by Efron et al. [2004] to examine the so-called LARS algorithm, is used for the performance evaluation. This dataset includes ten predictors, age, sex, body mass index, average blood pressure, and six blood serum measurements, measured from $n=442$ diabetes patients. The target variable of interest is a quantitative measure of disease progression one year after baseline. In the R-package care (Zuber and Strimmer. [2017]) a standardized version (all variables have zero mean and variance equal to one) of the diabetes dataset is available. This version is used in our experiments.

The hyperparameters of the proposed Bayesian approach are specified as follows: In an empirical Bayes manner, the parameters $p_{1},...,p_{p}$ which belong to the prior of $\mathcal{A}$ and represent the a-priori belief in the importance of the predictors $x_{1},...,x_{p}$ are assigned appropriate multiples (they must sum up to one) of the absolute values of the regression coefficients of the ridge regularized model which also determines the penalization factors of the adaptive lasso model. Moreover, the parameters $\tilde{p}_{1},...,\tilde{p}_{p}$ (see Section 2.3) are specified identically to the parameters $p_{1},...,p_{p}$. As already mentioned in Section 2.1, one could define the mapping $\tilde{p}$ as the probability mass function of a zero truncated binomial distribution with size parameter equal to the number of predictors $p$ and the second parameter chosen in such a way that the mean of the distribution equals the a-priori most probable model size. Thus, a zero truncated binomial prior is assigned to the unknown model size $k$. However, in general it is somewhat cumbersome to determine the second parameter in such a way that the distribution has a given mean. Note that the expected value of a zero truncated binomial distribution with size parameter $p$ and second parameter $q$ equals $pq/(1-(1-q)^{p})$. Further, note that by the Abel-Ruffini theorem there is no algebraic solution to the general polynomial equations of degree five or higher with arbitrary coefficients. A classical binomial distribution with size parameter $p$ and second parameter equal to $\mu/p$ has mean $\mu$. Thus, a zero truncated binomial distribution with these parameters will result in a distribution with expected value not too far away from $\mu$. Therefore, defining $\tilde{p}$ as probability mass function of such a zero truncated distribution is a good alternative and our final choice in the performance evaluation of the diabetes dataset. Now the only question remaining is how to specify $\mu$. Since there is no a-priori knowledge available one could determine a useful value of $\mu$ via cross-validation in an empirical Bayes manner. However, to save time we empirically tried some possible values for $\mu$ and very soon ended up with the specification $\mu=9.9$ (which results in a good performance as one can see later on in this section). The parameters $a$ and $b$ are both set to $0.001$ which will also hold for the remaining experimental studies, as already reported in Section 2. Finally, the tuning parameters (see Section 2.3) $p_{h},\varepsilon_{\sigma}$ and $\varepsilon_{g}$ are respectively set to the values $0.4,0.1$ and $60$.

Performing a $5$-fold cross-validation as described at the beginning of Section 3 results in Table 1. The proposed approach achieves the lowest MMSE as well as the lowest MMAD and thus performs better than all methods under comparison.

3.1.2 The prostate cancer dataset

In this section the performance of the seven considered methods is compared on the prostate cancer data provided by Stamey et al. [1989]. This dataset examines the relationship between the level of prostate specific antigen (lpsa) and eight clinical measures from $99$ patients who were about to receive a radical prostatectomy. The eight measurments are log cancer volume (lcavol), log prostate weight (lweight), age, log benign prostatic hyperplasia amount (lbph), seminal vesicle invasion (svi), log capsular penetration (lcp), Gleason score (gleason), and percentage Gleason scores 4 or 5 (pgg45). This dataset is included e.g. in the R-package lasso2 (Lokhorst et al. [2014]). The predictors are standardized, consequently we only have to zero-center the target variable lpsa to make the data compatible with the model proposed by us.

The hyperparameters of our new Bayesian approach are specified as follows: The parameters $p_{1},...,p_{p}$ and $\tilde{p}_{1},...,\tilde{p}_{p}$ are specified the same way as in Section 3.1.1. Again, $\tilde{p}$ is defined as probability mass function of a zero truncated binomial distribution with size parameter $p$ and second parameter $\mu/p$, but now $\mu$ is specified as the number of nonzero regression coefficients from the corresponding lasso model. The remaining parameters are also defined as in Section 3.1.1.

Thus, performing a $5$-fold cross-validation as described at the beginning of Section 3 results in Table 2.

Again, our method achieves the lowest MMSE, while some other methods go along with better MMADs on the prostate cancer dataset.

3.1.3 The ozone dataset

In this section the ozone datset originally studied by Breiman and Friedman [1985] is used for the performance comparison. The data was measured in the area of Upland, CA, east of Los Angeles on 330 days in 1976. The target variable is the daily maximum one-hour-average ozone reading (parts per million) and the predictors are:

Temperature (degrees F)

• 2.

Inversion base height (feet)

• 3.

Daggett pressure gradient (mmHg)

• 4.

Visibility (miles)

• 5.

Vandenburg 500 millibar height (m)

• 6.

Humidity (%)

• 7.

Inversion base temperature (degrees F)

• 8.

Wind speed (mph)

The dataset is available in the R-package gclus (Hurley [2012]). We decide to examine the regression model including all linear, quadratic and two-way interactions, resulting in $44$ possible predictors. Moreover, we zero center and standardize (variance equal to one) all variables previous to the training of the diverse models. We want to mention that in contrast to a statement at the beginning of Section 3 in this performance evaluation for all Bayesian models $100,000$ MCMC samples are drawn and not $50,000$ or $70,000$ as stated there. This is the only experimental study where this statement is violated.

The hyperparameters of the newly proposed Bayesian approach are specified as follows: We assign to each of the parameters $p_{1},...,p_{p}$ and $\tilde{p}_{1},...,\tilde{p}_{p}$ the value $1/p$, i.e. assume equal a-priori importance of the predictors. The function $\tilde{p}$ is defined in such a way, that it maps a value $k\in\left\{1,...,p\right\}$ to a value proportional to the third power of the probability $P(x=k)$, where $x$ denotes a random variable which is zero truncated binomial distributed with size parameter $p$ and second parameter $7/p$. We empirically discover that taking the power of three can be a helpful step to guarantee fast convergence of the proposed MH algorithm (see Section 2.3) in case of an increasing number of predictors and a prior based on a binomial distribution. By choosing the constant of proportionality in the definition of $\tilde{p}$ in such a way, that $\tilde{p}$ defines a probability distribution, the resulting distribution will have a smaller variance than the generating binomial distribution. Finally, the tuning parameters (see Section 2.3) $p_{h},\varepsilon_{\sigma}$ and $\varepsilon_{g}$ are set to the values $0.5,0.1$ and $60$, respectively.

Performing a $5$-fold cross-validation as described at the beginning of Section 3 results in Table 3. The new approach achieves the lowest MMAD and a MSE comparable to the other methods .

3.2 Simulation studies

In this section, on the basis of simulated data corresponding to two different artificial models, the performance of our new method is compared with the performance of other well-established methods. The artificial models on purpose include many correlated predictors, from which only a small subset is predictive, i.e. has regression coefficients different from zero. This should reflect difficult variable selection problems that more and more organizations have top face in their daily life nowadays and in the future. From both artificial models multiple datasets are simulated. In particular, $100$ datasets are sampled, each according to three different settings. First the sampling takes place with $n=50$ training observations, then with $n=100$ training observations and finally with $n=200$ training observations. The number of test observations is always the same and given by the value $n_{te}=200$. Moreover, all the observations are drawn independently from each other. For each sampled dataset the entire target vector (includes target values from both the training set and the testing set) is standardized via division by the sample standard deviation of the training target values. This simplifies finding a good choice of the tuning parameters for the MH algorithm proposed. The scaling allows for choosing them equal to specifications used in Section 3.1. The predictors are not scaled at all since they are sampled from distributions with mean zero and standard deviation one.

3.2.1 Simulation study 1

In this simulation study data are generated from the model

[TABLE]

with $(\beta_{2},\beta_{11},\beta_{21},\beta_{51},\beta_{71},\beta_{81})=(-2.5,-2,-1.5,1.5,2,2.5)/\sqrt{3}$ and the remaining coefficients equal to zero. This model is inspired by those used to evaluate the performance of the spike-and-slab lasso by Ročková and George [2018].

The hyperparameters of the proposed Bayesian approach are specified as follows: The parameters $p_{1},...,p_{p}$ are defined as multiples of the absolute sample correlations between the training responses $\mathbf{y}$ and the training predictors $\mathbf{x}_{1},...,\mathbf{x}_{p}$. Note, that the multiplicative factor has to be defined in such a way that the parameters sum up to one. Moreover, the parameters $\tilde{p}_{1},...,\tilde{p}_{p}$ are specified identical to the parameters $p_{1},...,p_{p}$. The function $\tilde{p}$ maps each value $k\in\left\{1,...,p\right\}$ to a value proportional to the third power of the probability $P(x=k)$, where $x$ denotes a random variable which is zero truncated binomially distributed with size parameter $p$ and second parameter $2.5/p$. Finally, the tuning parameters $p_{h},\varepsilon_{\sigma}$ and $\varepsilon_{g}$ are set to the values $0.5,0.1$ and $60$, respectively.

Aggregating the accuracy measures obtained by training the models to be compared, as described at the beginning of Section 3, results in Table 4. The proposed approach achieves the lowest MMSE as well as the lowest MMAD and thus performs better than all other methods.

In the Figures 4-4 boxplots of the MSEs and MADs obtained from the $100$ simulated datasets with $n=50$ and $n=100$ training observations can be found. These plots again illustrate that our new model performs significantly better than the other models considered. Moreover, one can see that the implementation of the Bayesian adaptive lasso has some serious problems with these simulated data, due to numerical instabilities.

To provide some more insights into the new method proposed we consider a representative dataset simulated from model $(61-65)$ with $n=100$ observations. In Figure 6 one can see the Markov chain of the model size $k$ (i.e. the number of non-zero regression coefficients) produced from our MH algorithm. Further, in Figure 6 the relative frequencies of a given regression coefficient being non-zero are plotted for the reduced Markov chain. By the latter one we mean the Markov chain remaining after the deletion step performed to obtain i.i.d. samples. The relative frequencies can be interpreted as the a posteriori probabilities of the predictors being included in the model. All the truly non-zero regression coefficients are colored blue, while the remaining coefficients are colored red. One can see that the proposed approach selects the true data generating model pretty well. Finally, we want to report that the acceptance rate of the proposed MH algorithm is given by the value $0.20672$, which is a good acceptance rate indicating that the algorithm converges fast.

3.2.2 Simulation study 2

In this simulation study data are generated from the model

[TABLE]

with $(\beta_{2},\beta_{11},\beta_{21},\beta_{51},\beta_{71},\beta_{81})=(-2.5,-2,-1.5,1.5,2,2.5)/\sqrt{3}$ and the remaining coefficients equal to zero. Simulating from this model results in more difficult datasets in context of variable selection than simulating from the model defined in Section 3.2.1. There are $400$ additional potential predictors and the pairwise correlations between them are increased from $0.6$ to $0.75$.

The hyperparameters of the newly proposed Bayesian approach are specified, in an empirical Bayes manner, as follows: The parameters $p_{1},...,p_{p}$ are assigned appropriate multiples (they must sum up to one) of the absolute values of the regression coefficients of the ridge regularized model, which also determines the penalization factors of the adaptive lasso model. The parameters $\tilde{p}_{1},...,\tilde{p}_{p}$ are chosen analogously, except that the third power of the ridge coefficients is chosen. The function $\tilde{p}$ maps each value $k\in\left\{1,...,p\right\}$ to a value proportional to the third power of the probability $P(x=k)$, where $x$ denotes a random variable which is zero truncated binomially distributed with size parameter $p$ and second parameter $1.75/p$. Finally, the tuning parameters $p_{h},\varepsilon_{\sigma}$ and $\varepsilon_{g}$ are set to the values $0.5,0.1$ and $60$, respectively.

Due to numerical issues in the corresponding R-implementation the Bayesian adaptive lasso has to be excluded from the model comparison. Aggregating the accuracy measures obtained by training the remaining models results in Table 5. The proposed approach achieves the lowest MMSE as well as the lowest MMAD in the settings with $n=100$ and $n=200$ training observations. In the setting with $n=50$ training observations it achieves a comparable accuracy. The Figures 10-10 show boxplots of the MSEs and MADs obtained from the $100$ simulated datasets with $n=50$ and $n=100$ training observations. In the setting $n=100$ our method outperforms the other ones. The median of the MADs as well as the median of the MSEs lies below all corresponding boxes of the other methods.

4 Conclusion

In this article we have presented a novel Bayesian approach to cope with the problem of variable selection in the multiple linear regression model with dependent predictors. To make our method robust to challenges such as multicollinearity, or the number of observations being smaller than the number of potential predictors, a $g$-prior with an additional ridge parameter was assigned to the unknown regression coefficients. While other authors apply Gibbs sampling to simulate from the joint posterior distribution we presented an intelligent Metropolis Hastings algorithm for this task. Thus, we do not need to compute the conditional distributions needed for the Gibbs sampling algorithm, which often are not available in an analytically closed form, and have to be approximated. Further, we have shown that the proposed Bayesian approach is consistent in terms of model selection under some nominal assumptions. Experimental studies with three different real-world datasets as well as simulated datasets from two artificial models demonstrated the good performance of the presented method. In particular, on the simulated datasets our new approach performed significantly better than all other well-established methods.

Acknowledgments

The publication is one of the results of the project iDev40 (www.idev40.eu). The iDev40 project has received funding from the ECSEL Joint Undertaking (JU) under grant agreement No 783163. The JU receives support from the European Union’s Horizon 2020 research and innovation programme. It is co-funded by the consortium members, grants from Austria, Germany, Belgium, Italy, Spain and Romania.

Appendix A Proof of Theorem 1

Note that for showing $(14)$ it suffices to show that

[TABLE]

holds, where $M_{\mathcal{A}}$ denotes the true model and $M_{\mathcal{A}^{\prime}}$ denotes another model different from the true one. The sets $\mathcal{A}$ and $\mathcal{A}^{\prime}$, respectively, denote the index sets of the active predictors (predictors probably being different from zero) of the two models $M_{\mathcal{A}},M_{\mathcal{A}^{\prime}}$ and thus uniquely determine them. In a first step of the proof we compute the joint posterior of $\mathcal{A}$ and $g$ up to a constant of proportionality. Simple linear algebra and the facts that the densities of multivariate normal distributions and inverse gamma distributions integrate to one reveal for $n>\zeta$:

[TABLE]

where

[TABLE]

[TABLE]

and $k$ denotes the cardinality of $\mathcal{A}$. In principle, the required marginal posterior of $\mathcal{A}$ can be found by integrating out $g$ in equation (A.2). However, there exists no closed-form solution of the corresponding integral. To overcome this problem, a lower and an upper bound of the marginal posterior are derived. The numerator of the quotient in (A.1) is then replaced by the upper bound and the denominator by the lower bound such that the resulting quotient is an upper bound of the original one. Showing that the upper bound of the original quotient converges to zero will immediately imply that (A.1) holds. To determine these bounds, at first, a property of the Gaussian hypergeometric function given by

[TABLE]

and convergent for $|z|<1$ with $\tau>\beta>0$ and for $|z|=1$ only if $\tau>\alpha+\beta$ and $\beta>0$ is derived. It can be easily validated that for $1>z\in\mathbb{R}$ and $\tau>\beta>0$ the following equation holds:

[TABLE]

Note that the second equality holds by identity $15.3.6$ in Abramowitz and Stegun [1964]. Above equation directly translates to the following identity which will later on be used to evaluate an upper and a lower bound of $p(\mathcal{A}|\mathbf{y},\mathbf{X})$:

[TABLE]

Determination of the upper bound:

At first we observe that:

[TABLE]

This inequality holds since it can be reduced to the well known inequality $1+x\leq\operatorname{exp}(x)~{}\forall x\in\mathbb{R}$. Moreover, we observe that:

[TABLE]

Using the inequalities (A.7) and (A.8) an upper bound of $p(\mathcal{A}|\mathbf{y},\mathbf{X})$ can be derived:

[TABLE]

Note that (A.12) is a direct application of (A.6) to (A.11). Equality (A.13) holds by identity $15.3.3$ in Abramowitz and Stegun [1964]. The last equality holds, since $2F_{1}(0,\beta,\tau,z)$ equals the area under the density of a beta distribution with parameters $\beta$ and $\tau-\beta$ and thus is given by the value $1$. Applying Stirling’s formula

[TABLE]

to (A.14), which was also used in the model consistency proof of Wang et al. [2015], allows another approximate simplification of the right-hand side expression of (A.14) for large $n$ through:

[TABLE]

Note that the approximation sign in Strirling’s formula above means that the ratio of the two sides converges to one as $x$ goes to infinity.

Determination of the lower bound:

Observing that $(g+1)^{r}\geq g^{r}\geq 0$ for $g\geq 0,r\geq 0$ and further observing that $\operatorname{exp}[g^{-1}(1-R_{\mathcal{A}}^{2})^{-1}]\geq 1+g^{-1}(1-R_{\mathcal{A}}^{2})^{-1}~{}\forall g\in\mathbb{R}$ since $\operatorname{exp}(x)\geq 1+x~{}\forall x\in\mathbb{R}$, a lower bound for $p(\mathcal{A}|\mathbf{y},\mathbf{X})$ is given by:

[TABLE]

Note that (A.21) holds since the function to integrate in (A.20) is proportional to the density of an inverse gamma distribution. Moreover, $c_{1}$ and $c_{2}$ denote constants which are independent of $n$.

Using the bounds (A.15) and (A.22), an upper bound for the quotient in (A.1) can be obtained for large $n$:

[TABLE]

Note that $c_{3}$ denotes a constant independent of $n$.

To finally prove model selection consistency, two cases are considered:

1. (i)

If $M_{\mathcal{A}}\nsubseteq M_{\mathcal{A}^{\prime}}$ the probability limit of the upper bound (A.25) for $n\rightarrow\infty$ evaluates as follows:

[TABLE]

Note that equation (A.27) holds by the equations $(15)$ and $(17)$. Further, note that the term $(\sigma^{2}/(\sigma^{2}+b_{\mathcal{A}^{\prime}}))^{n/2}$ is element of the interval $(0,1)$ and thus converges exponentially fast to zero. Therefore, the other factor of the product corresponding to equation (A.27) does not influence the limit at all, since it is polynomial in $n$. 2. (ii)

If $M_{\mathcal{A}}\subseteq M_{\mathcal{A}^{\prime}}$ according to the work of Fernández et al. [2001] the following holds:

[TABLE]

where $s$ has a \mbox{\large\chi}^{2} distribution with $k^{\prime}-k$ degrees of freedom and $\overset{D}{\longrightarrow}$ denotes convergence in distribution. Thus, for $n\rightarrow\infty$ the limit of the upper bound (A.25) evaluates as follows:

[TABLE]

Note that $k^{\prime}>k$ since $M_{\mathcal{A}}\neq M_{\mathcal{A}^{\prime}}$.

Finally, in both cases the upper bound of the quotient in (A.1) converges to [math] and thus also the quotient itself since it is greater or equal to [math] and the proof is completed.

$\square$