Simulation-based Regularized Logistic Regression

Robert B. Gramacy, Nicholas G. Polson

Introduction

Large scale logistic regression has numerous modern day applications from text classification to genetics. We develop a flexible framework for maximum likelihood, maximum a posteriori, and full Bayesian posterior inference for regularized models. Our motivations stem from a desire to find common ground between point estimation in “large- $p$ ” settings (Krishnapuram et al.,, 2005; Genkin et al.,, 2007), where $p$ is the number of predictors, and full Bayesian inference for “small- $p$ ” (Holmes and Held,, 2006; Frühwirth-Schnatter and Frühwirth,, 2007; Frühwirth-Schnatter et al.,, 2009; Fahrmeir et al.,, 2010; Frühwirth-Schnatter and Frühwirth,, 2010). Collecting such distinct methods into a unifying framework facilitates a number of novel enhancements including posterior inference for the amount of regularization, and an efficient handling of binomial data.

We start by framing a typical regularized optimization criteria as a powered-up posterior, or power-posterior (Friel and Pettitt,, 2008), with a shrinkage prior such as the lasso (Tibshirani,, 1996). We then show how inference may proceed by employing two (heretofore unrelated) data augmentation schemes: one for the powered-up logistic likelihood; and the other for the prior. The combined effect is a fully Gibbs MCMC sampler which, among other advantages, allows estimators previously requiring custom algorithms to be calculated via a single simulated annealing (Kirkpatrick et al.,, 1983) scheme.

Our approach offers a fully probabilistic alternative by viewing the objective function (1) as a (log) posterior distribution whose maximum a posteriori (MAP) estimator coincides with $\hat{\beta}$ . A multiplicity parameter $\kappa$ can then be introduced to help find the MAP via simulation. Our key insight, which makes the simulation efficient, is that the logistic likelihood component of the posterior can be written hierarchically using $z$ –distributions (Barndorff-Nielsen et al.,, 1982), leading to a data augmentation scheme that generalizes that of Holmes and Held, (2006) [HH hereafter]. Combining this with a standard data augmentation for the prior yields a highly blocked Gibbs MCMC algorithm for logistic regression. $Z$ -distributions also suggest a new representation of the likelihood that is equivalent (to HH) but requires $n$ fewer latent variables. Finally, we recognize that $\kappa$ has a secondary use for binomial data (multiple $y$ observed for each $x$ ) which otherwise would require more latent variables.

A distinctive feature of our framework is how it deals with the amount of regularization, $\nu$ , which is traditionally chosen by cross validation (CV). As an alternative, we may extend the hierarchical model to include a prior for $\nu$ so that the marginal likelihood can be computed and used to set $\nu=\hat{\nu}$ , or to integrate $\nu$ out. Posterior expectations, thus obtained, can give superior point–estimators for $\beta$ in large- $p$ linear regression contexts (Hans,, 2009), and we show how this extends to logistic regression.

The rest of the paper is outlined as follows. Section 2 provides our data augmentation strategies for sparse high dimensional logistic regression, and Section 3 develops an MCMC scheme for estimation. Section 4 illustrates our approach with empirical comparisons to modern competitors. Finally, Section 5 concludes with simple extensions and directions for future research. An supporting R package, reglogit, is available on CRAN.

Regularized logistic regression via power-posteriors

The central problem is to find the MLE, MAP, or posterior mean estimator in logistic regression. To do this, consider the following power-posterior distribution inspired by Eq. (1):

Observe that the likelihood–prior combination below yields Eq. (2) via Bayes’ rule.

The following subsections provide data augmentation schemes for this likelihood and prior. They primarily concentrate on the $\alpha=1$ case, i.e., the double–exponential or lasso prior, although results are developed in generality when possible. Section 5 briefly touches on the simpler $\alpha=2$ case.

