# Variational Inference for Sparse Poisson Regression

Mitra Kharabati<EMAIL_ADDRESS>Dpt. of Statistics, School of
Mathematics, Statistics and Computer Science, University of Tehran, Tehran,
Iran  Morteza Amini Corresponding author, e-mail<EMAIL_ADDRESS>Dpt.
of Statistics, School of Mathematics, Statistics and Computer Science,
University of Tehran, Tehran, Iran

###### Abstract

We have utilized the non-conjugate VB method for the problem of the sparse
Poisson regression model. To provide an approximated conjugacy in the model,
the likelihood is approximated by a quadratic function, which provides the
conjugacy of the approximation component with the Gaussian prior to the
regression coefficient. Three sparsity-enforcing priors are used for this
problem. The proposed models are compared with each other and two frequentist
sparse Poisson methods (LASSO and SCAD) to evaluate the prediction
performance, as well as, the sparsing performance of the proposed methods.
Throughout a simulated data example, the accuracy of the VB methods is
computed compared to the corresponding benchmark MCMC methods. It can be
observed that the proposed VB methods have provided a good approximation to
the posterior distribution of the parameters, while the VB methods are much
faster than the MCMC ones. Using several benchmark count response data sets,
the prediction performance of the proposed methods is evaluated in real-world
applications.

##  1 Introduction

Poisson regression (see for example Hilbe, 2014; Hayat and Higginsm, 2014, and
references therein) and its special cases, have been applied to many problems
in a wide area of science (see Legisso et al., 2023; Fantin et al., 2023;
Solomon et al., 2022; Winter and Bürkner, 2021; Chau et al., 2018; Coxe et
al., 1992; Alonso et al., 2015, among others). The Bayesian estimation of the
parameters of the Poisson regression model is considered by many researchers.
The MCMC sampling methods, including the Metropolis-Hastings algorithm
(Hastings, 1970) and Gibbs sampler (Geman and Geman, 1984) are used as common
methods for estimation of the posterior distribution in these models. El-
Sayyad (1973) developed the Bayesian analysis of the Poisson regression model.
Christiansen and Morris (1997) proposed a Poisson model with a non-informative
prior distribution for the dummy values and non-reversible gamma distributions
for the individual parameters. Hierarchical Bayesian Poisson regression models
are considered by Kim et al. (2013). There are many other works including
D’Angelo et al. (2022); Tomal et al. (2022); Zhang and Ye (2014); Koç et al.
(2013); Khafrim et al. (2008); Viallefont et al. (2002); Mohammadi et al.
(2016); Mzoughia et al. (2018); Dagne (2010); Sellers and Premeaux (2021);
Sellers et al. (2012); Xie et al. (2014) and Consul and Famoye (1992) where
the Bayesian inference is applied to Poisson regression model and its
varieties, using MCMC sampling methods.

A challenging problem in all regression models is variable selection, that is
to choose a subset of covariates that are the most effective on the response
variables. The sparse estimation is one of the most common solutions to this
problem. In the frequentist methods, the sparsity of the coefficients of the
regression model is obtained using sparsity penalty functions. The least
absolute shrinkage and selection operator (LASSO- Tibshirani, 1996), elastic
net (Li and Lin, 2010), smoothly clipped absolute deviation (SCAD- Fan and Li,
2001), mini-max concave penalty (MCP- Zhang, 2010) are the most well-known
penalized regression models for sparse estimation of the regression
parameters. The sparse Poisson regression model is considered by many
researchers, including Ahmed (2014); Algamal (2015); Algamal and Lee (2015);
Blazére et al. (2014); Famoye and Rothe (2003); Guastavino and Benvenuto
(2019); Igarashi et al. (2018); Hossain and Ahmed (2012); Ivanoff et al.
(2016); Jia et al. (2013, 2019); Noori Asl et al. (2022); Peng (2023); Saishu
et al. (2021); Yousif and Gatea (2021); Zou and Li (2008); Li and Cevher.
(2015) among others.

From a Bayesian point of view, the sparsity of the regression parameters is
achieved by using sparsity-enforcing priors. Some examples of such priors are
Laplace (Seeger, 2008), horseshoe (Chan and Vasconcelos, 2009), Student’s t
(Tipping, 2001), and spike-and-slab (Mitchell and Beauchamp, 1988; Geweke, J.,
1996; George and McCulloch, 1997) priors. A huge number of researches are
devoted to sparse Bayesian regression models. Some examples are Picard and
Cook (1984); Shao (1993); Xu and Ghosh (2015); Richardson et al. (2010); Shin
and Liu (2022); O’hara and Sillanpää (2009); Griffin and Brown (2013); Bondell
and Reich (2012); Ročková and George (2014); Ji and Schmidler (2013);
Nishimura and Suchard (2022); Polson et al. (2014); Seeger (2008); Carvalho et
al. (2009); Mitchell and Beauchamp (1988); Geweke, J. (1996); George and
McCulloch (1997).

The Bayesian sparse regression has also been of interest to many researchers.
Almost all of the papers in this field use the MCMC sampling for estimation of
the parameters of the sparse Poisson regression model (Gajda et al., 2010;
Dvorzak and Wagner, 2016; Powers et al., 2010; Benson and Friel, 2021; Tang
and Martin, 2023; Jiang, 2007; Shankar Raman, 2012; Crowther et al., 2012; Bar
et al., 2020; Raman and Roth, 2009; Bai, 2020; Jeong and Ghosal, 2021). Some
of these works consider the generalized linear models (GLM), which include
Poisson, logistic, and many other non-Gaussian regression models.

While the MCMC method is the most popular method for approximation of the
posterior distribution in Bayesian inference, it suffers a high computational
cost especially when the dimension of the regression parameters is high. There
are a variety of methods for approximation of the posterior distribution with
a lower computational cost, including integrated nested Laplace approximation
(INLA- Rue et al., 2010), mean field variational Bayes (VB- Blei et al.,
2017), expectation propagation (EP- Minka, 2001, 2013) and message passing
(MP- Winn et al., 2005). INLA and VB are the most frequently used methods in
this category. The VB method is applied in many problems, especially the
Bayesian regression models (see Šmídl and Quinn, 2006; Minka, 2000; Wand et
al., 2011; Broderick et al., 2013; Drugowitsch, 2013; Goldsmith et al., 2011;
Pham et al., 2013; Titsias and Lázaro-Gredilla, 2011, and the references
therein). The VB method is also applied to the sparse linear regression models
(Ormerod et al., 2017; Ray and Szabó, 2022; Titsias and Lázaro-Gredilla, 2011;
Huang et al., 2016).

The VB method is also applied to the Poisson regression model by a few
authors. Luts and Wand (2015) developed VB for Bayesian semi-parametric
regression with count data, using Poisson and Negative Binomial regression
models. Chen et al. (2010) applied the VB for the point process generalized
linear models in the neural spike trains analysis.

The method used in Luts and Wand (2015) was based on the non-conjugate VB
technique. In this method an assumed density filtering technique is used along
with the mean filed VB for approximation of the posterior density. Precisely
speaking, since the approximated component for the coefficients of the
regression model is not conjugate with its Gaussian prior, we assume a
Gaussian component and optimize the mean and the covariance function by
minimizing the Kullback-Leibler divergence between the unnormalized posterior
and its approximation. Another method in non-conjugate VB is the approximation
method, proposed by Jaakkola and Jordan (2000), which was successfully applied
to the logistic regression problem. In this method, the likelihood function is
approximated by a quadratic function, which provides the conjugacy of the
approximation component with the Gaussian prior of the regression coefficient.
The method of Jaakkola and Jordan (2000) has been applied to the sparse
logistic regression model by Zhang et al. (2019).

In this paper, we have used the VB method for the problem of the sparse
Poisson regression model. A similar approximation method to that of Jaakkola
and Jordan (2000) is used to provide conjugacy in the model approximation.
Three sparsity-enforcing priors are used for this problem and the VB
approximated formulas are calculated for all three models. The proposed models
are compared with each other and two frequentist sparse Poisson methods (LASSO
and SCAD) to evaluate the prediction performance, as well as, the sparsing
performance of the proposed methods. Throughout a simulated data example, the
accuracy of the VB methods is computed compared to the corresponding benchmark
MCMC methods. It can be observed that the proposed VB methods have provided a
good approximation to the posterior distribution of the parameters, while the
VB methods are much faster than the MCMC ones. Using several benchmark count
response data sets, the prediction performance of the proposed methods is
evaluated in real-world applications.

The rest of the paper is organized as follows. In Section 2, the mean filed
variational Bayes method is introduced. Three variational Bayesian sparse
Poisson regression models are then proposed in Section 3. The derivation of
the components as well as the evidence lower bound (ELBO) for these methods
are given in Appendix A. In Section 4, the posterior predictive density
function is approximated for the prediction of the response variable using new
samples of the covariates. Section 5 presents a simulation study for
evaluating the prediction and sparsing performance of the three proposed
methods and comparing them with two frequentist sparse Poisson regression
methods (LASSO and SCAD Poisson regression models). A single replication of
the simulation study is also used for computation of the accuracy of each
approximation component compared to its corresponding benchmark MCMC method
and also a visual comparison of the MCMC empirical posteriors with their
approximated counterparts. Some of these plots are given in Appendix B for the
case of briefness. Finally, in Section 6, the test set prediction performance
of the proposed methods is compared with those of the frequentist
alternatives, using 6 benchmark count response data sets. Some concluding
remarks are given in Section 7.

##  2 Mean field variational Bayes

Let $x$ be a vector of observed data, and $\theta$ be a parameter vector with
joint distribution $P(x,\theta)$. In the Bayesian inference framework, the
inference about $\theta$ is done based on the posterior distribution
$p(\theta|x)=p(x,\theta)/p(x)$ where $p(x)=\int p(x,\theta)d\theta$.

Variational Bayes (VB) is a method for finding an approximate distribution
$q(\theta)$ of the posterior distribution $p(\theta|x)$, by minimizing the
Kullback-Leibler divergence ${\rm KL}\left[q(\theta)||p(\theta|x)\right]$ as a
measure of closeness (Tran et al., 2023). In the mean-field VB we assume that
parameter vector $\theta$ is divided into $M$ partitions
$\\{\theta_{1},\ldots,\theta_{M}\\}$ and we want to approximate $p(\theta|x)$
by

$q(\theta)=\prod_{j=1}^{M}q(\theta_{j}).$

.

The best VB approximation $q^{\ast}(\cdot)$ is then obtained as

$q^{\ast}=\arg\min_{q}{\rm KL}\left[q(\theta)||p(\theta|x)\right],$

where

${\rm KL}\left[q(\theta)||p(\theta|x)\right]=\int
q(\theta)\log\frac{q(\theta)}{p(\theta|x)}d\theta=\log p(x)-\int
q(\theta)\log\frac{p(x,\theta)}{q(\theta)}d\theta,$ (1)

which results in

$\log
q(\theta_{j})=E_{-\theta_{j}}\left[\log{P}(X,\theta)\right],j=1,\ldots,M.$ (2)

It is clear from (1) that minimizing KL is equivalent to maximizing the
evidence lower bound (ELBO), defined as

${\rm ELBO}(q)=\int q(\theta)\log\frac{p(x,\theta)}{q(\theta)}d\theta$

ELBO is equal to the model log-likelihood when the KL divergence is zero,
which means a perfect fit. When the fit is not perfect ${\rm
ELBO}[q(\theta)]<\log p(X)$.

##  3 Variational Bayesian sparse Poisson regression models

In the following subsections, we consider three different sparsity enforcing
priors for the coefficients of the Poisson regression model. Then, we apply an
approximation method similar to that of Jaakkola and Jordan (2000) to the
likelihood function, to provide approximated conjugacy for the regression
coefficient vector. Finally, we propose the VB components for the parameters
of each model. The computation process of these components is given in
Appendix A.

###  3.1 Laplace prior

Suppose that $y_{i}\sim\text{Poiss}(\lambda_{i}=\exp(X_{i}\beta)),\quad
i=1,\ldots,n$, independently, where
$\beta=(\beta_{0},\beta_{1},\ldots,\beta_{p})^{\top}$. A sparsity enforcing
prior for $\beta$ is the Laplace prior (Seeger, 2008). Thus, we consider the
independent priors $\beta_{j}\sim\;\text{Laplace}(0,\lambda),\;j=1,\ldots,p$.
It is well-known that this prior is identical with the hierarchical model

$\beta_{j}\sim\;\mathcal{N}(0,\tau_{j}),\quad\tau_{j}\sim\;\text{exp}\left(\frac{\lambda^{2}}{2}\right),\quad
j=1,\ldots,p.$

Adding the intercept model, the Bayesian sparse Poisson regression model
considered here is as follows

$\displaystyle y_{i}\sim$
$\displaystyle\;\text{Poiss}(\lambda_{i}=\exp(\beta_{0}+X_{i}\beta)),\quad
i=1,\ldots,n$ $\displaystyle\beta_{j}\sim$
$\displaystyle\;\mathcal{N}(0,\tau_{j}),\quad\tau_{j}\sim\;\text{exp}\left(\frac{\lambda^{2}}{2}\right),\quad
j=1,\ldots,p,$ $\displaystyle\beta_{0}\sim$
$\displaystyle\;\mathcal{N}(0,\sigma^{2}),\;\;\;\;\lambda\sim\text{Gamma}(\nu,\delta)$
$\displaystyle\;\sigma^{2}\sim$ $\displaystyle\text{Inv-
Gamma}\left(\frac{1}{2},\frac{1}{a}\right),\;\;a\sim\text{Inv-
Gamma}\left(\frac{1}{2},\frac{1}{A}\right)$ (3)

The directed acyclic graph (DAG) for model (3) is given in Figure 1, which
shows the relation between data elements (squares), model parameters
(circles), and fixed hyper-parameters (polygons).

Figure 1: The directed acyclic graph (DAG) for model (3). Data elements are
shown by squares, model parameters are shown by circles and fixed hyper-
parameters are shown by polygons.

Because of the $\exp(\cdot)$ in the likelihood function of model (3), the VB
component for $\beta$ is not conjugate. Similar to the idea proposed by
Jaakkola and Jordan (2000) (for the logistic regression model), we consider
the following inequality, for $\xi\leq x,$

