The data augmentation algorithm

In modern datasets arising from diverse applications such as genetics, medical science, and other scientific disciplines, the number of covariates are often larger than the number of observed samples. The ordinary least squares method for estimating the regression coefficients in a linear regression is not applicable in such situations. On the other hand, penalized regression methods such as the least absolute shrinkage and selection operator (lasso) (Tibshirani,, 1996) are applicable as they can simultaneously induce shrinkage and sparsity in the estimation of the regression coefficients.

Consider the standard linear model

where $Z\in\mathbb{R}^{m}$ is the vector of responses, $\mu\in\mathbb{R}$ is the intercept, $1_{m}$ is the $m\times 1$ vector of 1’s, $W$ is the $m\times p$ (standardized) covariate matrix, $\beta=(\beta_{1},\dots,\beta_{p})\in\mathbb{R}^{p}$ is the vector of regression coefficients, and $\sigma^{2}$ is the variance parameter. The lasso is based on $L_{1}$ norm regularization, and it estimates $\beta$ by solving

for some shrinkage parameter $\lambda\in\mathbb{R}$, where $z$ is the vector of observed responses, and $\tilde{z}=z-{\bar{z}}1_{m}$. The lasso estimate can be interpreted as the posterior mode when, conditional on $\sigma^{2}$, the regression parameters $\beta_{j}$’s have independent and identical Laplace priors (Tibshirani,, 1996). Indeed, following Park and Casella, (2008), we consider the hierarchical model

where $\pi(\sigma^{2})$ is the prior density of $\sigma^{2}$. Since the columns of $W$ are centered, standard calculations show that

Thus, from (2.1) and (2.3), the joint posterior density of $(\beta,\sigma^{2})$ is

From (2.1) and (2.4) it follows that the mode of the (conditional on $\sigma^{2}$) posterior density of $\beta$ is the lasso estimate.

The posterior density (2.4) is intractable. Introducing augmented variables $y=(y_{1},\dots,y_{p})$ with $y_{i}>0$ for all $i$, Park and Casella, (2008) consider the following joint density of $(\beta,\sigma^{2},y)$

Replacing $t$ with $\beta_{j}/\sigma$, $s$ with $y_{j}$, and $a$ with $\lambda$ in the following representation of the Laplace density as a scale mixture of normals (Andrews and Mallows,, 1974)

we see that

that is, the $(\beta,\sigma^{2})-$ marginal density of the joint posterior density $\pi(\beta,\sigma^{2},y|z)$ given in (2.1) is the target posterior density (2.4). Thus, from Section 1, if sampling from the two conditional densities $\pi(\beta,\sigma^{2}|y,z)$ and $\pi(y|\beta,\sigma^{2},z)$ is straightforward, then we can construct a valid DA algorithm for (2.4).

From (2.1), we have

Recall that the Inverse-Gaussian $(\kappa,\psi)$ density is given by

Thus from (2.7), it follows that

In order to derive the conditional density $\pi(\beta,\sigma^{2}|y,z)$ we assume that apriori $\sigma^{2}\sim\mbox{Inverse-Gamma}(\alpha,\xi)$ for some $\alpha\geq 0$ and $\xi\geq 0$. The improper prior $\pi(\sigma^{2})=1/\sigma^{2}$ used in Park and Casella, (2008) is obtained by replacing $\alpha=0,\xi=0$. Note that a draw from $\pi(\beta,\sigma^{2}|y,z)$ can be made by first drawing from $\pi(\sigma^{2}|y,z)$ followed by a draw from $\pi(\beta|\sigma^{2},y,z)$ (Rajaratnam et al.,, 2019). Denoting the $p\times p$ diagonal matrix with diagonal elements $(y_{1},\dots,y_{p})$ by $D_{y}$, we have

Thus, $\beta|\sigma^{2},y,z\sim N_{p}((W^{\top}W+D_{y}^{-1})^{-1}W^{\top}\tilde{z},\sigma^{2}(W^{\top}W+D_{y}^{-1})^{-1})$. Also, from (2.1), we have

Thus, $\sigma^{2}|y,z\sim\mbox{Inverse-Gamma}((n-1+2\alpha)/2,(\tilde{z}^{\top}(I-W(W^{\top}W+D_{y}^{-1})^{-1}W^{\top})\tilde{z}+2\xi)/2).$ Hence, a single iteration of the DA algorithm uses the following two steps to move from $(\beta,\sigma^{2})$ to $(\beta^{\prime},{\sigma^{2}}^{\prime})$.