Extending a well-known technique for generating logistic regression (e.g., Andrews and Mallows,, 1974; Holmes and Held,, 2006), we represent the powered-up likelihood (3) for $\beta$ as a marginal quantity obtained after integrating over latent variables $(z,\lambda)$ , where $z=(z_{1},\ldots,z_{n})$ and $\lambda=(\lambda_{1},\ldots,\lambda_{n})$ . That is, each element of the product of independent logistic terms can be written as a two-dimensional integral:

This suggests a hierarchical representation in terms of latent variables, $z_{i}$ for each $y_{i}$ , mixed over $\lambda_{i}$ . It remains to determine the appropriate form of $p_{\kappa}(z_{i}|\beta,\lambda_{i},y_{i})$ and $p_{\kappa}(\lambda_{i})$ so that $(1+e^{-y_{i}x_{i}^{\top}\beta})^{-\kappa}=\int\!\int p_{\kappa}(z_{i}|\beta,\lambda_{i},y_{i})p_{\kappa}(\lambda_{i})\,d\lambda_{i}\,dz_{i}$ .The notation reserves $\pi(\cdot)$ for the marginal posterior $\beta$ as a visual queue for the quantity of primary interest. All other probability densities use $p(\cdot)$ , including the joint for latent $(z,\lambda)$ and all priors.

Our key result, generalizing HH, relies on a scale mixture representation of $z$ –distributions (Barndorff-Nielsen et al.,, 1982). These are characterized by their pdf as:

where $q_{a,b}(\lambda)$ is a Polya distribution, i.e., an infinite sum of exponentials:

and the weights $w_{k}$ are determined via $\delta=(a+b)/2$ and $\theta=(a-b)/2$ as

The (powered up) logistic function may be represented as follows.

If $z\sim Z(1,\kappa,1,yx^{\top}\beta)$ , then $F_{Z}(z)=1-(1+e^{z-yx^{\top}\beta})^{-\kappa}$ , giving $1-F_{Z}(0)=(1+e^{-yx^{\top}\beta})^{-\kappa}$ . In other words,

establishing the outer integration, over $z$ , in Eq. (9). Applying the representation in Eq. (5) yields the desired result. ∎

The statistical implication of this is a hierarchical model which we summarize in the following corollary.

The conditional distribution $p_{\kappa}(z_{i}|\beta,\lambda_{i},y_{i})$ and the mixing distribution $q_{1,\kappa}(\lambda_{i})$ imply that the latent $z_{i}$ follow

where $\mathcal{N}^{+}$ is the normal distribution truncated to the positive real line.

When $\kappa>1$ , the asymmetry of the $z$ –distribution makes it harder to extract $y_{i}$ from $y_{i}x_{i}^{\top}\beta+\frac{1}{2}(1-\kappa)\lambda_{i}$ , the mean of the truncated normal in Eq. (11). In Section 3.3, we indirectly suggest that one can interpret $\kappa y_{i}$ as a binomial response when $\kappa$ is an integer.

Theorem 9 shows how components of the powered-up logistic likelihood can be represented hierarchically by the cdf of $z$ –distributions. We therefore call that multiplicity extension to HH the cdf representation. However, further inspection reveals that it is possible to eliminate an integral in Eq. (4) and thus $n$ latent variables, and use the representation

which avoids integrating over $z_{i}$ . Instead, set them to zero (and $\mu_{i}=y_{i}x_{i}^{\top}\beta$ ) and directly obtain $(1+e^{y_{i}x_{i}^{\top}\beta})^{-\kappa}$ . By analogy, we call this a pdf representation as it involves evaluating a particular $z$ -density function. This simple representation is problematic, however, since the Polya mixing density $q_{0,\kappa}$ is improper. In particular, note that $\psi_{0}=0$ , resulting in a infinite weight in the generative formulation (8).

Fortunately, a similar representation may be generated