$e^{x}\geq
e^{\xi}\left[(1-\xi)(1+x)+\frac{x^{2}}{2}+\frac{\xi^{2}}{2}\right].$ (4)

The inequality (4) is exact when $\xi=x$. Thus, we let
$\xi_{i}=E_{q}(X_{i}\beta)=X_{i}\mu_{\beta}$, for $i=1,\ldots,n$, where
$\mu_{\beta}$ is given in (6).

The $q$ density is assumed to be

$q(\beta_{0},\beta,\tau,\lambda,\sigma^{2},a)=q(\beta_{0})q(\beta)q(\tau)q(\lambda)q(\sigma^{2})q(a),$

where $\tau=(\tau_{1},\ldots,\tau_{p})^{\top}$.

Define $M_{\xi}=e^{\xi}(1-\xi)$ and
$S_{X}^{\xi}=\sum_{i=1}^{n}e^{\xi_{i}}X_{i}^{T}X_{i}$ Using inequality (4), we
have

$\displaystyle\log P(y|\beta)=$
$\displaystyle-\sum_{i=1}^{n}\exp(X_{i}\beta)+\sum_{i=1}^{n}y_{i}X_{i}\beta-\sum_{i=1}^{n}\log
y_{i}!$
$\displaystyle\leq-\sum_{i=1}^{n}e^{\xi_{i}}\left[(1-\xi_{i})(1+X_{i}\beta)+\frac{\xi_{i}^{2}}{2}+\frac{1}{2}\beta^{T}X_{i}X_{i}^{T}\beta\right]$
$\displaystyle+\sum_{i=1}^{n}y_{i}X_{i}\beta-\sum_{i=1}^{n}\log y_{i}!$
$\displaystyle=\log\tilde{P}(y|\beta,\xi)$ (5)

To provide conjugacy in the VB component $q(\beta)$, we replace $P(y|\beta)$
by $\tilde{P}(y|\beta,\xi)$, in the computation of this element.

We denote this method by Laplace-VB. In Appendix A, we have proved that the
optimal posterior densities for the parameters are as follows

$q(\beta)={\cal N}(\beta;\mu_{\beta},\Sigma_{\beta})$

where ${\cal N}(x;\mu,\Sigma)$ stands for the multivariate normal density
function and

$\Sigma_{\beta}=\left[S_{X}^{\xi}+\text{diag}E(\tau^{-1})\right]^{-1},\quad\mu_{\beta}=\Sigma_{\beta}X^{T}(y-M_{\xi}),$
(6)

and thus, $E(\beta_{j}^{2})=\mu_{\beta_{j}}^{2}+\Sigma_{\beta_{j}}=d_{jj}$,

$q(\lambda)=\text{Gamma}(\lambda;p+\nu,\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})),$

and $E(\lambda)=\frac{p+\nu}{\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})}$,
for $j=1,\ldots,p$,

$q(\tau_{j})=\text{Gen-Inv-
Gauss}\Big{(}\tau_{j};\frac{1}{2},E(\lambda),E(\beta_{j}^{2})\Big{)},$

independently, where

$\text{Gen-Inv-
Gauss}\Big{(}x;p_{*},a_{*},b_{*},\Big{)}=\frac{(a_{*}/b_{*})^{p_{*}/2}}{2K_{p}(\sqrt{a_{*}b_{*}})}x^{p_{*}-1}e^{-(a_{*}x+b_{*}/x)/2},$

is the pdf of generalized inverse Gaussian distribution, where $K_{p}(\cdot)$
is the modified Bessel function of the second kind. Hence,

$E(\tau_{j})=\frac{\sqrt{E(\beta_{j}^{2})}K_{3/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}{\sqrt{E(\lambda)}K_{1/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})},$
$E(\tau_{j}^{-1})=\frac{\sqrt{E(\lambda)}K_{3/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}{\sqrt{E(\beta_{j}^{2})}K_{1/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}-\frac{1}{E(\beta_{j}^{2})},$

and

$E(\log\tau_{j})=\log\sqrt{\frac{E(\beta_{j}^{2})}{E(\lambda)}}+\frac{\partial}{\partial
p}\log K_{p}(\sqrt{E(\lambda)E(\beta_{j}^{2})})|_{p=1/2},$
$q(\tau_{0})=\text{Inv-Gamma}(\tau_{0};1,\frac{1}{2}d_{00}+E(a^{-1})),$

and thus $E(\tau_{0}^{-1})=\dfrac{1}{\frac{1}{2}d_{00}+E(a^{-1})}$ and
$E(\log\tau_{0})=-\psi(1)+\log(\frac{1}{2}d_{00}+E(a^{-1}))$,

$q(a)=\text{Inv-Gamma}(a;1,E(\tau_{0}^{-1})+\dfrac{1}{A}),$

and therefore $E(a^{-1})=\dfrac{1}{E(\tau_{0}^{-1})+\dfrac{1}{A}}$ and $E(\log
a)=-\psi(1)+\log(E(\tau_{0}^{-1})+\dfrac{1}{A}))$.

The ELBO is calculated as follows

$\displaystyle\log p(y;q)=$ $\displaystyle-
M_{\xi}^{T}(1+X_{i}\mu_{\beta})-\sum_{i=1}{n}\frac{\xi_{i}^{2}}{2}e^{\xi_{i}}-\frac{1}{2}tr(S_{X}^{\xi}D_{\beta})+y^{T}X\mu_{\beta}\allowdisplaybreaks$
$\displaystyle-\sum_{i=1}^{n}\log\Gamma(y_{i}+1)+y_{i}\sigma(C_{i})(1-\sigma(C_{i}))+\frac{1}{2}E(\log\tau_{0})\allowdisplaybreaks$
$\displaystyle-\frac{1}{2}\sum_{j=1}^{p}E(\log\tau_{j})-\frac{1}{2}E(\tau_{0}^{-1})E(\beta_{0}^{2})-\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}^{-1})E(\beta_{j}^{2})\allowdisplaybreaks$
$\displaystyle-\frac{1}{2}\log(\frac{1}{2}d_{00}+E(a^{-1}))+\frac{1}{4}\sum_{j=1}^{p}\dfrac{\mu_{j}^{\ast}}{\lambda_{j}^{\ast}}\allowdisplaybreaks$
$\displaystyle+(\nu-1)\left[\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))\right]-\delta\dfrac{p+\nu}{\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})}\allowdisplaybreaks$
$\displaystyle-\frac{1}{2}\log(\frac{1}{A}+E(\tau_{0}^{-1}))-\frac{3}{2}\log(\frac{1}{2}d_{00}+E(a^{-1}))-E(a^{-1})E(\tau_{0}^{-1})\allowdisplaybreaks$
$\displaystyle+\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))-\frac{1}{2}E(\lambda)E(\tau_{j})\allowdisplaybreaks$
$\displaystyle+\frac{3}{2}\Big{(}\psi(1)-\log\big{(}\frac{1}{A}+E(\tau_{0}^{-1})\big{)}\Big{)}-\dfrac{\frac{1}{A}}{\frac{1}{A}+E(\tau_{0}^{-1})}+\frac{1}{2}\log|\Sigma_{\beta}|\allowdisplaybreaks$
$\displaystyle-(p+\nu-1)\left[\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))\right]\allowdisplaybreaks$
$\displaystyle-(p+\nu)\log\big{(}\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})\big{)}+\log(\frac{1}{2}d_{00}+E(a^{-1}))\allowdisplaybreaks$
$\displaystyle-\frac{1}{2}E(\log\tau_{j})+\frac{1}{2}E(\tau_{j}^{-1})d_{jj}+\frac{1}{2}E(\lambda)E(\tau_{j})\allowdisplaybreaks$
$\displaystyle-\frac{1}{2}\log(\lambda_{j}^{\ast})+\log\big{(}E(\tau_{0}^{-1})+\frac{1}{A}\big{)}+\text{Const.}$
(7)

Algorithm 1 presents the algorithm for the Laplace-VB method.

Algorithm 1 Laplace-VB method for sparse Poisson regression model.

$\bullet$ Set proper initial values for the hyper-parameters of $q(\cdot)$
functions,

$\bullet$ Set $\epsilon$ equal to an arbitrary small value,

while The absolute relative change in ELBO is greater than $\epsilon$ do

$\bullet$ update
$\Sigma_{\beta}\leftarrow\left[S_{X}^{\xi}+\text{diag}E(\tau^{-1})\right]^{-1},$
and $\mu_{\beta}\leftarrow\Sigma_{\beta}X^{T}(y-M_{\xi}),$

$\bullet$ update
$E(\beta_{j}^{2})\leftarrow\mu_{\beta_{j}}^{2}+\Sigma_{\beta_{j}}=d_{jj}$,

$\bullet$ update
$E(\lambda)\leftarrow\frac{p+\nu}{\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})}$

$\bullet$ update
$E(\tau_{j})\leftarrow\frac{\sqrt{E(\beta_{j}^{2})}K_{3/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}{\sqrt{E(\lambda)}K_{1/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})},$

$\bullet$ update
$E(\tau_{j}^{-1})\leftarrow\frac{\sqrt{E(\lambda)}K_{3/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}{\sqrt{E(\beta_{j}^{2})}K_{1/2}(\sqrt{E(\lambda)E(\beta_{j}^{2})})}-\frac{1}{E(\beta_{j}^{2})},$

$\bullet$ update
$E(\tau_{0}^{-1})\leftarrow\dfrac{1}{\frac{1}{2}d_{00}+E(a^{-1})}$

$\bullet$ update
$E(a^{-1})\leftarrow\dfrac{1}{E(\tau_{0}^{-1})+\dfrac{1}{A}}$,

$\bullet$ calculate ELBO from (7)

end while

###  3.2 Continuous spike and slab prior

One of the most famous sparsity enforcing priors is the spike and slab prior
(Mitchell and Beauchamp, 1988; Geweke, J., 1996; George and McCulloch, 1997),
which is a mixture of a point mass $\delta$ and a continuous prior, which is
usually considered to be Gaussian. Let $Z=(Z_{1},\ldots,Z_{p})^{\top}$ be a
vector of the latent binary vector. Then, the spike and slab prior to the
vector of regression coefficient $\beta$ are defined as

$P(\beta)=\prod_{j=1}^{p}\left[Z_{j}\mathcal{N}(\beta_{j};0,\sigma^{2})+(1-Z_{j})\delta(\beta_{j})\right],$
(8)

where $\delta(\beta_{j})=1,$ if $\beta_{j}=0$ and $=0$, otherwise.

Continuous relaxations of (9), with $\delta(\cdot)$ replaced by a peaked
continuous density, is considered by many authors (see e.g., George and
McCulloch, 1993; Ishwaran and Rao, 2003, 2005; Ročková, 2018, among others).
One of the most famous variations of the continuous spike and slab prior is as
follows

$P(\beta)=\prod_{j=1}^{p}\left[Z_{j}\mathcal{N}(\beta_{j};0,\sigma^{2})+(1-Z_{j})\mathcal{N}(\beta_{j};0,c\sigma^{2})\right],$
(9)

in which the constant $c$ is small enough to enforce sparsity to $\beta_{j}$.
Here, we take $c=0.001$.

Thus, the model is assumed to be as follows

$\displaystyle y_{i}\sim$
$\displaystyle\text{Poiss}(\lambda_{i}=\exp(X_{i}\beta)),\quad i=1,\ldots,n$
$\displaystyle p(\beta)=$
$\displaystyle\prod_{j=0}^{p}\left[Z_{j}\mathcal{N}(0,\sigma^{2})+(1-Z_{j})\mathcal{N}(0,c\sigma^{2})\right],$
$\displaystyle Z_{0}=1,$ $\displaystyle{\rm w.p.1}$ $\displaystyle Z_{j}\sim$
$\displaystyle\text{Ber}(\pi_{j}),\;\;\pi_{j}\sim\text{Beta}(\theta_{1},\theta_{2})\quad
j=1,\ldots,p,$ $\displaystyle\sigma^{2}\sim$ $\displaystyle\text{Inv-
Gamma}(\dfrac{1}{2},\dfrac{1}{a}),\;\;a\sim\text{Inv-
Gamma}\left(\dfrac{1}{2},\dfrac{1}{A}\right).$ (10)

The corresponding DAG for model (10) is presented in Figure 2.

Figure 2: The directed acyclic graph (DAG) for model (10). Data elements are
shown by squares, model parameters are shown by circles and fixed hyper-
parameters are shown by polygons.

Again, we use the approximation given in (5) to provide conjugacy in VB
components. Here, we let $\xi_{i}=X_{i}\mu_{\beta}$, for $i=1,\ldots,n$, with
$\mu_{\beta}$ is given in (11).

This method is denoted by CS-VB. The optimal posterior densities for the
parameters are computed in Appendix A and are as follows

$q(\beta)={\cal N}(\beta;\mu_{\beta},\Sigma_{\beta})$

where

$\Sigma_{\beta}=\left[S_{X}^{\xi}+E(\sigma^{-2})\left[\text{diag}(P)+c^{-1}(\text{I}-\text{diag}(P))\right]\right]^{-1}$

and

$\mu_{\beta}=\Sigma_{\beta}X^{T}(y-M_{\xi}),$ (11) $q(\sigma^{2})=\text{Inv-
Gamma}\left(\sigma^{2};\dfrac{p-1}{2},\frac{1}{2}\text{tr(diag}(P)D_{\beta})+\dfrac{1}{2c}\text{tr(diag}(1-P)D_{\beta})+E(a^{-1})\right),$
$q(a)=\text{Inv-Gamma}\left(a;1,E(\sigma^{-2})+\dfrac{1}{A}\right),$
$q(Z_{j})=\text{Ber}(Z_{j};P_{j}),\quad j=1,\ldots,p,$

where
$P_{j}=\sigma\left(\psi(\theta_{1}+P_{j})-\psi(\theta_{2}-P_{j}+1)-\frac{1}{2}E(\sigma^{-2})d_{jj}(1-\frac{1}{c})\right),$
and

$q(\pi_{j})=\text{Ber}(\pi_{j};P_{j}+\theta_{1},\theta_{2}-P_{j}+1).$

The ELBO is also obtained as follows

$\displaystyle\log p(y;q)=$ $\displaystyle-
M_{\xi}X\mu_{\beta}-\sum_{i=1}{n}(M_{\xi_{i}}-e^{\xi_{i}}\frac{\xi_{i}}{2})^{2}-\frac{1}{2}tr(S_{X}^{\xi}D_{\beta})+y^{T}X\mu_{\beta}\allowdisplaybreaks$
$\displaystyle-\frac{p}{2}\log
2\pi-\frac{1}{2}\log(1-\sum_{j=1}^{p}P_{j})-\frac{1}{2}(1-\frac{1}{c})E(\sigma^{-2})\sum_{j=1}^{p}P_{j}d_{jj}-\frac{1}{2c}E(\sigma^{-2})\sum_{j=1}^{p}d_{jj}\allowdisplaybreaks$
$\displaystyle-\frac{p}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}})\frac{1}{2}(\psi(1)-\log(\frac{1}{A}+E(\sigma^{-2})))-\log\Gamma(\frac{1}{2})\allowdisplaybreaks$
$\displaystyle-\frac{3}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}})-\dfrac{E(\sigma^{-2})-\frac{1}{A}}{\frac{1}{A}+E(\sigma^{-2})}+\sum_{j=1}^{p}P_{j}(\psi(\alpha^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))\allowdisplaybreaks$
$\displaystyle+\sum_{j=1}^{p}(1-P_{j})(\psi(\beta^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))\sum_{j=1}^{p}\log\dfrac{\Gamma(\theta_{1})\Gamma(\theta_{2})}{\Gamma(\theta_{1}+\theta_{2})}\allowdisplaybreaks$
$\displaystyle+(\theta_{1}-1)\sum_{j=1}^{p}(\psi(\alpha^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))+(\theta_{2}-1)\sum_{j=1}^{p}(\psi(\beta^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))\allowdisplaybreaks$
$\displaystyle+\frac{3}{2}\Big{(}\psi(1)-\log\big{(}\frac{1}{A}+E(\sigma^{-2})\big{)}\Big{)}+\frac{1}{2}\log|\Sigma_{\beta}|-\alpha_{\sigma^{2}}\log\beta_{\sigma^{2}}\allowdisplaybreaks$
$\displaystyle+\log\Gamma(\alpha_{\sigma^{2}})+(\alpha_{\sigma^{2}}-1)\left[\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}}\right]+\alpha_{\sigma^{2}}-E\Big{[}\sum_{j=1}^{p}P_{j}\log(P_{j})\allowdisplaybreaks$
$\displaystyle+(1-P_{j})\log(1-P_{j})\Big{]}+\log\Big{(}\frac{1}{A}+E(\sigma^{-2})\Big{)}+\text{Const.}$
(12)