Although the lasso estimator has been extensively used in applications as diverse as agriculture, genetics, and finance, it may perform unsatisfactorily if the predictors are highly correlated. For example, if there is a group structure among the variables, lasso tends to select only one variable from each group. Zou and Hastie, (2005) proposed the Elastic Net (EN) to achieve better performance in such situations. The EN estimator is obtained by solving

where $\lambda_{1}$ and $\lambda_{2}$ are tuning parameters. From (2.9) we see that the elastic net uses both an $L_{1}$ penalty as in lasso and an $L_{2}$ penalty as in the ridge regression. Following Kyung et al., (2010) and Roy and Chakraborty, (2017), we consider the hierarchical Bayesian EN model

where $\pi(\sigma^{2})$ is the prior density of $\sigma^{2}$. From (2.1) and (2.3) it follows that the joint posterior density of $(\beta,\sigma^{2})$ is

From (2.9) and (2.11) we see that the mode of the (conditional on $\sigma^{2}$) posterior density of $\beta$ is the EN estimate. Since the density (2.11) is intractable, using augmented variables $y=(y_{1},\dots,y_{p})$, following Kyung et al., (2010) and Roy and Chakraborty, (2017), we consider the joint density of $(\beta,\sigma^{2},y)$ given by

From (2.6) we see that

From (2.1) using similar calculations as before, we see that the conditional distributions of $1/y_{j}$’s are the same as (2.8) with $\lambda=\lambda_{1}$. As before, we assume that apriori $\sigma^{2}\sim\mbox{Inverse-Gamma}(\alpha,\xi)$ for some $\alpha\geq 0$ and $\xi\geq 0$. Denoting the $p\times p$ diagonal matrix with diagonal elements $(1/(\lambda_{2}+1/y_{1}),\dots,1/(\lambda_{2}+1/y_{p}))$ by $\tilde{D}_{y}$, we have

Thus, $\beta|\sigma^{2},y,z\sim N_{p}((W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top}\tilde{z},\sigma^{2}(W^{\top}W+\tilde{D}_{y}^{-1})^{-1})$. Next, doing similar calculations as for the Bayesian lasso, we see that $\sigma^{2}|y,z\sim\mbox{Inverse-Gamma}((n-1+2\alpha)/2,(\tilde{z}^{\top}(I-W(W^{\top}W+\tilde{D}_{y}^{-1})^{-1}W^{\top})\tilde{z}+2\xi)/2).$ Hence, a single iteration of the DA algorithm for the Bayesian EN uses the following two steps to move from $(\beta,\sigma^{2})$ to $(\beta^{\prime},{\sigma^{2}}^{\prime})$.

Different Bayesian generalized lasso methods, such as the Bayesian group lasso, the Bayesian sparse group lasso, and the Bayesian fused lasso models have been proposed in the literature to handle the situations where the covariates are known to have some structures, for example, when they form groups or are ordered in some way (Kyung et al.,, 2010; Xu and Ghosh,, 2015). Introducing appropriate augmented variables $y$, DA algorithms for these models can also be constructed (Kyung et al.,, 2010; Jin and Tan,, 2021).

Logistic regression model is likely the most popular model for analyzing binomial data. Let $(Z_{1},Z_{2},\dots,Z_{m})$ be a vector of independent Binomial random variables, and $w_{i}$ be the $p\times 1$ vector of known covariates associated with $Z_{i}$ for $i=1,\dots,m$. Let $\beta\in\mathbb{R}^{p}$ be the vector of unknown regression coefficients. Assume that $Z_{i}\sim$ Binomial $(\ell_{i},F(w_{i}^{\top}\beta))$ where $F(t)=e^{t}/(1+e^{t})$ is the cumulative distribution function of the standard logistic random variable. Denoting the observed responses by $z=(z_{1},z_{2},\dots,z_{m})$, the likelihood function for the logistic regression model is

Consider a Bayesian analysis that employs the following Gaussian prior for $\beta$

where $\mu_{0}\in\mathbb{R}^{p}$ and $Q$ is a $p\times p$ positive definite matrix. Then the intractable posterior density of $\beta$ is given by

where

is the marginal density of $z$.

There have been many attempts (Holmes and Held,, 2006; Frühwirth-Schnatter and Frühwirth,, 2010) to produce an efficient DA algorithm for the logistic regression model, that is (2.14), without much success until recently, when Polson et al., (2013) (denoted as PS&W hereafter) proposed a new DA algorithm based on the Pólya-Gamma (PG) distribution. A random variable $\theta$ follows a PG distribution with parameters $a>0,\,b\geq 0$, if