which involves a proper Polya mixing density as long as $(a,b)>0$ and $a+b=\kappa$ . In Section 3.1 & 4, we show how the extra $e^{a\mu}\equiv e^{ay_{i}x_{i}^{\top}\beta}$ poses no problem for efficient inference, and that $(a=\frac{1}{2},b=\kappa-\frac{1}{2})$ works well in practice. But first, we complete the power-posterior specification with a family of regularization priors on $\beta$ .

2 Prior regularization

Box and Tiao, (1973) provide a general discussion of (14) in the linear regression context. Some notable special cases in the recent literature on sparse logistic regression include the following: when $\nu=1$ , $\alpha=1$ , and $\sigma_{j}=\lambda_{j}$ it is the Laplace prior used in Genkin et al., (2007); when $\alpha=2,\sigma_{j}=1$ and $\nu=\sigma^{2}$ it is the Gaussian prior, and when $\alpha=2,\sigma_{j}=1$ and $\nu^{-1}=\lambda$ it is the Laplace prior from Krishnapuram et al., (2005).The $\lambda_{j}$ and $\lambda$ variables correspond to the shrinkage parameters so named in our references. They should not be confused with the latent $\lambda_{i}$ used in our hierarchical likelihood representation. Inference for $\nu$ in these cases typically proceeds by CV, or by inspecting the paths of $\hat{\beta}_{\nu}$ solutions for varying $\nu$ . Assessing the uncertainty in estimators $\hat{\beta}_{\hat{\nu}}$ on the final choice of $\hat{\nu}$ can pose difficulties.

The prior in Eq. (14)—for the purposes of efficient inference [Section 3]—is an adaptation of a scale mixture of normals result from West, (1987) to account for $\kappa$ . Specifically,

Simulation-based logistic regression

We develop a Gibbs sampling algorithm [Section 3.2] for sampling the augmented power-posterior $p_{\kappa}(\beta,z,\omega,\lambda,\nu|y,\sigma^{2})$ , for any $\kappa$ . We first derive the relevant posterior conditionals [Section 3.1], treating cdf and pdf representations in turn. When $\kappa=1$ the marginal samples of $\beta$ summarize the posterior distribution of the main parameters of interest. Obtaining the MAP or MLE requires an inhomogeneous Markov chain [Section 3.2]. Finally, we describe how a vectorized $\kappa$ can facilitate efficient Bayesian binomial regression [Section 3.3].

To begin, consider the latent $z$ and $\lambda$ variables in the cdf and pdf representations, in turn, followed by the coefficients $\beta$ and corresponding regularization prior parameters $(\omega,\nu)$ .

By construction [Eq. (11) of Corollary 1], the posterior full conditional for the latents, $p_{\kappa}(z_{i}|\beta,\lambda_{i},y_{i})$ , is a truncated (non-negative) normal distribution. Obtaining samples, independently for $i=1,\dots,n$ , is straightforward following the methods of Robert, (1995).

Sampling from the full conditional $p_{\kappa}(\lambda_{i}|\beta,z_{i},y_{i})$ is complicated by the infinite sum in the expression for the prior (6), which precludes a naïve approach via truncation since certain combinations of $\lambda_{i}$ and $b\equiv\kappa$ can give highly inaccurate, even negative, evaluations. HH derive an expression for this conditional when $\kappa=1$ and provide a rejection sampling algorithm by squeezing (Devroye,, 1986). Although adaptable for general $\kappa$ , we prefer a Rao–Blackwellized approach. Interchanging the order of integration in Eq. (9) suggests a corollary to Theorem 9 that is helpful in constructing a Metropolis–Hastings (MH) scheme for obtaining $\lambda_{i}$ draws.

The following is an alternate integral representation of the logistic function

where $\Phi$ is the cdf of the standard normal distribution.

Proposals $\lambda_{i}^{\prime}\sim q_{1,\kappa}(\lambda)$ can then be accepted via MH with probability $\min\{1,A_{i}\}$ where

Good proposals may be obtained by truncating the sum in Eq. (8) at $K=100$ for $\kappa=b=1$ , with improvements for larger $\kappa$ . Direct sampling is also possible (e.g., Weron,, 1996).