Algorithm 2 presents the algorithm for the Laplace-VB method.

Algorithm 2 CS-VB method for sparse Poisson regression model.

$\bullet$ Set proper initial values for the hyper-parameters of $q(\cdot)$
functions,

$\bullet$ Set $\epsilon$ equal to an arbitrary small value,

while The absolute relative change in ELBO is greater than $\epsilon$ do

$\bullet$ update
$\Sigma_{\beta}\leftarrow\left[S_{X}^{\xi}+E(\sigma^{-2})\left[\text{diag}(P)+c^{-1}(\text{I}-\text{diag}(P))\right]\right]^{-1}$

$\bullet$ update $\mu_{\beta}\leftarrow\Sigma_{\beta}X^{T}(y-M_{\xi}),$

$\bullet$ update
$E(\sigma^{-2})\leftarrow{\dfrac{p-1}{2}}\left(\frac{1}{2}\text{tr(diag}(P)D_{\beta})+\dfrac{1}{2c}\text{tr(diag}(1-P)D_{\beta})+E(a^{-1})\right)^{-1},$

$\bullet$ update $E(a^{-1})=\left(E(\sigma^{-2})+\dfrac{1}{A}\right)^{-1},$

$\bullet$ update
$P_{j}\leftarrow\sigma\left(\psi(\theta_{1}+P_{j})-\psi(\theta_{2}-P_{j}+1)-\frac{1}{2}E(\sigma^{-2})d_{jj}(1-\frac{1}{c})\right),$

$\bullet$ calculate ELBO from (12)

end while

###  3.3 Bernoulli sparsity enforcing prior

To enforce sparsity to the regression coefficient vector $\beta$, Zhang et al.
(2019) used the product of a diagonal matrix $\Gamma$ and the vector $\beta$,
in which the diagonal elements of $\Gamma$ are Bernoulli distributed. When the
$j$th diagonal element of $\Gamma$, is zero, the effect of $\beta_{j}$ is
removed from the linear regression function. Zhang et al. (2019) used this
prior for the logistic regression model.

Thus, the Bayesian sparse Poisson model with Bernoulli prior is as follows

$\displaystyle y_{i}\sim$
$\displaystyle\text{Poiss}(\lambda_{i}=\exp(X_{i}\Gamma\beta)),\quad
i=1,\ldots,n$ $\displaystyle X_{i}=$
$\displaystyle(1,X_{i1},\ldots,X_{ip})^{T},\quad\Gamma=\text{diag}(1,\gamma_{1},\ldots,\gamma_{p})$
$\displaystyle\beta\sim$ $\displaystyle\;\mathcal{N}_{p+1}(0,A^{-1})$
$\displaystyle A=$
$\displaystyle\;\text{diag}(\alpha_{0},\alpha_{1},\ldots,\alpha_{p}),\;\;\alpha_{j}\sim\text{Gamma}(a_{j},b_{j}),\;j=0,\ldots,p$
$\displaystyle\gamma_{j}\sim$
$\displaystyle\text{Ber}(\pi_{j}),\;\;\pi_{j}\sim\text{Beta}(\theta_{1},\theta_{2})\quad
j=1,\ldots,p,$ (13)

The corresponding DAG for model (13) is presented in Figure 3.

Figure 3: The directed acyclic graph (DAG) for model (13). Data elements are
shown by squares, model parameters are shown by circles and fixed hyper-
parameters are shown by polygons.

Similar to the approximation given in (5), we find an upper bound for the
likelihood function as follows to provide conjugacy in VB components

$\displaystyle\log P(y_{i};\beta,\gamma)=$
$\displaystyle-\exp(X_{i}\Gamma\beta)+y_{i}X_{i}\Gamma\beta-\log y_{i}!$
$\displaystyle\log P(y_{i}|\beta,\gamma)\leq$
$\displaystyle-e^{\xi_{i}}\Big{[}(1-\xi_{i})(1-X_{i}\Gamma\beta)+\dfrac{\xi_{i}^{2}}{2}$
$\displaystyle+\frac{1}{2}\beta^{T}\Gamma
X_{i}^{T}X_{i}\Gamma\beta\Big{]}+y_{i}X_{i}\Gamma\beta+\log y_{i}!$
$\displaystyle=$ $\displaystyle\log\tilde{P}(y_{i}|\beta,\gamma,\xi_{i})$

Here, we let $\xi_{i}=X_{i}{\rm diag}(P)\mu_{\beta}$, for $i=1,\ldots,n$, with
$\mu_{\beta}$ and $P$ given in (14) and (15), respectively.

We denote this method by Bernoulli-VB. In Appendix A, we prove that the
optimal posterior densities for the parameters are as follows

$q(\beta)=\mathcal{N}(\beta;\mu_{\beta},\Sigma_{\beta})$

where

$\displaystyle\Sigma_{\beta}=\left[S_{X}^{\xi}\odot\Omega+\text{diag}\Big{(}E(\alpha)\Big{)}\right]^{-1},$
$\mu_{\beta}=\Sigma_{\beta}\text{diag}(P)X^{T}(y-M_{\xi})$ (14)

in which $\Omega=PP^{T}+\text{diag}(P)(I_{p+1}-\text{diag}(P)),$

$\displaystyle q(\gamma_{j})=$ $\displaystyle\text{Ber}(\gamma_{j};P_{j})$

where

$P_{j}=\sigma\left((y-M_{\xi})^{T}X_{j}\beta_{j}-\frac{1}{2}S_{X_{jj}}^{\xi}d_{jj}-\frac{1}{2}\sum_{i\neq
j}P_{i}S_{X_{ij}}^{\xi}d_{ij}+E(\log\pi_{j})-E(\log(1-\pi_{j}))\right),$ (15)

and $P=(1,P_{1},\ldots,P_{p}),$

$q(\alpha_{j})=\text{Gamma}\left(\alpha_{j};a_{j}+1,b_{j}+\frac{d_{jj}}{2}\right),$

thus, $E(\alpha_{j})={(a_{j}+1)}/{(b_{j}+{d_{jj}}/{2})},$ and

$q(\pi_{j})=\text{Beta}(\pi_{j};P_{j}+\theta_{1},\theta_{2}-P_{j}+1),$

and hence,
$E(\log\pi_{j})=\psi(P_{j}+\theta_{1})-\psi(\theta_{1}+\theta_{2}+1),$ and
$E(\log(1-\pi_{j}))=\psi(\theta_{2}-P_{j}+1)-\psi(\theta_{1}+\theta_{2}+1).$

Furthermore, ELBO is obtained as follows

$\displaystyle\log p(y;q)=$ $\displaystyle-
M_{\xi}XP\mu_{\beta}-\frac{1}{2}\text{tr}(\Sigma_{\beta}+\mu_{\beta}\mu_{\beta}^{T})(S_{X}^{\xi}\odot\Omega)+\mu_{\beta}^{T}\Sigma_{\beta}^{-1}\mu_{\beta}\allowdisplaybreaks$
$\displaystyle-\sum_{i=1}^{n}\big{(}e^{\xi_{i}}+\xi_{i}-\frac{\xi_{i}^{2}}{2}\big{)}+\frac{1}{2}\Big{[}\sum_{j=1}^{p}\psi(\alpha_{j})-\sum_{j=1}^{p}\log
b_{j}\allowdisplaybreaks$
$\displaystyle-\text{tr}\big{[}(\Sigma_{\beta}+\mu_{\beta}\mu_{\beta}^{T})\text{diag}(\frac{a_{j}}{b_{j}})\big{]}\Big{]}+\sum_{j=1}^{p}\big{[}P_{j}\log\pi_{j}+(1-P_{j})\allowdisplaybreaks$
$\displaystyle\log(1-\pi_{j})\big{]}+\sum_{j=1}^{p}(\alpha_{0}+1)\left[\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right]\allowdisplaybreaks$
$\displaystyle+\sum_{j=1}^{p}(\beta_{0}-1)\left[\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right]+\frac{1}{2}\log\big{(}\det(\Sigma_{\beta})\big{)}\allowdisplaybreaks$
$\displaystyle-\sum_{j=1}^{p}P_{j}(\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast}))\allowdisplaybreaks$
$\displaystyle-\sum_{j=1}^{p}(1-P_{j})(\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast}))-\sum_{j=1}^{p}\log\Gamma(\alpha_{j}^{\ast}-\beta_{j}^{\ast})\allowdisplaybreaks$
$\displaystyle+\sum_{j=1}^{p}P_{j}\log
P_{j}+(1-P_{j})\log(1-P_{j})\allowdisplaybreaks$
$\displaystyle+\log\Gamma(\alpha_{j}^{\ast})+\log\Gamma(\beta_{j}^{\ast})-(\alpha_{j}^{\ast}-1)(\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast}))\allowdisplaybreaks$
$\displaystyle-(\beta_{j}^{\ast}-1)(\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast}))+\sum_{j=1}^{p}\big{[}(a_{0}-1)(\psi(a_{0})\allowdisplaybreaks$
$\displaystyle-\log(b_{j}+\frac{d_{jj}}{2}))-b_{0}\frac{a_{j}}{b_{j}+\frac{d_{jj}}{2}}\big{]}+\text{Const.}$
(16)

Algorithm 3 presents the algorithm for the Laplace-VB method.

Algorithm 3 Bernoulli-VB method for sparse Poisson regression model.

$\bullet$ Set proper initial values for the hyper-parameters of $q(\cdot)$
functions,

$\bullet$ Set $\epsilon$ equal to an arbitrary small value,

while The absolute relative change in ELBO is greater than $\epsilon$ do

$\bullet$ update
$\Sigma_{\beta}\leftarrow\left[S_{X}^{\xi}\odot\Omega+\text{diag}\Big{(}E(\alpha)\Big{)}\right]^{-1},$
$\Omega=PP^{T}+\text{diag}(P)(I_{p+1}-\text{diag}(P)),$

$\bullet$ update
$\mu_{\beta}\leftarrow\Sigma_{\beta}\text{diag}(P)X^{T}(y-M_{\xi}),$

$\bullet$ update $E(\alpha_{j})={(a_{j}+1)}/{(b_{j}+{d_{jj}}/{2})},$

$\bullet$ update
$E(\log\pi_{j})=\psi(P_{j}+\theta_{1})-\psi(\theta_{1}+\theta_{2}+1),$ and
$E(\log(1-\pi_{j}))=\psi(\theta_{2}-P_{j}+1)-\psi(\theta_{1}+\theta_{2}+1),$

$\bullet$ update
$P_{j}=\sigma\left((y-M_{\xi})^{T}X_{j}\beta_{j}-\frac{1}{2}S_{X_{jj}}^{\xi}d_{jj}-\frac{1}{2}\sum_{i\neq
j}P_{i}S_{X_{ij}}^{\xi}d_{ij}+E(\log\pi_{j})-E(\log(1-\pi_{j}))\right),$

$\bullet$ calculate ELBO from (16)

end while

##  4 Posterior predictive mass function

The posterior predictive mass function (ppmf) of the response variable
associated with the covariates of a new sample $x_{0}$ given the train data
set $(X,y)$ is obtained as follows

$\displaystyle p(y_{0}|x_{0},X,y)$ $\displaystyle=\int
p(y_{0}|x_{0},\theta)p(\theta|X,y)d\theta$
$\displaystyle\approx\int{p}(y_{0}|x_{0},\theta)q(\theta)d\theta,$

where $\theta=\beta$ for Laplace-VB and CS-VB methods, and
$\theta=(\gamma,\beta)^{\top}$ for Bernoulli-VB method.

For Laplace-VB and CS-VB methods, we have