where $g_{i}\overset{iid}{\sim}$Gamma$(a,1)$. From Wang and Roy, 2018c , the pdf for PG$(a,b),a>0,\,b\geq 0$ is

for $\theta>0,$ where the hyperbolic cosine function $\cosh(t)=(e^{t}+e^{-t})/2$. We denote this as $\theta\sim$PG$(a,b)$.

Let $y=(y_{1},y_{2},...y_{n})$, $k_{i}=z_{i}-\ell_{i}/2,\,i=1,...,n$ and $p(y_{i}\mid\ell_{i},0)$ be the pdf of $y_{i}\sim$PG$(\ell_{i},0)$. Define the joint posterior density of $\beta$ and $y$ given $z$ as

By Theorem 1 in Polson et al., (2013), it follows that the $\beta$-marginal density of the joint posterior density $\pi(\beta,y\mid z)$ is the target density $\pi(\beta\,|\,z)$ given in (2.14), that is,

PS&W’s DA algorithm for the logistic regression model is based on the joint density $\pi(\beta,y\mid z)$. We now derive the conditional densities, $\pi(\beta|y,z)$ and $\pi(y|\beta,z)$. From (2.17) we see that

and thus from (2.16) we have $y_{i}|\beta,z\overset{ind}{\sim}\text{PG}\left(\ell_{i},\left|w_{i}^{\top}\beta\right|\right),$ for $i=1,\dots,m$. Let $W$ denote the $m\times p$ design matrix with $i$th row $w_{i}^{\top}$ and $Y$ be the $m\times m$ diagonal matrix with $i$th diagonal element $y_{i}$. From (2.17) and (2.13), it follows that the conditional density of $\beta$ is

where $\kappa=\left(\kappa_{1},\dots,\kappa_{m}\right)^{\top}$. Thus the conditional distribution of $\beta$ is multivariate normal. In particular,

Thus, a single iteration of PS&W’s DA algorithm uses the following two steps to move from $\beta$ to $\beta^{\prime}$.

PS&W’s DA algorithm is also applicable to the situation when $Q=0$ in (2.13), that is, $\pi(\beta)\propto 1$, the improper uniform prior on $\beta$ provided the posterior density $\pi(\beta|z)$ is proper, that is $c(z)$ defined in (2.15) is finite. Polson et al., (2013) demonstrate superior empirical performance of their DA algorithm over some other DA and MCMC algorithms for the logistic regression model.

Generalized linear mixed models (GLMMs) are often used for analyzing correlated binary observations. The random effects in the linear predictor of a GLMM can accommodate for overdispersion as well as dependence among correlated observations arising from longitudinal or repeated measures studies. Let $(Z_{1},Z_{2},\dots,Z_{m})$ denote the vector of Bernoulli random variables. Let $w_{i}$ and $v_{i}$ be the $p\times 1$ and $q\times 1$ known covariates and random effect design vectors, respectively associated with the $i$th observation $Z_{i}$ for $i=1,\dots,m$. Let $\beta\in\mathbb{R}^{p}$ be the unknown vector of regression coefficients and $u\in\mathbb{R}^{q}$ be the random effects vector. The probit GLMM connects the expectation of $Z_{i}$ with $w_{i}$ and $v_{i}$ using the probit link function, $\Phi^{-1}(\cdot)$ as

where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal random variable.

Assume that we have $r$ random effects with $u=(u_{1}^{\top},\dots,u_{r}^{\top})^{\top}$, where $u_{j}$ is a $q_{j}\times 1$ vector with $q_{j}>0$, $q_{1}+\cdots+q_{r}=q$, and $u_{j}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}N(0,\Lambda_{j}\otimes R_{j})$ where the low-dimensional covariance matrix $\Lambda_{j}$ is unknown and must be estimated, and the structured matrix $R_{j}$ is usually known. Here, $\otimes$ indicates the Kronecker product. Denoting $\Lambda=(\Lambda_{1},\dots,\Lambda_{r})$, the probit GLMM is given by

Let $z=(z_{1},z_{2},\dots,z_{m})^{\top}$ be the observed Bernoulli response variables. Note that, the likelihood function for $(\beta,\Lambda)$ is

which is not available in closed form. Here, $A(\Lambda)=\oplus_{j=1}^{r}\Lambda_{j}\otimes R_{j}$, with $\oplus$ indicating the direct sum, and $\phi_{q}(s;a,B)$ denotes the probability density function of the $q-$dimensional normal distribution with mean vector $a$, covariance matrix $B$ and evaluated at $s$.