arXiv:2606.24357v3 Announce Type: replace-cross Abstract: Covariate selection in Generalized Linear Models (GLMs) is a fundamental problem in statistics, as including irrelevant predictors might lead to overfitting and poor interpretability, while omitting relevant ones might result in biased estimates. Most Bayesian approaches to variable selection -- including spike-and-slab priors and continuous shrinkage priors -- have key limitations, e.g., (i) are based on non fully conjugate formulations, (ii) are restricted to a linear model, or (iii) lack posterior consistency guarantees for the variable selection procedure and model parameters. In this work, we propose a fully Bayesian hierarchical and conjugate framework for covariate selection in GLMs, applicable to any distribution in the exponential family, based on modeling a binary inclusion indicator that directly encodes covariate inclusion in the linear predictor. In our approach, variable selection and parameter estimation are performed simultaneously, incorporating both sources of uncertainty in posterior inference. Consequently, our methodology provides a valid post-model Bayesian selection procedure. We present theoretical guarantees of the proposed fully conjugate Bayesian variable selection for GLMs, establishing posterior consistency of both the inclusion indicators and the active regression coefficients. We derive an efficient Gibbs Sampling algorithm with a corresponding R package implementation. We validate the proposed method on synthetic and real-world datasets, demonstrating competitive predictive and inferential performance.