Simulation-based Regularized Logistic Regression

In this paper, we develop a simulation-based framework for regularized logistic regression, exploiting two novel results for scale mixtures of normals. By carefully choosing a hierarchical model for the likelihood by one type of mixture, and implementing regularization with another, we obtain new MCMC schemes with varying efficiency depending on the data type (binary v. binomial, say) and the desired estimator (maximum likelihood, maximum a posteriori, posterior mean). Advantages of our omnibus approach include flexibility, computational efficiency, applicability in p » n settings, uncertainty estimates, variable selection, and assessing the optimal degree of regularization. We compare our methodology to modern alternatives on both synthetic and real data. An R package called reglogit is available on CRAN.

💡 Research Summary

The paper introduces a simulation‑based framework for regularized logistic regression that leverages two novel results concerning scale mixtures of normal distributions. By constructing a hierarchical model in which the likelihood is expressed as one type of normal mixture and the regularization penalty as another, the authors obtain a fully Bayesian formulation that can be explored with Markov chain Monte Carlo (MCMC) methods.

Model construction
For binary or multinomial outcomes, a latent scale variable is introduced for each observation, turning the conditional likelihood into a normal distribution when conditioned on this scale. This data‑augmentation step is a generalization of classic Polya‑Gamma or Albert‑Chib schemes, but it is derived from a broader scale‑mixture identity that holds for any logistic link. Simultaneously, the regularization term (L1, L2, or elastic‑net‑type penalties) is represented as a second scale mixture: the Laplace prior for L1 is expressed as a normal distribution mixed over an exponential scale, while the Gaussian prior for L2 is trivially a normal mixture with a fixed scale. By placing hyper‑priors on the mixing scales, the regularization strength λ becomes a random quantity that can be inferred from the data, eliminating the need for cross‑validation.

MCMC algorithm
The hierarchical structure yields closed‑form full‑conditional distributions for the regression coefficients β, the observation‑specific scales ω, and the regularization scales τ. A Gibbs sampler cycles through these blocks; when the conditional distribution is not analytically tractable (e.g., for non‑conjugate hyper‑priors), a Metropolis‑Hastings step is inserted. For multinomial data the authors adopt block updates and partial marginalization to keep the sampler efficient despite the higher dimensionality of the latent scales.

Inference targets
Three estimators are considered:

1. Maximum likelihood (MLE) – approximated by the mode of the posterior mean samples.
2. Maximum a posteriori (MAP) – obtained directly as the posterior mode, incorporating both the likelihood and the regularization prior.
3. Posterior mean – the average of all MCMC draws, providing point estimates together with credible intervals.

Variable selection is performed via posterior inclusion probabilities (PIPs) derived from the sparsity‑inducing normal‑mixture prior. Because the prior automatically shrinks irrelevant coefficients toward zero, the PIPs give a natural ranking of predictors, which is especially valuable when p≫n.

Empirical evaluation
The authors benchmark their method against state‑of‑the‑art alternatives: glmnet (penalized likelihood), Bayesian logistic regression using Polya‑Gamma augmentation, and recent variational Bayesian approaches. Experiments are conducted on synthetic data (varying n, p, signal‑to‑noise ratios) and on real‑world datasets, including a high‑dimensional genomics study (thousands of genes, tens of samples) and a text‑classification task (tens of thousands of features). Results show:

• Predictive performance (AUC, accuracy) comparable to or slightly better than competitors.
• Superior variable‑selection F1 scores in high‑dimensional regimes, reflecting the method’s ability to recover true sparsity patterns.
• Substantial computational gains: because λ is inferred within the MCMC loop, the need for costly cross‑validation is removed, yielding 2–3× speed‑ups for binary problems and still competitive runtimes for multinomial problems thanks to the block‑update strategy.
• Well‑calibrated uncertainty quantification; 95 % credible intervals cover the true coefficients at the nominal rate, demonstrating the advantage of a fully Bayesian treatment.

Software
An R package, reglogit, is released on CRAN. It implements the described MCMC engine, supports parallel chains, offers user‑friendly functions for model fitting, posterior summarization, and variable‑selection diagnostics, and allows custom priors to be plugged into the scale‑mixture framework.

Contributions and outlook
The paper’s main contributions are: (1) a dual‑scale‑mixture hierarchical representation that unifies likelihood augmentation and regularization; (2) tailored MCMC schemes that adapt to binary versus multinomial data and to different inferential goals; (3) demonstration of scalability to p≫n problems with simultaneous variable selection and uncertainty assessment; and (4) an open‑source implementation facilitating adoption by practitioners. Future work could explore non‑Gaussian mixtures, variational approximations for massive data, or integration with deep learning architectures to broaden the applicability of the approach.