$p(y_{0}|x_{0},X,y)\approx\frac{1}{y_{0}!}\int\exp\\{-\exp\\{x_{0}\beta\\}\\}\exp\\{y_{0}x_{0}\beta\\}{\cal
N}(\beta;\mu_{\beta},\Sigma_{\beta})\;d\beta$

Using re-paramerization $\lambda_{0}=\exp\\{x_{0}\beta\\}$, we see that

$\displaystyle p(y_{0}|x_{0},X,y)$ $\displaystyle\approx\frac{1}{y_{0}!}\int
e^{-\lambda_{0}}\lambda_{0}^{y_{0}}{\cal
LN}(\lambda_{0};x_{0}\mu_{\beta},x_{0}\Sigma_{\beta}x_{0}^{\top})\;d\lambda_{0}.$
(17)

The predictive probability mass function in (17) can be obtained using a uni-
variate numerical integration, and might be used for obtaining a point
predictor (e.g. most likely predictor) and a prediction interval (e.g. the
most likely prediction interval) for $y_{0}$ given $x_{0}$ and $(X,y)$.

For Bernoulli-VB method, we have

$p(y_{0}|x_{0},X,y)\approx\prod_{j=1}^{p}\sum_{\gamma_{j}=0}^{1}{\rm
Ber}(\gamma_{j};P_{j})\frac{1}{y_{0}!}\int\exp\\{-\exp\\{x_{0}\Gamma\beta\\}\\}\exp\\{y_{0}x_{0}\Gamma\beta\\}{\cal
N}(\beta;\mu_{\beta},\Sigma_{\beta})\;d\beta$

Using re-paramerization $\lambda_{0\gamma}=\exp\\{x_{0}\Gamma\beta\\}$, we see
that

$\displaystyle p(y_{0}|x_{0},X,y)$
$\displaystyle\approx\prod_{j=1}^{p}\sum_{\gamma_{j}=0}^{1}{\rm
Ber}(\gamma_{j};P_{j})\frac{1}{y_{0}!}\int
e^{-\lambda_{0\gamma}}\lambda_{0\gamma}^{y_{0}}{\cal
LN}(\lambda_{0\gamma};x_{0}\Gamma\mu_{\beta},x_{0}\Gamma\Sigma_{\beta}\Gamma
x_{0}^{\top})\;d\lambda_{0\gamma}.$ (18)

Since, for the Bernoulli-VB method, the parameters $P_{j}$ take only 0 and 1
values, the formulla (18) can be simply written as

$\displaystyle p(y_{0}|x_{0},X,y)$ $\displaystyle\approx\frac{1}{y_{0}!}\int
e^{-\lambda_{0}}\lambda_{0}^{y_{0}}{\cal
LN}(\lambda_{0};x_{0}P\mu_{\beta},x_{0}P\Sigma_{\beta}Px_{0}^{\top})\;d\lambda_{0}.$
(19)

##  5 Simulation study

In order to examine the performance of the proposed models, we have conducted
a simulation study as follows. In each of $N=1000$ replications, we have
generated samples of size $n=30$, form a sparse Poisson regression model,
whose vector of coefficients $\beta_{p\times 1}$, with $p=5$, is first
generated from $N(1,1)$ and then 30 percent of randomly chosen elements of
$\beta$ are set to zero. The rows of the covariate matrix $X$ are generated
independently from a multivariate normal distribution with zero mean and the
variance-covariance matrix with elements $\Sigma_{ij}=0.3^{|i-j|}$.

Three proposed VB models (Laplace, CS, Bernoulli) are fitted to each generated
data set. Furthermore, two frequentist alternatives, LASSO and SCAD-penalized
Poisson regression models are considered for comparison. The sparsity
parameters of both LASSO and SCAD models are optimized based on the Corrected
Akaike Information Criterion

${\rm AICc}=-\log\widehat{L}(\beta|x)+2{\rm df}+2{\rm df}({\rm df}+1)/(n-{\rm
df}-1),$

where $\widehat{L}(\beta|x)$ is the estimated likelihood of the model and
${\rm df}$ is the degree of the freedom of the model.

For the CS-VB and Laplace-VB models, the $j$th covariate is selected, if
$J_{j}=1$, where

$J_{j}=I\left(\frac{P_{j}}{\sqrt{P_{j}(1-P_{j})}}>z_{1-0.05/p}\right),\;j=1,\ldots,p,$
(20)

for the CS-VB method, and

$J_{j}=I\left(\frac{\mu_{\beta_{j}}}{\sqrt{\Sigma_{\beta_{jj}}}}>z_{1-0.05/p}\right),\;j=1,\ldots,p,$
(21)

for Laplace-VB model, where $z_{\alpha}$ stands for the $\alpha$th quantile of
the standard normal distribution.

To compare the performance of the 3 proposed models with each other and with
LASSO and SCAD-penalized Poisson regression models, the regression coefficient
relative error is computed as follows

${\rm
CRE}=\frac{\frac{1}{N}\sum_{t=1}^{N}(\widehat{\beta}-\beta)^{\top}(\widehat{\beta}-\beta)}{\frac{1}{N}\sum_{t=1}^{N}\beta^{\top}\beta}.$

Furthermore, the train-set and test-set relative errors are obtained,
respectively, as

${\rm TRRE}=\frac{\frac{1}{N}\sum_{t=1}^{N}(\widehat{y}_{\rm train}-{y}_{\rm
train})^{\top}(\widehat{y}_{\rm train}-{y}_{\rm
train})}{\frac{1}{N}\sum_{t=1}^{N}{y}_{\rm train}^{\top}{y}_{\rm train}},$

and

${\rm TSRE}=\frac{\frac{1}{N}\sum_{t=1}^{N}(\widehat{y}_{\rm test}-{y}_{\rm
test})^{\top}(\widehat{y}_{\rm test}-{y}_{\rm
test})}{\frac{1}{N}\sum_{t=1}^{N}{y}_{\rm test}^{\top}{y}_{\rm test}},$

where ${y}_{\rm train}=({y}_{1},\ldots,{y}_{n})^{\top}$ and ${y}_{\rm
test}=({y}^{\prime}_{1},\ldots,{y}^{\prime}_{n^{\prime}})^{\top}$, are the
train-set (80% of the full sample) and test-set (20 % of the full sample)
response samples, respectively.

The false negative rate and the false positive rate are also defined as
follows