Empirically, the MH acceptance rate is high ( $>$ 90%) for $\kappa=1$ because posterior is similar to the prior ( $q_{1,1}$ ). Therefore the MH scheme may be preferable to the rejection/squeezing method of HH who report acceptance rates as low as 25%. Both rates decline as $\kappa$ is increased, but the MH rate is still above $1\%$ for $\kappa=20$ . A good rule of thumb is to thin $\lceil\kappa\rceil$ draws for each draw saved, which is reasonable from a computational standpoint as sampling from $q_{a,b}$ is fast. Even when thinning more than 10-fold, the MH sampler is competitive to HH/Devroye in terms of sheer speed. The MH requires two $\Phi$ evaluations, a few arithmetic operations, and two square roots. HH/Devroye, by contrast, can perform dozens (or more) expensive operations such as pow before the squeeze is made.Finally, drawing $\lambda_{i}$ unconditional on $z_{i}$ yields lower autocorrelation in the overall joint MCMC sampling scheme.

The pdf representation is simpler since $z_{i}$ is set to zero. Proposed $\lambda_{i}^{\prime}\sim q_{a,b}$ may be accepted or rejected via MH by exchanging a cdf for a pdf in Eq. (16) and replacing $\frac{1}{2}(1-\kappa)$ with $\frac{1}{2}(a-b)$ . Another feature that works well for the pdf representation is an adaptation of the slice sampler of Godsill, (2000). Given $\lambda_{i}$ , the next sample $\lambda_{i}^{\prime}$ may be obtained via an auxiliary uniform random variable as follows. Let $\phi_{i}\equiv\phi\{(-y_{i}x_{i}^{\top}\beta+\frac{1}{2}(a-b)\lambda_{i})/\sqrt{\lambda_{i}}\}$ , where $\phi$ is the pdf of a standard normal distribution. Then sample

where the second step is facilitated by accept/rejects following random draws from the Polya mixing density. Although more automatic in that it does not require thinning, we show in Section 4.2 that the MH scheme is faster overall. The two methods behave similarly when $\kappa$ gets large, causing the rate of rejections/required thinning to increase.

Regularized regression coefficient parameters (β,ω,ν)𝛽𝜔𝜈(\beta,\omega,\nu)

The full conditional distribution of each latent $\omega_{j}$ is proportional to the integrand of Eq. (15). When $\alpha=1$ we have the following adaptation of a standard result.

From the integrand in Eq. (15) with $\alpha=1$ we have

Our IG priors for $\nu$ are both conditionally conjugate. An IG prior for $\nu^{2}$ and the representation in Eq. (15) gives

extending the analysis of Park and Casella, (2008).

2 Gibbs sampling and annealing for point estimators

A full Gibbs sampling algorithm for both cdf and pdf representations is outlined in Figure 1.

The samples on output may be used to approximate expectations under the power-posterior distribution with multiplicity $\kappa$ . If $\kappa=1$ then these are samples from a well-defined posterior distribution which may be used, e.g., to approximate the posterior mean of $\beta$ or provide samples from the posterior predictive distribution. Both take into account the full the uncertainties of all parameters (including $\nu$ ) into account—a feature unique to full Bayesian analysis.

Settings of $\kappa>0$ are useful for finding other popular estimators via simulated annealing (SA). In our context, SA establishes an inhomogeneous Markov chain over a sequence of power-posteriors, starting with $\kappa=1$ and then increasing according to a pre-determined schedule. Except when Gibbs sampling is possible for all $\kappa$ (as for our power-posterior), it is usually difficult to ensure that the Markov chain mixes well, particularly when $\kappa$ increases. A pragmatic approach starts at $\kappa\approx 1$ , and systematically makes modest increases in $\kappa$ until Monte Carlo variation in the power-posterior expectations of the quantities of interest is below a pre-determined threshold. Each annealing iteration is initialized with the last value $\beta^{(S)}$ , $\nu^{(S)}$ , $\lambda^{(S)}$ and $z^{(S)}$ , from the previous iteration, thereby stitching the inhomogeneous Markov chains together. The chain for each $\kappa$ must have enough iterations to establish convergence to its particular power-posterior.

