Conjugate Bayes for probit regression via unified skew-normals

Regression models for dichotomous data are ubiquitous in statistics. Besides being useful for inference on binary responses, such representations are also relevant building-blocks in more complex formulations, covering density regression, nonparametric classification and graphical models. Within a Bayesian framework, inference typically proceeds by updating the Gaussian priors for the regression coefficients with the likelihood induced by a probit or logit model for the observed binary responses. The apparent absence of conjugacy in this updating has motivated several computational methods, including Markov Chain Monte Carlo (MCMC) routines and algorithms which approximate the posterior. Despite being routinely implemented, current data augmentation MCMC methods tend to face mixing or time-efficiency issues in imbalanced high-dimensional data sets, whereas approximate routines fail to capture the skewness and heavy tails typically observed in the posterior distribution. This article shows that the posterior for the coefficients of a probit regression coincides with a unified skew-normal, under Gaussian priors. Due to this result, it is possible to derive explicitly the posterior distribution along with the predictive probability mass function of the response data and the marginal likelihood for model selection. A class of conjugate priors for Bayesian probit regression, improving flexibility in prior specification without affecting tractability in posterior inference, is also provided along with novel sampling and inference methods.