${\rm FNR}=\frac{\\#\\{j;\;1\leq j\leq p,\;\beta_{j}\neq
0,\;\hat{\beta}_{j}=0\\}}{\\#\\{j;\;1\leq j\leq p,\;\beta_{j}\neq 0\\}},$
${\rm FPR}=\frac{\\#\\{j;\;1\leq j\leq p\;\beta_{j}=0,\;\hat{\beta}_{j}\neq
0\\}}{\\#\\{j;\;1\leq j\leq p,\;\beta_{j}=0\\}},$

where $\\#A$ stands for the cardinality of the set $A$. These criteria are
used to examine the sparsing performance of the competitive methods.

The computation time of 3 proposed models as well as the LASSO and SCAD-
penalized Poisson regression models are also computed in each iteration.

For the Laplace-VB method, we set the hyper-priors $\nu=0.001$, and
$\delta=0.001$. The hyper-parameters $\theta_{1}$ and $\theta_{2}$ for the CS-
VB and Bernoulli-VB methods are both set to $3/2$. Also all $a_{j}$ and
$b_{j}$, hyper-parameters for the Bernoulli-VB method are set to 0.01. The
hyper-parameter $A$ is set to 0.01 for the Laplace-VB and CS-VB methods.

Figures 5 and 6 show the box plots of all aforementioned criteria for all
competitive methods. The Bernoulli-VB and CS-VB methods are the best methods
with respect to the CRE criterion. The Laplace-VB method performs as well as
the LASSO method. Based on the TRRE, the Bernoulli-VB method is as good as the
LASSO, while the CS-VB and Laplace-VB perform better than the SCAD method. The
TSRE for all three proposed VB methods are better than LASSO and SCAD. Note
that for the computation of TSRE values for the LASSO and SCAD methods, the
prediction of the response is computed as

$\hat{y}_{i}=\exp(X_{\rm test}\hat{\beta}_{\rm sparse}),$

in which $\hat{\beta}_{\rm sparse}$ is the sparse estimator obtained by the
LASSO and SCAD methods, while the TSRE values for the VB methods are computed
using the ppmf functions (17) for the CS-VB and Laplace-VB and the ppmf
function (19) for the Bernoulli-VB method. The mode of (17) is considered as
the point predictor of the response value. Furthermore, it can be seen from
Figure 5 that the TSRE values are smaller for Laplace-VB compared to CS-VB and
Bernoulli-VB, on average.

Considering the FPR, the CS-VB method is the best method, while the Bernoulli-
VB fails. The FPR of the Laplace-VB is a little better than the LASSO method.
The FNR of the CS-VB and Laplace-VB are as good as SCAD, while the Bernoulli-
VB is comparable with LASSO. Although the computation time of the VB methods
is a little more than the LASSO and SCAD methods, this comparison is not fair,
since the LASSO and SCAD methods use the C++ programming for accelerating the
computation, while we have done all computations in the R software. So, it
seems that the VB methods are very fast, as expected.

Table 1 presents the average coverage probabilities of the highest posterior
density (HPD) confidence intervals for the regression coefficients for three
VB methods at level 0.95. As it can be seen from Table 1, the average coverage
probabilities of all regression coefficients, except the intercept, are
greater than the nominal value. A sample of HPD intervals in a single
iteration of the simulation study, for a zero coefficient, based on three VB
methods is given in Figure 4.

Figure 4: A sample of HPD intervals in a single iteration of the simulation study, for a zero coefficient, based on Bernoulli-VB (left) CS-VB (middle) and Laplace-VB (right) methods. Table 1: The average coverage probabilities of the HPD intervals for the regression coefficients for three VB methods at level 0.95. | Intercept | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$
---|---|---|---|---|---
Bernoulli | 0.92 | 0.99 | 0.97 | 0.97 | 0.98
CS | 0.95 | 0.95 | 0.97 | 0.96 | 0.99
Laplace | 0.88 | 0.98 | 0.94 | 0.97 | 0.97
Table 2: Accuracy values for the parameters and hyper-parameters of the VB against an MCMC benchmark for different priors. | Parameters |  |  |  |
---|---|---|---|---|---
Method | Intercept | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$
Bernoulli | 95.82 | 96.18 | 94.96 | 29.54 | 31.60
CS | 96.34 | 94.42 | 93.59 | 64.20 | 74.32
Laplace | 97.27 | 97.42 | 97.12 | 90.42 | 92.70
Method |  | $\gamma_{1}\;(z_{1})$ | $\gamma_{2}\;(z_{2})$ | $\gamma_{3}\;(z_{3})$ | $\gamma_{4}\;(z_{4})$
Bernoulli |  | 100 | 100 | 53.80 | 66.40
CS |  | 100 | 100 | 58.03 | 52.15
Method |  | $\pi_{1}$ | $\pi_{2}$ | $\pi_{3}$ | $\pi_{4}$
Bernoulli |  | 98.37 | 98.31 | 78.59 | 84.11
CS |  | 95.24 | 97.10 | 81.75 | 77.51
Method | $\alpha_{0}$ | $\alpha_{1}$ | $\alpha_{2}$ | $\alpha_{3}$ | $\alpha_{4}$
Bernoulli | 84.14 | 81.37 | 81.85 | 33.92 | 43.8
Method | $\sigma^{2}$ | $\tau_{1}$ | $\tau_{2}$ | $\tau_{3}$ | $\tau_{4}$
Laplace | 30.56 | 60.07 | 63.14 | 66.3 | 65.27

Figure 5: The simulation results: the coefficient relative error (top), the
train relative error (middle), and the test relative error (bottom) for
different methods.

Figure 6: The simulation results: the FPR (top), the FNR (middle), and the
computation time (bottom) for different methods.

Figure 7: Approximate posterior density functions of the first three
regression coefficients for different models.

Figure 8: Approximate posterior density functions of the last two regression
coefficients for different models.

In order to examine the accuracy of the approximation of the posterior
distribution by VB methods, the accuracy measure (see Luts and Wand, 2015) is
computed for each parameter component $\theta_{j},\;j=1,\ldots,k,$ as follows

${\rm
accuracy}(q)=100\left(1-\frac{1}{2}\int_{-\infty}^{\infty}|q(\theta_{j})-p(\theta_{j}|y)|d\theta\right)\%,$

in a single iteration of the simulation study, where $p(\theta_{j}|y)$ is the
kernel density estimator using the MCMC sample. For the Bernoulli components
$q(\gamma_{j})$ and $q(z_{j})$ for Bernoulli-VB and CS-VB methods, the
accuracy is computed as follows

${\rm
accuracy}(q)=100\left(1-\frac{1}{2}\left(\left|q(\theta_{j}=1)-\frac{1}{T}\sum_{t=1}^{T}I(\theta_{jt}=1)\right|+\left|q(\theta_{j}=0)-\frac{1}{T}\sum_{t=1}^{T}I(\theta_{jt}=0)\right|\right)\right)\%,$

where $\theta_{j1},\ldots,\theta_{jT}$ are MCMC samples. The accuracy values
are given in Table 2. Furthermore, the $q$-functions as well as the kernel
density estimator of the marginal posteriors of the regression coefficients
are given in Figures 7, 8, and the corresponding plots for the other
parameters are given in Figures 11-16 of the Appendix-B. For the MCMC methods,
we have used the packages R2OpenBUGS and rjags in R software version 4.2.1,
OpenBUGS version 3.2.3, and JAGS version 4.3.1.

As one can see from Table 2, the accuracy of the posterior approximation for
the regression coefficients is high, except for $\beta_{3}$ by Bernoulli-VB
and CS-VB methods, and for $\beta_{4}$ for Bernoulli-VB method. Both these
coefficients have zero values (see Figure 8), and both Bernoulli-VB and CS-VB
methods under-estimate the posterior densities in such a situation, as one can
see from Figure 8. Indeed, for the zero coefficients, the Laplace-VB method
provides the most accurate posterior approximation. However, in general,
Bernoulli-VB and CS-VB have lower coefficient relative error than Laplace-VB
(see Figure 5) and the FPR-FNR performance of the CS-VB is better than that of
the Laplace-VB (see Figure 6). Table 2, also presents the accuracy for the
parameters $\gamma$ and $z$ of the VB against an MCMC benchmark for Bernoulli
and CS priors. As one can see from these values, the accuracy for the non-zero
coefficients is 100%, while the accuracy of the zero coefficients is low for
both Bernoulli-VB and CS-VB methods. The accuracy for the $\pi$ hyper-
parameter is high for both Bernoulli-VB and CS-VB methods. Regarding the
variance hyper-parameters of the Bernoulli-VB and Laplace-VB methods, we can
see from Tables 2 that for the Bernoulli-VB method, the low accuracies are for
the variances associated with zero coefficients, while the accuracies of the
variance hyper-parameters for the Laplace-VB method are all low.

##  6 Benchmark real data analysis

In order to examine the prediction performance of the proposed VB methods, we
have considered 6 benchmark real data sets with a count response variable. The
data set affairs (Fair, 1978) which is stored in a data frame in the R package
COUNT, is a data set with 601 observations and 18 variables, including,
naffairs (the number of affairs within the last year), and 17 covariates about
children, marriage happiness, being religious, and age levels. Greene (2003),
modeled this data using Poisson regression, although given the amount of over-
dispersion in the data, employing a negative binomial model is an appropriate
strategy. The variable naffairs is considered as the response variable. The
top left graph of Figure 9 shows the histogram of the response variable
naffairs. The bike-sharing data set contains the hourly and daily count of
rental bikes between the years 2011 and 2012 in the Capital bike-share system
with the corresponding weather and seasonal information. This data set
includes 731 observations and 14 variables (except date and instant number),
including cnt (the response variable), and 13 covariates including season,
year, month, holiday, weekday, working day, 5 weather variables, casual, and
registered. The top middle graph in Figure 9 presents the histogram of the
response variable for the bike-sharing data set. The fishing data (Zuur et
al., 2013) are originally adapted from Bailey et al. (2009), who were
interested in how certain deep-sea fish populations were impacted when
commercial fishing began in locations with deeper water than in previous
years. The observations are recorded from 147 sites for 7 variables including
totabund (total fish counted per site) and 6 covariates including depth, area,
foliage density, catch site, year, and period. Here, we have only considered 3
covariates density, mean-depth, and swept-area. The top right plot in Figure 9
shows the histogram for this data set. The three remaining data sets,
azcabgptca, azdrg112, and azpro, from package COUNT, are related to samples
from 1991 (azcabgptca and azpro) and 1995 (azdrg112) Arizona hospital
cardiovascular patient files, to receive one of two standard cardiovascular
procedures: CABG and PTCA. The data set azcabgptca has 1959 observations on 6
variables, azdrg112 has 1,798 observations on 4 variables and azpro has 3589
observations on 6 variables. The variable los (length of hospital stay) is the
response variable in all three data sets, and the covariates are procedure
type, sex, age, and some other covariates. The bottom plots of Figure 9 show
the histograms of the los response variable in azcabgptca, azdrg112, and
azpro, respectively, from left to right.

Table 3: Test relative error means (standard deviation) for 10 random partitions of 6 benchmark real data sets. Data set | LASSO | SCAD | Bernoulli-VB | CS-VB | LAPLACE-VB
---|---|---|---|---|---
affairs | 0.865 (0.071) | 0.868 (0.078) | 0.878 (0.065) | 0.875 (0.065) | 0.882 (0.068)
bike sharing | 0.053 (0.005) | 0.057 (0.007) | 0.058 (0.008) | 0.059 (0.009) | 0.058 (0.009)
fishing | 0.512 (0.608) | 0.519 (0.624) | 0.518 (0.621) | 0.570 (0.750) | 0.517 (0.620)
azcabgptca | 0.534 (0.049) | 0.535 (0.049) | 0.546 (0.052) | 0.547 (0.045) | 0.546 (0.052)
azdrg112 | 0.851 (0.038) | 0.851 (0.038) | 0.854 (0.031) | 0.858 (0.029) | 0.859 (0.029)
azpro | 0.624 (0.048) | 0.624 (0.048) | 0.626 (0.047) | 0.626 (0.047) | 0.626 (0.047)

We randomly partition the observations of each data set to train set (80 %)
and test set (20 %) and replicate the random partitioning 10 times to compute
the test set relative prediction error and compare it for all considered
competitors. Table 3 presents the means and standard deviation of test set
relative errors for 10 random partitions of all data sets, for three VB
methods, obtained from ppmfs in (17) and (19), as well as those for LASSO and
SCAD Poisson regression methods. As seen from Table 3, the proposed VB methods
perform as well as the frequentist sparse Poisson regression methods for all 6
benchmark data sets.

Figure 9: Histograms of the responses for the benchmark real data sets. Top
left; affairs, top middle; bike sharing, top right; fishing, bottom left;
azcabgptca, bottom middle; azdrg112, bottom right; azpro.

As an example, the posterior approximated densities of the regression
coefficients, and the corrresponding HPD intervals, for the fishing data set,
are plotted in Figure 10.

Figure 10: The posterior approximated densities of the regression
coefficients, and the corrresponding HPD intervals, for the fishing data set.
Columns from left to right are associated with Bernoulli-VB, CS-VB, and
Laplace-VB, and rows from top to bottom are associated with the intercept and
3 regression coefficients density, mean-depth, and swept-area.

##  7 Concluding remarks

The VB methods proposed in this paper provide a fast Bayesian inference about
the parameters of the sparse Poisson regression model. These methods are
significantly faster than the MCMC methods. The proposed VB methods are used
for simultaneous estimation and variable selection in the Poisson regression
model. It is observed that the proposed Bayesian inference performs better
than the two well-known frequentist methods, LASSO and SCAD penalized Poisson
regression models, in general, while there are some advantages of the Bayesian
inference, compared to the frequentist ones, including the assessment of the
highest posterior density intervals, posterior distributions of the
parameters, and the predictive probability mass function of the response,
given observations. It is observed that the proposed VB methods provide a good
approximation of the marginal posterior densities. A limitation of the
proposed models is their sensitivity to the assumption of the Poisson
distribution of the response, which fails when there is over-dispersion and
under-dispersion in the response observations. This issue is considered by
authors in their future research. The codes and data for this paper are
available online on github (soon).

## References

  * Ahmed (2014) Ahmed, S. E.(2014). Estimation Strategies in Poisson Regression Models. Penalty, Shrinkage, and Pretest Strategies: Variable Selection and Estimation, 101-115.
  * Algamal (2015) Algamal, Z. Y. (2015). Penalized Poisson regression model using adaptive modified elastic net penalty. _Electronic Journal of Applied Statistical Analysis_ , 8(2), 236-245.
  * Algamal and Lee (2015) Algamal, Z. Y., & Lee, M. H. (2015). Adjusted adaptive lasso in high-dimensional Poisson regression model. _Modern Applied Science_ , 9(4), 170.
  * Alonso et al. (2015) Alonso, J. B., Achcar, J. A., & Hotta, L. K. (2010). Climate changes and their effects in the public health: use of Poisson regression models. _Pesquisa Operacional_ , 30, 427-442.
  * Angers and Biswas (2003) Angers, J. F., & Biswas, A. (2003). A Bayesian analysis of zero-inflated generalized Poisson model. _Computational statistics & data analysis_, 42(1-2), 37-46.
  * Bai (2020) Bai, R. (2020). Spike-and-slab group lasso for consistent estimation and variable selection in non-Gaussian generalized additive models._arXiv preprint arXiv_ :2007.07021.
  * Bailey et al. (2009) Bailey, D. M., Collins, M. A., Gordon, J. D., Zuur, A. F., & Priede, I. G. (2009). Long-term changes in deep-water fish populations in the northeast Atlantic: a deeper reaching effect of fisheries? _Proceedings of the Royal Society B: Biological Sciences_ , 276(1664), 1965-1969.
  * Bar et al. (2020) Bar, H. Y., Booth, J. G., & Wells, M. T. (2020). A scalable empirical Bayes approach to variable selection in generalized linear models. _Journal of Computational and Graphical Statistics_ , 29(3), 535-546.
  * Benson and Friel (2021) Benson, A., & Friel, N. (2021). Bayesian inference, model selection and likelihood estimation using fast rejection sampling: the Conway-Maxwell-Poisson distribution. _Bayesian Analysis_ , 16(3), 905-931.
  * Blazére et al. (2014) Blazére, M., Loubes, J. M., & Gamboa, F. (2014). Oracle inequalities for a group lasso procedure applied to generalized linear models in high dimension. _IEEE Transactions on Information Theory_ , 60(4), 2303-2318.
  * Blei et al. (2017) Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. (2017). Variational inference: A review for statisticians. _Journal of the American statistical Association_ , 112(518), 859-877.
  * Bondell and Reich (2012) Bondell, H. D., & Reich, B. J. (2012). Consistent high-dimensional Bayesian variable selection via penalized credible regions. _Journal of the American Statistical Association_ , 107(500), 1610-1624.
  * Broderick et al. (2013) Broderick, T., Boyd, N., Wibisono, A., Wilson, A. C., & Jordan, M. I. (2013). Streaming variational bayes. _Advances in neural information processing systems_ , 26.
  * Carvalho et al. (2009) Carvalho, C. M., Polson, N. G., & Scott, J. G. (2009). Handling sparsity via the horseshoe. _Journal of Machine Learning Research_ W&CP, 5, 73-80.
  * Chan and Vasconcelos (2009) Chan, A. B., & Vasconcelos, N. (2009). Bayesian Poisson regression for crowd counting._In 2009 IEEE 12th international conference on computer vision_ 545-551
  * Chau et al. (2018) Chau, A. M. H., Lo, E. C. M., Wong, M. C. M., & Chu, C. H. (2018). Interpreting Poisson regression models in dental caries studies. _Caries Research_ , 52(4), 339-345.
  * Chen and Ibrahim (2000) Chen, M. H., & Ibrahim, J. G. (2000). Bayesian predictive inference for time series count data. _Biometrics_ , 56(3), 678-685.
  * Chen et al. (2010) Chen, Z., Kloosterman, F., Wilson, M. A., & Brown, E. N. (2010, March). Variational Bayesian inference for point process generalized linear models in neural spike trains analysis. _In 2010 IEEE International Conference on Acoustics, Speech and Signal Processing_ 2086-2089.
  * Christiansen and Morris (1997) Christiansen, C. L., & Morris, C. N. (1997). Hierarchical Poisson regression modeling. _Journal of the American Statistical Association_ , 92(438), 618-632.
  * Consul and Famoye (1992) Consul, P., & Famoye, F. (1992). Generalized Poisson regression model. _Communications in Statistics-Theory and Methods_ , 21(1), 89-109.
  * Coxe et al. (1992) Coxe, S., West, S. G., & Aiken, L. S. (2009). The analysis of count data: A gentle introduction to Poisson regression and its alternatives._Journal of personality assessment_ , 91(2), 121-136.
  * Crowther et al. (2012) Crowther, M. J., Riley, R. D., Staessen, J. A., Wang, J., Gueyffier, F., & Lambert, P. C. (2012). Individual patient data meta-analysis of survival data using Poisson regression models. _BMC Medical Research Methodology_ , 12(1), 1-14.
  * Dagne (2010) Dagne, G. A. (2010). Bayesian semiparametric zero-inflated Poisson model for longitudinal count data. _Mathematical biosciences_ , 224(2), 126-130.
  * D’Angelo et al. (2022) D’Angelo, L., & Canale, A. (2022). Efficient posterior sampling for Bayesian Poisson regression. _Journal of Computational and Graphical Statistics_ , 1-10.
  * Drugowitsch (2013) Drugowitsch, J. (2013). Variational Bayesian inference for linear and logistic regression. _arXiv preprint arXiv:1310.5438_.
  * Dvorzak and Wagner (2016) Dvorzak, M., & Wagner, H. (2016). Sparse Bayesian modelling of underreported count data. _Statistical Modelling_ , 16(1), 24-46.
  * El-Sayyad (1973) El-Sayyad, G. M. (1973). Bayesian and classical analysis of Poisson regression._Journal of the Royal Statistical Society: Series B (Methodological)_ , 35(3), 445-451.
  * Fair (1978) Fair, R. C. (1978). A theory of extramarital affairs. _Journal of political economy_ , 86(1), 45-61.
  * Fan and Li (2001) Fan, J., & Li, R. (2001). Variable selection via nonconcave penalized likelihood and its oracle properties. _Journal of the American statistical Association_ , 96(456), 1348-1360.
  * Fantin et al. (2023) Fantin, R., Barboza-Solís, C., Hildesheim, A., & Herrero, R. (2023). Excess mortality from COVID 19 in Costa Rica: a registry based study using Poisson regression. _The Lancet Regional Health–Americas_ , 20.
  * Famoye and Rothe (2003) Famoye, F., & Rothe, D. E. (2003). Variable selection for Poisson regression model. _Journal of Modern Applied Statistical Methods_ , 2(2), 11.
  * Gajda et al. (2010) Gajda, D., Guihenneuc-Jouyaux, C., Rousseau, J., Mengersen, K., & Nur, D. (2010). Use in practice of importance sampling for repeated MCMC for Poisson models._Electronic Journal of Statistics_ ,4 (2010) 361-383
  * Geman and Geman (1984) Geman, S. and D. Geman. (1984). Stochastic relaxation, Gibbs distribution and Bayesian restoration of images._IEE Transactions on Pattern Analysis and Machine Intelligence_ , 6: 721–741.
  * George and McCulloch (1993) George, E. I., & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. _Journal of the American Statistical Association,_ 88(423), 881-889.
  * George and McCulloch (1997) George, E. I., & McCulloch, R. E. (1997). Approaches for Bayesian variable selection. _Statistica Sinica_ , 7(2),339-373.
  * Geweke, J. (1996) Geweke, J.(1996). Variable selection and model comparison in regression._Bayesian statistics_ , 5, 609-620.
  * Goldsmith et al. (2011) Goldsmith, J., Wand, M. P., & Crainiceanu, C. (2011). Functional regression via variational Bayes. _Electronic journal of statistics_ , 5, 572.
  * Greene (2003) Greene, W. H. (2003). _Econometric analysis_. Pearson Education India.
  * Griffin and Brown (2013) Griffin, J. E., & Brown, P. J. (2013). Some priors for sparse regression modelling. _Bayesian Analysis_ 8(3), 691-702.
  * Guastavino and Benvenuto (2019) Guastavino, S., & Benvenuto, F. (2019). A consistent and numerically efficient variable selection method for sparse Poisson regression with applications to learning and signal recovery._Statistics and Computing_ , 29(3), 501-516.
  * Hastings (1970) Hastings, W. K. (1970). Monte Carlo sampling methods using Markov chains and their applications._Biometrika_ ,97-109
  * Hayat and Higginsm (2014) Hayat, M. J., & Higgins, M. (2014). Understanding Poisson regression. _Journal of Nursing Education_ , 53(4), 207-215.
  * Hilbe (2014) Hilbe, J. M. (2014). _Modeling count data_. Cambridge University Press.
  * Hossain and Ahmed (2012) Hossain, S., & Ahmed, E. (2012). Shrinkage and penalty estimators of a Poisson regression model. _Australian & New Zealand Journal of Statistics_, 54(3), 359-373.
  * Huang et al. (2016) Huang, X., Wang, J., & Liang, F. (2016). A variational algorithm for Bayesian variable selection. _arXiv preprint arXiv_ , 1602.07640.
  * Igarashi et al. (2018) Igarashi, Y., Takenaka, H., Nakanishi-Ohno, Y., Uemura, M., Ikeda, S., & Okada, M. (2018). Exhaustive search for sparse variable selection in linear regression. _Journal of the Physical Society of Japan_ , 87(4), 044802.
  * Ishwaran and Rao (2003) Ishwaran, H., & Rao, J. S. (2003). Detecting differentially expressed genes in microarrays using Bayesian model selection. _Journal of the American Statistical Association_ , 98(462), 438-455.
  * Ishwaran and Rao (2005) Ishwaran, H., & Rao, J. S. (2005). Spike and slab variable selection: frequentist and Bayesian strategies. _Annals of Statistics_ , 33(2): 730-773.
  * Ivanoff et al. (2016) Ivanoff, S., Picard, F., & Rivoirard, V. (2016). Adaptive Lasso and group-Lasso for functional Poisson regression. _The Journal of Machine Learning Research_ , 17(1), 1903-1948.
  * Jaakkola and Jordan (2000) Jaakkola, T. S., & Jordan, M. I. (2000). Bayesian parameter estimation via variational methods. _Statistics and Computing_ , 10, 25-37.
  * Jeong and Ghosal (2021) Jeong, S., & Ghosal, S. (2021). Posterior contraction in sparse generalized linear models. _Biometrika_ , 108(2), 367-379.
  * Ji and Schmidler (2013) Ji, C., & Schmidler, S. C. (2013). Adaptive Markov chain Monte Carlo for Bayesian variable selection. _Journal of Computational and Graphical Statistics_ , 22(3), 708-728.
  * Jia et al. (2013) Jia, J., Rohe, K., & Yu, B. (2013). The lasso under Poisson-like heteroscedasticity. Statist. _Sinica_ , 23(1):99-118.
  * Jia et al. (2019) Jia, J., Xie, F., & Xu, L. (2019). Sparse Poisson regression with penalized weighted score function. _Electronic Journal of Statistics_ , 13 (2) 2898-2920.
  * Jiang (2007) Jiang, W. (2007). Bayesian variable selection for high dimensional generalized linear models: convergence rates of the fitted densities. _Ann. Statist_. 35(4): 1487-1511.
  * Khafrim et al. (2008) Khafrim, S., Kazemnejad, A., & Eskandari, F. (2008). Hierarchical Bayesian analysis of bivariate Poisson regression model._World Applied Sciences Journal_ , 4, 667-675.
  * Kim et al. (2013) Kim, S., Chen, Z., Zhang, Z., Simons-Morton, B. G., & Albert, P. S. (2013). Bayesian hierarchical Poisson regression models: an application to a driving study with kinematic events. _Journal of the American Statistical Association_ , 108(502), 494-503.
  * Koç et al. (2013) Koç, H., Cengiz, M. A., & Koç, T. (2013). Convergence diagnostics in MCMC for zero-inflated Poisson models for air pollution data. _Journal of Advanced Scientific Research_ , 4(04), 12-15.
  * Li and Cevher. (2015) Y. H. Li & V. Cevher. (2015). Consistency of $\ell_{1}$-regularized maximum-likelihood for compressive Poisson regression. _International Conference on Acoustics, Speech, and Signal Processing, South Brisbane, QLD, Australia_ , 2015, pp. 3606-3610
  * Li and Lin (2010) Li, Q., & Lin, N. (2010). The Bayesian elastic net. _Bayesian Analysis_ ,151-170
  * Luts and Wand (2015) Luts, J., & Wand, M. P. (2015). Variational inference for count response semiparametric regression. _Bayesian analysis_ , 991-1023.
  * Legisso et al. (2023) Legisso, T. Z., Mamo, B. G., Bimrew, A. M., & Fikadu, T. (2023). Blood Pressure Examination Habit and Its Determinants Among Civil Servants in Arba Minch Town: A Cross-Sectional Study–Using Hurdle Poisson Regression Model. _Integrated Blood Pressure Control_ , 1-9.
  * Minka (2000) Minka, T. P. (2000). _Variational Bayes for mixture models: Reversing em_. Technical report, MIT.
  * Minka (2001) Minka, T. P. (2001). _A family of algorithms for approximate Bayesian inference._ Doctoral dissertation, Massachusetts Institute of Technology.
  * Minka (2013) Minka, T. P. (2013). Expectation propagation for approximate Bayesian inference. _arXiv preprint arXiv_ :1301.2294.
  * Mitchell and Beauchamp (1988) Mitchell, T.,& Beauchamp, J. (1988). Bayesian variable selection in linear regression. _Journal of the American Statistical Association_ , 83(404), 1023-1032.
  * Mohammadi et al. (2016) Mohammadi, T., Kheiri, S., & Sedehi, M. (2016). Analysis of blood transfusion data using bivariate zero-inflated poisson model: a Bayesian approach. _Computational and Mathematical Methods in Medicine_ , 7878325
  * Mzoughia et al. (2018) Mzoughia, M. B., Borle, S., & Limam, M. (2018). A MCMC approach for modeling customer lifetime behavior using the COM-Poisson distribution._Applied Stochastic Models in Business and Industry_ , 34(2), 113-127.
  * Nishimura and Suchard (2022) Nishimura, A., & Suchard, M. A. (2022). Prior-preconditioned conjugate gradient method for accelerated Gibbs sampling in “large n, large p” Bayesian sparse regression. _Journal of the American Statistical Association_ , 1-14.
  * Noori Asl et al. (2022) Noori Asl, M., Bevrani, H., & Arabi Belaghi, R. (2022). Penalized and ridge-type shrinkage estimators in Poisson regression model. _Communications in Statistics-Simulation and Computation_ , 51(7), 4039-4056.
  * O’hara and Sillanpää (2009) O’hara, R. B., & Sillanpää, M. J. (2009). A review of Bayesian variable selection methods: what, how and which. _Bayesian Analysis_. 4(1): 85-117.
  * Ormerod et al. (2017) Ormerod, J. T., You, C., & Müller, S. (2017). A variational Bayes approach to variable selection _Electronic Journal of Statistics_ , 11 (2), 3549-3594.
  * Peng (2023) Peng, L. (2023). Oracle inequalities for weighted group Lasso in high-dimensional Poisson regression model. _Communications in Statistics-Theory and Methods_ , 1-27.
  * Picard and Cook (1984) Picard, R. R., & Cook, R. D. (1984). Cross-validation of regression models. _Journal of the American Statistical Association_ , 79(387), 575-583.
  * Pham et al. (2013) Pham, T. H., Ormerod, J. T., & Wand, M. P. (2013). Mean field variational Bayesian inference for nonparametric regression with measurement error. _Computational Statistics & Data Analysis_, 68, 375-387.
  * Polson et al. (2014) Polson, N. G., Scott, J. G., & Windle, J. (2014). The Bayesian Bridge. _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 76(4), 713-733.
  * Powers et al. (2010) Powers, S., Gerlach, R., & Stamey, J. (2010). Bayesian variable selection for Poisson regression with underreported responses. _Computational statistics & data analysis_, 54(12), 3289-3299.
  * Raman and Roth (2009) Raman, S., & Roth, V. (2009). Sparse Bayesian regression for grouped variables in generalized linear models. In _Joint Pattern Recognition Symposium_. Berlin, Heidelberg: Springer Berlin Heidelberg, 242-251
  * Ray and Szabó (2022) Ray, K., & Szabó, B. (2022). Variational Bayes for high-dimensional linear regression with sparse priors. _Journal of the American Statistical Association_ , 117(539), 1270-1281.
  * Ročková (2018) Ročková, V. (2018). Bayesian estimation of sparse signals with a continuous spike-and-slab prior. _Annals of Statistics,_ 46(1): 401-437
  * Ročková and George (2014) Ročková, V., & George, E. I. (2014). EMVS: The EM approach to Bayesian variable selection. _Journal of the American Statistical Association_ , 109(506), 828-846.
  * Richardson et al. (2010) Richardson, S., Bottolo, L., & Rosenthal, J. S. (2010). Bayesian models for sparse regression analysis of high dimensional data. _Bayesian statistics_ , 9, 539-569.
  * Rue et al. (2010) Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 71(2), 319-392.
  * Saishu et al. (2021) Saishu, H., Kudo, K., & Takano, Y. (2021). Sparse Poisson regression via mixed-integer optimization. _Plos one_ , 16(4), e0249916.
  * Seeger (2008) Seeger, M. W. (2008). Bayesian inference and optimal design for the sparse linear model. _The Journal of Machine Learning Research_ , 9, 759-813.
  * Sellers et al. (2012) Sellers, K. F., Borle, S., & Shmueli, G. (2012). The COM-Poisson model for count data: a survey of methods and applications. _Applied Stochastic Models in Business and Industry_ , 28(2), 104-116.
  * Sellers and Premeaux (2021) Sellers, K. F., & Premeaux, B. (2021). Conway-Maxwell-Poisson regression models for dispersed count data. _Wiley Interdisciplinary Reviews: Computational Statistics_ , 13(6), e1533
  * Shankar Raman (2012) Shankar Raman, S. (2012). _Bayesian grouped variable selection_. Doctoral dissertation, University of Basel.
  * Shao (1993) Shao, J. (1993). Linear model selection by cross-validation. Journal of the American Statistical Association, 88(422), 486-494.
  * Shin and Liu (2022) Shin, M., & Liu, J. S. (2022). Neuronized priors for Bayesian sparse linear regression. _Journal of the American Statistical Association_ , 117(540), 1695-1710.
  * Šmídl and Quinn (2006) Šmídl, V., & Quinn, A. (2006). _The variational Bayes method in signal processing. Springer Science & Business Media._
  * Solomon et al. (2022) Solomon, E. T., Gari, S. R., & Alemu, B. M. (2022). Prevalence and risk factors of diarrhea among children less than five years of age in the rural suburbs of Dire Dawa, Eastern Ethiopia; Robust Poisson Regression Analysis. _African Health Sciences_ , 4(4), 653-663.
  * Tang and Martin (2023) Tang, Y., & Martin, R. (2023). Empirical Bayes inference in sparse high-dimensional generalized linear models. _arXiv preprint arXiv_ :2303.07854.
  * Tibshirani (1996) Tibshirani, R. (1996). Regression shrinkage and selection via the lasso. _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 58(1), 267-288.
  * Tipping (2001) Tipping, M. E. (2001). Sparse Bayesian learning and the relevance vector machine. _The Journal of Machine Learning Research_ , 1, 211-244.
  * Titsias and Lázaro-Gredilla (2011) Titsias, M., & Lázaro-Gredilla, M. (2011). Spike and slab variational inference for multi-task and multiple kernel learning. _dvances in neural information processing systems_ , 24.
  * Titsias and Lázaro-Gredilla (2011) Titsias, M., & Lázaro-Gredilla, M. (2014). Doubly stochastic variational Bayes for non-conjugate inference. _International conference on machine learning_ ,1971-1979
  * Tomal et al. (2022) Tomal, J. H., Khan, J. R., & Wahed, A. S. (2022). Weighted Bayesian Poisson Regression for The Number of Children Ever Born per Woman in Bangladesh. _Journal of Statistical Theory and Applications_ , 21(3), 79-105.
  * Tran et al. (2023) Tran, M. N., Tseng, P., & Kohn, R. (2023). Particle Mean Field Variational Bayes. _arXiv preprint arXiv_ ,2303.13930.
  * Viallefont et al. (2002) Viallefont, V., Richardson, S., & Green, P. J. (2002). Bayesian analysis of Poisson mixtures._Journal of nonparametric statistics_ , 14(1-2), 181-202.
  * Wand et al. (2011) Wand, M. P., Ormerod, J. T., Padoan, S. A., & Frühwirth, R. (2011). Mean field variational Bayes for elaborate distributions._Bayesian Analysis_ , 847-900
  * Winn et al. (2005) Winn, J., Bishop, C. M., & Jaakkola, T. (2005). Variational message passing. _Journal of Machine Learning Research_ , 6(4).
  * Winter and Bürkner (2021) Winter, B., & Bürkner, P. C. (2021). Poisson regression for linguists: A tutorial introduction to modelling count data with brms. _Language and Linguistics Compass_ , 15(11), e12439.
  * Yousif and Gatea (2021) Yousif, A. H., & Gatea, A. H. (2021). Atan regularized for the high dimensional Poisson regression model. _International Journal of Nonlinear Analysis and Applications_ , 12(Special Issue), 2197-2202.
  * Xie et al. (2014) Xie, F. C., Lin, J. G., & Wei, B. C. (2014). Bayesian zero-inflated generalized Poisson regression model: estimation and case influence diagnostics. _Journal of Applied Statistics_ , 41(6), 1383-1392.
  * Xu and Ghosh (2015) Xu, X., & Ghosh, M. (2015). Bayesian variable selection and estimation for group lasso. _Bayesian Analysis_ , 10(4): 909-936
  * Zhang et al. (2019) Zhang, C. X., Xu, S., & Zhang, J. S. (2019). A novel variational Bayesian method for variable selection in logistic regression models. _Computational Statistics & Data Analysis_, 133, 1-19.
  * Zhang (2010) Zhang, C. H. (2010). Nearly unbiased variable selection under minimax concave penalty. _Ann. Statist_ , 38(2): 894-942
  * Zhang and Ye (2014) Zhang, Y., & Ye, K. (2014). Bayesian D-optimal designs for Poisson regression models. _Communications in Statistics-Theory and Methods_ , 43(6), 1234-1247.
  * Zou and Li (2008) Zou, H., & Li, R. (2008). One-step sparse estimates in nonconcave penalized likelihood models. Annals of Statistics, 36(4), 1509.
  * Zuur et al. (2013) Zuur, A. F., Hilbe, J. M., & Leno, E. N. (2013). _A beginner’s guide to GLM and GLMM with R: A frequentist and Bayesian perspective for ecologists_ , 270 p Newburgh. UK: Highland Statistics Ltd.

## Appendix (A): Derivation of $q^{*}$ density functions

Here, we propose the details of the derivation of the VB components, as well
as, the computation of the ELBO, for the three proposed VB models.

### Laplace prior

For model (3), and using (2), we have

$\displaystyle\log q(\beta)=$ $\displaystyle
E_{-\beta}\left[\log\tilde{P}(y|\beta,\xi)+\log
P(\beta|\tau)\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle
E_{-\beta}\left[\sum_{i=1}^{n}M_{\xi_{i}}X_{i}\beta+\sum_{i=1}^{n}y_{i}X_{i}\beta\right.$
$\displaystyle-\left.\frac{1}{2}\sum_{i=1}^{n}e^{\xi_{i}}\beta^{T}X_{i}X_{i}^{T}\beta-\frac{1}{2}\beta^{T}\text{diag}(\tau_{j}^{-1})\beta\right]+\text{Const.}$
$\displaystyle=$
$\displaystyle(y-M_{\xi})^{T}X\beta-\frac{1}{2}\beta^{T}\left[\sum_{i=1}^{n}e^{\xi_{i}}X_{i}X_{i}^{T}+\text{diag}(E(\tau^{-1}))\right]\beta+\text{Const.,}$

which is the kernel of the Gaussian distribution. Also, we have

$\displaystyle\log q(\lambda)=$ $\displaystyle
E_{-\lambda}\left[\sum_{j=1}^{p}\log P(\tau_{j}|\lambda)+\log
P(\lambda)\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle
E_{-\lambda}\left[P\log\lambda-\frac{\lambda}{2}\sum_{j=1}^{p}\tau_{j}+(\nu-1)\log\lambda-\lambda\delta\right]+\text{Const.}$
$\displaystyle=$
$\displaystyle(p+\nu-1)\log\lambda-\lambda(\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})+\delta)+\text{Const.},$

which is the kernel of the gamma distribution. For the variance parameters
$\tau_{j}$, for $j=1,\ldots,p$, we can see that

$\displaystyle\log q(\tau_{j})=$ $\displaystyle\left[\log
P(\beta_{j}|\tau_{j})+\log P(\tau_{j})\right]+\text{Const.}$ $\displaystyle=$
$\displaystyle
E\left[-\frac{1}{2}\log(\tau_{j})-\frac{1}{2}\dfrac{\beta_{j}^{2}}{\tau_{j}}-\frac{\lambda}{2}\tau_{j}\right]+\text{Const.}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\log(\tau_{j})-\frac{1}{2}\dfrac{E[\beta_{j}^{2}]}{\tau_{j}}-\frac{E[\lambda]}{2}\tau_{j}+\text{Const.},$

which is the logarithm of the generalized inverse Gaussian density, and for
$\tau_{0}$, we have

$\displaystyle\log q(\tau_{0})=$ $\displaystyle E_{-\tau_{0}}\left[\log
P(\beta_{0}|\tau_{0})+\log P(\tau_{0}|a)\right]+\text{Const.}$
$\displaystyle=$ $\displaystyle
E(-\frac{1}{2}\log\tau_{0}-\frac{1}{2}\dfrac{\beta_{0}^{2}}{\tau_{0}}-\frac{3}{2}\log\tau_{0}-\dfrac{1}{a\tau_{0}})+\text{Const.}$
$\displaystyle=$
$\displaystyle-2\log\tau_{0}-\frac{1}{\tau_{0}}\left[\frac{1}{2}E(\beta_{0}^{2})+E(a^{-1})\right]+\text{Const.},$

Finally, we can see that

$\displaystyle\log q(a)=$ $\displaystyle E_{-a}\left[\log P(\tau_{0}|a)+\log
P(a)\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle
E_{-a}\left[-\dfrac{1}{2}\log a-\dfrac{1}{a\tau_{0}}+\dfrac{3}{2}\log
a-\dfrac{1}{Aa}\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle-2\log
a-\dfrac{1}{a}(E(\tau_{0}^{-1})+\dfrac{1}{A})+\text{Const.},$

which is the kernel of the inverse gamma distribution.

To obtain the ELBO, we compute the following terms

$\displaystyle E_{q}(\log\tilde{P}(y|\beta,\xi))=$ $\displaystyle-
M_{\xi}^{T}(1+X_{i}\mu_{\beta})-\sum_{i=1}{n}\frac{\xi_{i}^{2}}{2}e^{\xi_{i}}$
$\displaystyle-\frac{1}{2}{\rm tr}(S_{X}^{\xi}D_{\beta})+y^{T}X\mu_{\beta},$
$\displaystyle E(\log P(\beta|\tau))=$
$\displaystyle\frac{1}{2}E(\log\tau_{0})-\frac{1}{2}\sum_{j=1}^{p}E(\log\tau_{j})$
$\displaystyle-\frac{1}{2}E(\tau_{0}^{-1})E(\beta_{0}^{2})-\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}^{-1})E(\beta_{j}^{2}),$
$\displaystyle E(\log P(\lambda))=$ $\displaystyle(\nu-1)E(\log\lambda)-\delta
E(\lambda)$ $\displaystyle=$
$\displaystyle(\nu-1)\left[\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))\right]$
$\displaystyle-$
$\displaystyle\delta\dfrac{p+\nu}{\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})},$
$\displaystyle E(\log P(\tau_{0}|a))=$ $\displaystyle E(-\frac{1}{2}\log
a-\frac{3}{2}\log\lambda_{0}-\frac{1}{a\tau_{0}})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\log(\frac{1}{A}+E(\tau_{0}^{-1}))-\frac{3}{2}\log(\frac{1}{2}d_{00}+E(a^{-1}))$
$\displaystyle-E(a^{-1})E(\tau_{0}^{-1}),$ $\displaystyle E(\log
P(\tau_{j}|\lambda))=$ $\displaystyle
E\left[\log\lambda-\frac{1}{2}\lambda\tau_{j}\right]$ $\displaystyle=$
$\displaystyle\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))-\frac{1}{2}E(\lambda)E(\tau_{j}),$
$\displaystyle E_{q}(\log P(a))=$ $\displaystyle E_{q}(-\frac{3}{2}\log
a-\dfrac{1}{Aa})$ $\displaystyle=$
$\displaystyle\frac{3}{2}\Big{(}\psi(1)-\log\big{(}\frac{1}{A}+E(\tau_{0}^{-1})\big{)}\Big{)}-\dfrac{\frac{1}{A}}{\frac{1}{A}+E(\tau_{0}^{-1})},$
$\displaystyle-E(\log q(\beta))=+\frac{1}{2}\log|\Sigma_{\beta}|,$
$\displaystyle-E(\log q(\lambda))=$
$\displaystyle-(p+\nu-1)\left[\psi(p+\nu)-\log(\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j}))\right]$
$\displaystyle-(p+\nu)\log\big{(}\delta+\frac{1}{2}\sum_{j=1}^{p}E(\tau_{j})\big{)},$
$\displaystyle-E(\log q(\tau_{0}))=$ $\displaystyle
2E(\log(\tau_{0}))+E(\tau_{0}^{-1})\big{(}\frac{1}{2}d_{00}+E(a^{-1})\big{)}$
$\displaystyle-\log(\frac{1}{2}d_{00}+E(a^{-1}))$ $\displaystyle=$
$\displaystyle\log(\frac{1}{2}d_{00}+E(a^{-1})),$ $\displaystyle-E(\log
q(\tau_{j}))=$
$\displaystyle-\frac{1}{2}E(\log\tau_{j})+\frac{1}{2}E(\tau_{j}^{-1})d_{jj}$
$\displaystyle+\frac{1}{2}E(\lambda)E(\tau_{j})-\frac{1}{2}\log(\lambda_{j}^{\ast}),$