Annealed procedures such as ours present an MCMC alternative to EM-style algorithms. Importantly, SA is known to converge to the global optima in certain conditions (when $\kappa\rightarrow\infty$ ), whereas EM is only guaranteed to find a local optima. Although convergence for EM is usually quick, there are no guarantees that it will be so and indeed there are examples, particularly in high dimensional settings, where convergence can be arbitrarily slow. SA however, comes with the burden of choosing the schedule for increasing $\kappa$ . We have found that for our regularized logistic regression scheme, convergence is fast and mixing so good that short schedules such as $\kappa=1,5,10,20$ are a safe default [see Section 4]. Even jumping immediately to modest $\kappa\;(\approx 20$ ), skipping $\kappa=1$ , can very often yield cheap and accurate approximations.

3 Efficient handling of binomial data

Applications

The Pima Indian diabetes data [UCI Machine Learning Repository (Asuncion and Newman,, 2007)] includes outcomes for diabetes tests performed on $n=768$ women of Pima heritage with 8 real-valued predictors. Some of the predictors have many zeros, which may reasonably be interpreted as “missing” values. To remain consistent with the treatment of this data by HH, and other authors, we do not treat these values in any special way. The following analysis highlights properties of regularized estimators of $\beta=(\beta_{0}\equiv\mu,\beta_{1},\dots,\beta_{8})$ obtained with $\alpha=1$ , $\sigma_{j}=1$ for $j=1,\dots,8$ , and $T=1000$ samples from the resulting posterior (the first 100 as burn-in).

Figure 2 summarizes the marginal power posterior(s) for $\beta$ with boxplots. Three settings of $\kappa\in\{1,5,20\}$ (each panel) were used, and heavy regularization (fixing $\nu=6$ ) was applied. Only the first panel ( $\kappa=1$ ) summarizes samples from the true posterior. The $\kappa>1$ settings are useful for obtaining other estimators. The MLE, obtained from the glm command in R (R Development Core Team,, 2009), and the MAP as estimated from the sample(s), are also shown. Shrinkage is apparent in the divergence between the MAP and MLE values in all panels. Observe how the quartiles and outliers converge on the MAP as $\kappa$ is increased, reflecting higher confidence in the accurate estimation of those values. Convergence is particularly rapid for the intercept term, and the two coefficients with considerable mass near zero ( $\beta_{4}$ and $\beta_{5}$ ). These columns of $X$ have the highest concentration of “missing” values (30% and 49% respectively), so it is not surprising the that MAP estimator excludes them.

Figure 3 illustrates how mass concentrates on the MAP in two disparate cases for varying values of $\kappa$ . For $\beta_{2}$ (left panel), which is decidedly non-zero in the power posterior(s), the convergence to the MAP (apparently around $\beta_{2}=6$ ) is modest. In the case of $\beta_{4}$ (right panel) the convergence to the MAP (to zero) is more rapid as $\kappa$ is increased, allowing for confident variable de-selection in a way similar to the lasso for linear regression.

Finally, we consider the case where $\nu$ is also inferred by MCMC, jointly with the other parameters in the model. We use the IG prior on $\nu$ with $(r=2,d=0.1)$ , a typical default choice for linear regression (e.g., Gramacy and Pantaleo,, 2010).

Figure 4 shows the marginal posterior for $\nu$ under our settings of $\kappa$ . The rate of convergence is modest, with the spread of samples in the $\kappa=20$ case being only half that of the $\kappa=1$ case.

2 Comparing c/pdf representations on binomial data

Table 1 compares four different implementations of regularized binomial logistic regression ( $\alpha=1$ ) based on the output of 100 repeated experiments with $\sum n_{i}=2000$ (i.e., $m=100$ distinct $x_{i}$ predictors). The metrics for comparison are root mean squared error (RMSE) between the true and posterior mean $\beta$ s, and overall computing time of the respective MCMC samplers. In all cases, we use $T=1000$ MCMC rounds with MH sampling of $\lambda_{i}$ at thinning level(s) set by $\kappa^{\prime}$ (i.e., via $\kappa^{\prime}_{i}$ for each $\lambda_{i}$ ) as described in Section 3.1. The first 100 rounds were discarded as burn-in. The left table shows that there is no significant difference between the cdf and pdf representations, or between the flattened or multiplicity handling of binomial data, in terms of RMSE. The right table portrays a more interesting story in terms of CPU times. The many fewer latent variables needed by the multiplicity implementation leads to a much (9x) faster execution compared to flattening, with no cost in accuracy (via RMSE). In contrast, there is no speed gain to using $n$ fewer latent $z_{i}$ variables in the pdf representation.

Figure 5 illuminates the differences in behavior between the MH and slice sampler for the $\lambda_{i}$ draws (in the pdf representation). A particularly “sticky” case, as chosen from output of the experiment, had $\kappa^{\prime}_{i}=14$ . The top panel shows that many proposals from $q_{a,b}$ can be rejected under the MH ratio, even when the chain is automatically thinned. The bottom panel shows the chain obtained for the same $\lambda_{i}$ under the slice sampler, which never saves any rejected draws. However, this comes at the expense of many rejections in the inner–loop of the slice, resulting in a slow overall sampler. The median was four, but the mean was 81 owing to a heavy right-hand tail in the distribution of rejections whose central 95% quantile spanned to 114 and maximum reached 140,600. The overall MCMC scheme based on the slice sampler took four times longer than the one based on MH. Despite the absence of rejections, the mixing in slice sampler chain (assessed visually) was no better than MH. Indeed, their effective sample size due to autocorrelation (Kass et al.,, 1998) was nearly identical: 223 for slice sampling, and 221 for MH. Therefore, MH is recommended for speed considerations.

3 A simulated p≫nmuch-greater-than𝑝𝑛p\gg n experiment

We turn now to a predictive comparison of the methods of this paper, both fully Bayesian and full/joint MAP (including $\nu$ ), benchmarked against other modern approaches to regularized logistic regression. Consider a synthetic data experiment like the one in Section 4.2 except: $n_{i}=5$ for each of 20 unique predictors $x_{i}$ , so that $\sum n_{i}=100$ . Three variations on the data-generating $\beta$ vectors were used. In the first case $p=9$ and $\beta=(2,-3,0.74,-0.9,0,0,0,0)^{\top}$ ; in the second case $p=100$ , augmenting $\beta$ from the first case with 91 more zeros; and in the third $p=1000$ with 900 more zeros still. Each experiment involves a new random training design in the unit $p$ -cube. Random testing set are created similarly, except that $n_{i}^{\prime}=100$ so $\sum n_{i}^{\prime}=10000$ . The metrics of comparison are (approximated) expected log likelihood (ELL)Specifically, the average of $(1-p_{i})\log(1-\hat{p}_{i})+p_{i}\log\hat{p}_{i}$ over all testing locations $i$ , where $p_{i}$ and $\hat{p}_{i}$ are the true and estimated predictive probabilities of the first label, respectively. and misclassification rates.