and

$\displaystyle-E(\log q(a))=$
$\displaystyle-\log\big{(}\frac{1}{A}+E(\tau_{j}^{-1})\big{)}+2E(\log
a)-\log\big{(}\frac{1}{A}+E(\tau_{0})\big{)}$ $\displaystyle=$ $\displaystyle
2E(\log a)+E(a^{-1})\big{(}E(\tau_{0}^{-1})+\frac{1}{A}\big{)}$
$\displaystyle=$ $\displaystyle
2\log\big{(}E(\tau_{0}^{-1})+\frac{1}{A}\big{)}-\psi(1)+E(a^{-1})\big{(}E(\tau_{0}^{-1})+\frac{1}{A}\big{)}$
$\displaystyle=$
$\displaystyle\log\big{(}E(\tau_{0}^{-1})+\frac{1}{A}\big{)}.$

Adding the above terms, we reach in (7).

### Continuous spike and slab prior

For model (10), we can see that

$\displaystyle\log q(\beta)=$
$\displaystyle(y-M_{\xi})^{T}X\beta-\frac{1}{2}\beta^{T}S_{X}^{\xi}\beta-\frac{1}{2}E(\sigma^{-2})\beta^{T}\text{diag}(P)\beta$
$\displaystyle-\frac{1}{2}c^{-1}E(\sigma^{-2})\beta^{T}(\text{I}_{p}-\text{diag}(P))\beta+\text{Const.},$

which is again the kernel of the Gaussian distribution,

$\displaystyle\log q(\sigma^{2})=$ $\displaystyle E_{-\sigma^{2}}\left[\log
P(\beta|\sigma^{2})+\log P(\sigma^{2}|a)\right]+\text{Const.}$
$\displaystyle=$
$\displaystyle-\frac{p}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{2}}E\left[\text{tr(diag}(Z))\beta\beta^{T}\right]$
$\displaystyle-\dfrac{1}{2c\sigma^{2}}E\left[\text{tr(diag}(1-Z))\beta\beta^{T}\right]-\frac{1}{\sigma^{2}}E(a^{-1})$
$\displaystyle-\dfrac{1}{2}\log\sigma^{2}+\text{Const.}$ $\displaystyle=$
$\displaystyle-\dfrac{p+1}{2}\log\sigma^{2}-\dfrac{1}{\sigma^{2}}\left[\frac{1}{2}\text{tr(diag}(P)D_{\beta})\right.$
$\displaystyle+\left.\dfrac{1}{2c}\text{tr(diag}(1-P)D_{\beta})+E(a^{-1})\right]+\text{Const.},$

which is the logarithm of an inverse gamma density, and

$\displaystyle\log q(a)=$ $\displaystyle E_{-a}\left[\log P(\sigma^{2}|a)+\log
P(a)\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle
E_{-a}\left[\dfrac{1}{2}\log a-\dfrac{1}{a\sigma^{2}}+\dfrac{3}{2}\log
a-\dfrac{1}{Aa}\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle 2\log
a-\dfrac{1}{a}(E(\sigma^{-2})+\dfrac{1}{A})+\text{Const.},$

which corresponds with an inverse gamma distribution. Furthermore, for
$j=1,\ldots,p$,

$\displaystyle\log q(Z_{j})=$ $\displaystyle E_{-Z_{j}}\left[\log
P(\beta|\sigma^{2})+\log P(Z_{j})\right]+\text{Const.}$ $\displaystyle=$
$\displaystyle
E_{-Z_{j}}\Big{[}Z_{j}(-\dfrac{1}{2\sigma^{2}}\beta_{j}^{2})+(1-Z_{j})(-\dfrac{1}{2c\sigma^{2}}\beta_{j}^{2})$
$\displaystyle+Z_{j}\log\pi_{j}+(1-Z_{j})\log(1-\pi_{j})\Big{]}+\text{Const.},$

which is the kernel of a Bernoulli probability mass function,

$\displaystyle\log q(\pi_{j})=$ $\displaystyle E_{-\pi_{j}}\left[\log
P(Z_{j}|\pi_{j})+\log P(\pi_{j})\right]+\text{Const.}$ $\displaystyle=$
$\displaystyle
E\Big{[}Z_{j}\log\pi_{j}+(1-Z_{j})(1-\pi_{j})+(\theta_{1}-1)\log\pi_{j}$
$\displaystyle+(\theta_{2}-1)\log(1-\pi_{j})\Big{]}+\text{Const.}$
$\displaystyle=$
$\displaystyle(P_{j}+\theta_{1}-1)\log\pi_{j}+(\theta_{2}-P_{j})\log(1-\pi_{j})+\text{Const.},$

which is the beta distribution.

Now we compute the following terms for obtaining ELBO

$\displaystyle E_{q}(\log\tilde{P}(y|\beta,\xi))=$ $\displaystyle-
M_{\xi}X\mu_{\beta}-\sum_{i=1}^{n}\left(M_{\xi_{i}}-e^{\xi_{i}}\frac{\xi_{i}}{2}^{2}\right)$
$\displaystyle-\frac{1}{2}{\rm tr}(S_{X}^{\xi}D_{\beta})+y^{T}X\mu_{\beta},$
$\displaystyle E_{q}(\log P(\beta))=$ $\displaystyle
E_{q}\Big{[}\sum_{j=1}^{p}Z_{j}\big{(}-\frac{p}{2}\log
2\pi-\frac{1}{2}\log\sigma^{2}-\dfrac{1}{2\sigma^{2}}\beta_{j}^{2}\big{)}$
$\displaystyle-(1-Z_{j})\big{(}-\frac{p}{2}\log 2\pi-\frac{1}{2}\log
c\sigma^{2}-\dfrac{1}{2c\sigma^{2}}\beta_{j}^{2}\big{)}\Big{]}$
$\displaystyle=$ $\displaystyle\sum_{j=1}^{p}P_{j}\big{(}-\frac{p}{2}\log
2\pi-\frac{1}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}})-\frac{1}{2}E(\sigma^{-2})d_{jj}\big{)}$
$\displaystyle+(1-P_{j})\big{(}-\frac{p}{2}\log
2\pi-\frac{1}{2c}E(\sigma^{-2})d_{jj}-\frac{1}{2}\log c$
$\displaystyle-\frac{1}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}})\big{)}$
$\displaystyle=$ $\displaystyle-\frac{p}{2}\log
2\pi-\frac{1}{2}\log(1-\sum_{j=1}^{p}P_{j})-\frac{1}{2}(1-\frac{1}{c})E(\sigma^{-2})\sum_{j=1}^{p}P_{j}d_{jj}$
$\displaystyle-\frac{1}{2c}E(\sigma^{-2})\sum_{j=1}^{p}d_{jj}-\frac{p}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}}),$
$\displaystyle E(\log P(\sigma^{2}))=$ $\displaystyle E\Big{[}-\frac{1}{2}\log
a-\frac{3}{2}E(\log\sigma^{2})-\frac{1}{2}E(\sigma^{-2})-\log\Gamma(\frac{1}{2})\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{2}(\psi(1)-\log(\frac{1}{A}+E(\sigma^{-2})))-\log\Gamma(\frac{1}{2})$
$\displaystyle-\frac{3}{2}(\psi(\alpha_{\sigma^{2}})-\log\beta_{\sigma^{2}})-\dfrac{E(\sigma^{-2})}{\frac{1}{A}+E(\sigma^{-2})},$
$\displaystyle E(\log P(Z))=$
$\displaystyle\sum_{j=1}^{p}E\left[Z_{j}\log\pi_{j}+(1-Z_{j})\log(1-\pi_{j})\right]$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{p}P_{j}(\psi(\alpha^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))$
$\displaystyle+\sum_{j=1}^{p}(1-P_{j})(\psi(\beta^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j})),$
$\displaystyle E(\log P(\pi))=$
$\displaystyle\sum_{j=1}^{p}\log\dfrac{\Gamma(\theta_{1})\Gamma(\theta_{2})}{\Gamma(\theta_{1}+\theta_{2})}+(\theta_{1}-1)\sum_{j=1}^{p}(\psi(\alpha^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j}))$
$\displaystyle+(\theta_{2}-1)\sum_{j=1}^{p}(\psi(\beta^{\ast}_{j})-\psi(\beta^{\ast}_{j}+\alpha^{\ast}_{j})),$
$\displaystyle E_{q}(\log P(a))=$ $\displaystyle E_{q}(-\frac{3}{2}\log
a-\dfrac{1}{Aa})$ $\displaystyle=$
$\displaystyle\frac{3}{2}\Big{(}\psi(1)-\log\big{(}\frac{1}{A}+E(\sigma^{-2})\big{)}\Big{)}-\dfrac{\frac{1}{A}}{\frac{1}{A}+E(\sigma^{-2})},$
$\displaystyle-E(\log q(\beta))=+\frac{1}{2}\log|\Sigma_{\beta}|,$
$\displaystyle-E(\log q(\sigma^{2}))=$
$\displaystyle-\alpha_{\sigma^{2}}\log\beta_{\sigma^{2}}+\log\Gamma(\alpha_{\sigma^{2}})$
$\displaystyle+(\alpha_{\sigma^{2}}-1)E(\log\sigma^{2})+\beta_{\sigma^{2}}E(\sigma^{-2})$
$\displaystyle=$
$\displaystyle-\alpha_{\sigma^{2}}\log\beta_{\sigma^{2}}+\log\Gamma(\alpha_{\sigma^{2}})+(\alpha_{\sigma^{2}}-1)$
$\displaystyle\left[\psi(\alpha_{\sigma^{2}})-\beta_{\sigma^{2}}\right]+\alpha_{\sigma^{2}},$
$\displaystyle-E(\log q(Z))=$
$\displaystyle-E\left[\sum_{j=1}^{p}Z_{j}\log(P_{j})+(1-Z_{j})\log(1-P_{j})\right]$
$\displaystyle=$
$\displaystyle-E\left[\sum_{j=1}^{p}P_{j}\log(P_{j})+(1-P_{j})\log(1-P_{j})\right],$