Fully Bayesian posterior mean estimators (i.e., $\kappa=1$ ) are derived via priors/MCMC exactly as described in the preceding sections with $(100,1000)$ , $(500,1500)$ , $(1000,2000)$ burn-in and total MCMC rounds in each of the cases $p=9,100,1000$ , respectively. MAP estimators are found by running a $\kappa=10$ chain initialized at $(\beta,\lambda,\nu)$ -values from the $\kappa=1$ chain used for the mean estimators, except in the $p=1000$ case where $\nu$ was fixed to its posterior mean for reasons laid out in Section 3.2. Comparators include: the MLE obtained via the glm command in R; a binomial fit from the glmnet package (Friedman et al.,, 2010); and the estimator of Krishnapuram et al., (2005)This is equivalent to the Genkin et al., (2007) estimator but computationally less efficient. [“krish” for short]. The MLE was unstable in the $p=100$ & $1000$ cases, so these results were omitted. CV was used to choose the penalty parameter in the $p=9$ & $100$ cases for glmnet, via cv.glmnet. The same procedure gave fatal errors in the $p=1000$ case so we plugged in the estimate obtained from the corresponding $p=100$ run in for this final case. Reliably setting the penalty parameter for “krish”, via CV or otherwise, was too computationally intensive for the $p=100,1000$ cases so we picked a setting by hand using out-of-sample simulations from the $p=9$ case.

The results of the Monte Carlo experiment are summarized in Figure 6 by boxplots, and numerically. The best estimators have high ELL, low miss rates, and lower variability across the 100 repetitions. The fully Bayesian and “krish” methods are the best when $p=9$ (left-hand region of the boxplots and the top region of the tables). The former wins by ELL, having fewer low values, and the latter wins on miss rate, having more small ones. The “krish” method wins by both metrics on average, since it employs a fortuitously hand-chosen setting of the penalty parameter. The MLE is good on average, but has some extreme ELL and miss rate values. The glmnet and MAP estimators are positioned in between.

Distinctions in performance between the methods increases with $p$ . See the right-hand regions of the boxplots and the bottom regions of the tables. The “krish” method suffers from high variability due to the fixed choice of the penalty parameter. The glmnet variability is much lower, but there are many extreme outliers. Behavior in both $p=100$ and $1000$ cases is qualitatively similar for this estimator even though the former used CV to set the penalty parameter and the latter used the same fixed value. The MAP and fully Bayesian estimators have similar average behavior compared to other estimators, but with lower variability. Apparently, choosing the penalty parameter via the posterior offers the most stability in high dimensional settings. The fully Bayesian approach appears preferable to the MAP in all cases, but this distinction is harder to make out as $p$ increases.

4 Spam data with interactions

For a similar real-data experiment, consider the Spambase data set from UCI. It contains the binary classifications of 4601 emails based on 57 attributes which are treated as predictors. An interaction-expanded version of the predictor set contains approximately 1700 predictors. We performed a Monte Carlo experiment comprising of 20 random 5-fold CV training and testing sets using both the original and expanded predictors. Estimators were fit on the 100 training sets, and validated by misclassification rate on the testing ones. The Bayes estimators used (500,1500) MCMC (burn-in, total) rounds with the original 57 attributes, and (1000,2000) with the expanded set. The MAP and glmnet calculations were exactly as described for the $p=100$ case in Section 4.3 for the original predictor set, and like the $p=1000$ case for the expanded one. And “krish” was like $p=9$ and $p=100$ , respectively.

The results of the experiment are summarized in Figure 7. The first thing we notice is that, in contrast with the results in Section 4.3, the performance improves as the predictor set expands since some of the interaction terms make good predictors. The MLE is unstable, and so the regularized estimators offer an improvement even when the number of predictors is small relative to the number of instances. The Bayesian methods unilaterally outperform glmnet, and using the posterior to set the value of the regularization parameter is important in high dimensional settings. The “krish” estimator with fortuitous regularization is the best on the original predictor set, but worst on the expanded one where a revised setting of regularization could not be automated efficiently.

Discussion and extension

We provide a simulation-based approach to regularized logistic regression that facilitates a variety of inferential goals under a single framework. Most of the development of the methodology, and all of the applications, involved the $\alpha=1$ case. Everything extends to the ridge prior ( $\alpha=2$ ), i.e., an independent normal prior for each coefficient $\beta_{j}$ with variance $\sigma_{j}^{2}\nu^{2}/\kappa$ . Then, $p_{\kappa}(\omega_{j}|\beta,\nu)$ is a point mass at $\omega_{j}=1$ . Thus similar conjugacy results hold for the gamma prior on $\nu$ and $\nu^{2}$ .

From a computational perspective, our methods are competitive with the state-of-the art in un-regularized (and $\kappa=1$ ) contexts too. For example, we compared the efficiency of our methods to the “dRUM” MH sampler described by Frühwirth-Schnatter and Frühwirth, (2010). This method is attractive because it is fast and easy to implement. For example, on the Pima data it takes about 32s to generate 10,000 samples from the posterior which is about 7x faster than our pdf representation, which took 230s. However, the MH acceptance rate of the dRUM method was 46% which lead to an marginal ESS of 957 averaged over the nine $\beta_{j}$ coefficients. Our pdf representation had an average ESS that was about 5x better, at 4518. So the methods work out to have similar overall efficiences in that example. But in higher dimension like the 57-d spam data, our Gibbs sampling approach is much more attractive. The acceptance rate for dRUM was extremely low at 0.4%, which leads to ESSs that are essentially nil. Although our pdf representation is (again) 7x slower, faster convergence due to better movement in the chain leads to reasonable ESSs around 500.

There are several extensions of our methodology that readily present themselves. For example, handling polychotomous data (i.e., $>2$ classes) is straightforward. Following the setup in HH we may introduce $C$ collections of coefficients $\beta^{(1)},\dots,\beta^{(C)}$ for $C$ classes with the convention that $\beta^{(C)}=0$ so that logistic regression is recovered in the $C=2$ case. Then, we simply work with the conditional likelihoods $L(\beta^{(j)}|y,\beta^{(-j)})$ which turn out to have exactly the form of a logistic regression likelihood for the class indicator that each $y_{i}=j$ , independently for $i=1,\dots,n$ . If there are $n_{i}>1$ trials for predictors $x_{i}$ , then our algorithm for binomial logistic regression is applicable via a vectorized multiplicity parameter as described in Section 3.3. Extending the methods to ordinal responses is even easier. Johnson and Albert, (1999, Chapter 4) describe a Bayesian probit model which may be adapted for the logit case following either HH or our cdf representation. The pdf representation may not be readily applicable because the latent $z_{i}$ are useful for efficient sampling of the ordinal break points.

An further direction is to other classes of regularization priors. Implementing the Normal–Gamma extension (Griffin and Brown,, 2010) requires adding an extra (conjugate) parameter. A promising new approach is the horseshoe prior (Carvalho et al.,, 2010), which can be implemented with the addition of a slice sampler. Often variable selection is a primary goal of regularization, for which our methods would require further extension. For example, HH describe an approach to variable selection for logistic regression via Reversible Jump MCMC (Green,, 1995) which is adaptable to our framework. A similar regularized approach in a linear regression is provided byGramacy and Pantaleo, (2010). For variable selection for logistic regression using spike-and-slab priors, see Tüchler, (2008).

This research was partially funded by EPSRC grant EP/D065704/1 to RBG. The authors would like thank Matt Taddy for interesting discussions on the efficient handling of Binomial data, extensions to Multinomial regression, and EM code for the MAP estimator(s). We are grateful to two referees and an associate editor for valuable comments.

Appendix A Posterior conditional for β𝛽\beta in the pdf representation

For a particular $\lambda$ , i.e., ignoring the integral in Eq. (13), we have the following expression for the likelihood in vector/matrix form.

An expression for the posterior conditional for $\beta$ can then obtained by multiplying by the kernel of the MVN prior given $\omega$ , provided below Eq. (15), namely: $\exp\{-\frac{1}{2}\beta^{\top}\!(\frac{\kappa^{2}}{\nu^{2}}\Sigma^{-1}\Omega^{-1}\beta)\}$ . Combining the terms in the three exponents gives the following quadratic form.