and

$\displaystyle-E(\log q(a))=$
$\displaystyle-\log\Big{(}\frac{1}{A}+E(\sigma^{-2})\Big{)}+2E(\log
a)+\Big{(}\frac{1}{A}+E(\sigma^{-2})\Big{)}E(\frac{1}{a})$ $\displaystyle=$
$\displaystyle-\log\Big{(}\frac{1}{A}+E(\sigma^{-2})\Big{)}-2\big{(}\psi(1)-\log\left(\frac{1}{A}+E(\sigma^{-2})\right)$
$\displaystyle=$ $\displaystyle\log\Big{(}\frac{1}{A}+E(\sigma^{-2})\Big{)}.$

Summation of the above terms would result in (12).

### Bernoulli sparsity enforcing prior

For model (13), the VB elements are computed as follows.

$\displaystyle\log q(\beta)=$ $\displaystyle
E_{-\beta}\left[\log\tilde{P}(y;\beta,\gamma,\xi)+\log
P(\beta|A)\right]+\text{Const.}$ $\displaystyle=$ $\displaystyle
E_{-\beta}\left[(y-M_{\xi})^{T}X\Gamma\beta-\frac{1}{2}\tilde{P}(y;\beta,\gamma,\xi)-\frac{1}{2}\beta^{T}A\beta\right]+\text{Const.}$
$\displaystyle=$
$\displaystyle(y-M_{\xi})^{T}X\text{diag}(P)\beta-\frac{1}{2}\beta^{T}\left[S_{X}^{\xi}\odot\Omega+E(A^{-1})\right]\beta+\text{Const.},$

which results in a Gaussian element. Also, for $j=1,\ldots,p$, we have

$\displaystyle\log q(\gamma_{j})=$ $\displaystyle
E_{-\gamma_{j}}\left[\log\tilde{P}(y;\beta,\gamma,\xi)+\log(\gamma_{j}|\pi_{j})\right]+\text{Const.}$
$\displaystyle=$ $\displaystyle
E_{-\gamma_{j}}\Big{[}(y-M_{\xi})^{T}X\Gamma\beta-\frac{1}{2}\beta^{T}\Gamma
S_{X}^{\xi}\Gamma\beta$
$\displaystyle+\gamma_{j}\log\big{(}\dfrac{\pi_{j}}{1-\pi_{j}}\big{)}\Big{]}+\text{Const.}$
$\displaystyle=$ $\displaystyle
E_{-\gamma_{j}}\Big{[}(y-M_{\xi})^{T}X_{j}\gamma_{j}\beta_{j}-\frac{1}{2}\gamma_{j}S_{X_{jj}}^{\xi}\beta_{j}^{2}$
$\displaystyle-\frac{1}{2}\sum_{i\neq
j}\gamma_{i}\gamma_{j}S_{X_{ij}}^{\xi}\beta_{i}\beta_{j}+\gamma_{j}\big{(}\dfrac{\pi_{j}}{1-\pi_{j}}\big{)}\Big{]}+\text{Const.}$
$\displaystyle=$
$\displaystyle\gamma_{j}\Big{[}(y-M_{\xi})^{T}X_{j}\beta_{j}-\frac{1}{2}S_{X_{jj}}^{\xi}d_{jj}$
$\displaystyle-\frac{1}{2}\sum_{i\neq
j}P_{i}S_{X_{ij}}^{\xi}d_{ij}+\psi(\alpha_{j}^{\ast})-\psi(\beta_{j}^{\ast})\Big{]}+\text{Const.},$

which is the kernel of a Bernoulli distribution. Furthermore, for
$j=0,1,\ldots,p$,

$\displaystyle\log q(\alpha_{j})=$ $\displaystyle E_{-\alpha_{j}}\left[\log
P(\beta|A^{-1})+\log P(\alpha_{j})\right]+\text{Const.}$ $\displaystyle=$
$\displaystyle
a_{j}\log\alpha_{j}-(b_{j}+\frac{d_{jj}}{2})\alpha_{j}+\text{Const.},$

which corresponds with a gamma density.

$\displaystyle\log P(\pi_{j})=$ $\displaystyle E_{-\pi_{j}}\left[\log
P(\gamma_{j}|\pi_{j})+\log P(\pi_{j})\right]+\text{Const.}$ $\displaystyle=$
$\displaystyle
P(P_{j}+\theta_{1}-1)\log\pi_{j}+(\theta_{2}-P_{j})\log(1-\pi_{j})+\text{Const.},$

which is the logarithm of a beta density.

The ELBO in (16) is computed by computation of

$\displaystyle E(\log\tilde{P}(y|\beta,\gamma,\xi))=$ $\displaystyle
E_{q}\Big{[}M_{\xi}^{T}X\Gamma\beta-\frac{1}{2}\beta^{T}\Gamma
S_{X}^{\xi}\Gamma\beta$
$\displaystyle+y^{T}X\Gamma\beta-\sum_{i=1}^{n}(e^{\xi_{i}}+\xi_{i}-\frac{\xi_{i}^{2}}{2})\Big{]}$
$\displaystyle=$ $\displaystyle-
M_{\xi}XP^{\ast}\mu_{\beta}-\frac{1}{2}\text{tr}(\Sigma_{\beta}+\mu_{\beta}\mu_{\beta}^{T})(S_{X}^{\xi}\odot\Omega)$
$\displaystyle+\mu_{\beta}^{T}\Sigma_{\beta}^{-1}\mu_{\beta}-\sum_{i=1}^{n}\big{(}e^{\xi_{i}}+\xi_{i}-\frac{\xi_{i}^{2}}{2}\big{)},$
$\displaystyle E(\log P(\beta|\alpha))=$
$\displaystyle\frac{1}{2}E_{q}\Big{[}\sum_{j=1}^{p}\log\alpha_{j}+p\log(2\pi)-\beta^{T}A\beta\Big{]}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\Big{[}\sum_{j=1}^{p}\psi(\alpha_{j})-\sum_{j=1}^{p}\log
b_{j}+p\log(2\pi)$
$\displaystyle-\text{tr}\big{[}(\Sigma_{\beta}+\mu_{\beta}\mu_{\beta}^{T})\text{diag}(\frac{a_{j}}{b_{j}})\big{]}\Big{]},$
$\displaystyle E(\log P(\alpha))=$
$\displaystyle\sum_{j=1}^{p}E_{q}\left[a_{0}\log
b_{0}-\log\Gamma(a_{0})+(a_{0}-1)\log\alpha_{j}-b_{0}\alpha_{j}\right]$
$\displaystyle=$ $\displaystyle p\big{[}a_{0}\log
b_{0}-\log\Gamma(a_{0})\big{]}+\sum_{j=1}^{p}\big{[}(a_{0}-1)(\psi(a_{0})$
$\displaystyle-\log(b_{j}+\frac{d_{jj}}{2}))-b_{0}\frac{a_{j}}{b_{j}+\frac{d_{jj}}{2}}\big{]},$
$\displaystyle E(\log P(\gamma|\pi))=$
$\displaystyle\sum_{j=1}^{p}E_{q}\left[\gamma_{j}\log\pi_{j}+(1-\gamma_{j})\log(1-\pi_{j})\right]$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{p}\left[P_{j}\left[\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right]+(1-P_{j})\left[\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right]\right],$
$\displaystyle E(\log P(\pi))=$
$\displaystyle\sum_{j=1}^{p}(\alpha_{0}+1)\left[\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right]$
$\displaystyle+\sum_{j=1}^{p}(\beta_{0}-1)\left[\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})\right],$
$\displaystyle E(\log q(\beta))=$
$\displaystyle-\frac{1}{2}\big{(}\log\big{(}\det(\Sigma_{\beta})\big{)}+p\log(2\pi)+p\big{)},$
$\displaystyle E(\log q(\gamma))=$
$\displaystyle\sum_{j=1}^{p}E_{q}\left[\gamma_{j}\log
P_{j}+(1-\gamma_{j})\log(1-P_{j})\right]$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{p}P_{j}\log P_{j}+(1-P_{j})\log(1-P_{j}),$
$\displaystyle E(\log q(\pi))=$
$\displaystyle\sum_{j=1}^{p}\log\Gamma(\alpha_{j}^{\ast}+\beta_{j}^{\ast})-\log\Gamma(\alpha_{j}^{\ast})-\log\Gamma(\beta_{j}^{\ast})$
$\displaystyle+(\alpha_{j}^{\ast}-1)(\psi(\alpha_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast}))$
$\displaystyle+(\beta_{j}^{\ast}-1)(\psi(\beta_{j}^{\ast})-\psi(\alpha_{j}^{\ast}+\beta_{j}^{\ast})),$

and

$\displaystyle E(\log q(\alpha))=$
$\displaystyle\sum_{j=1}^{p}(a_{j}-1)\psi(a_{j})+\log
b_{j}-a_{j}-\log\Gamma(a_{j}).$

## Appendix (B): Additional plots

Here, we give some additional plots.

Figure 11: Posterior bar-plots of the Bernoulli mass functions of the $\gamma$
for Bernoulli product prior model.

Figure 12: P posterior bar-plots of the Bernoulli mass functions of the $z$
for continuous spike prior model.

Figure 13: Approximate posterior density functions of the beta posteriors of
the hyper-parameter $\pi$ for Bernoulli product prior and continuous spike
prior models.

Figure 14: Approximate (normalized) posterior density functions of the inverse
gamma posteriors of the variance hyper-parameters for Bernoulli product prior
model.

Figure 15: Approximate (normalized) posterior density function of the inverse
gamma posterior of the variance hyper-parameter for continuous spike prior
model.

Figure 16: Approximate (normalized) posterior density functions of the inverse
gamma and generalized inverse Gaussian posteriors of the variance hyper-
parameters for Laplace prior model.