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 methods are also fundamental building-blocks in more complex formulations, covering density regression, nonparametric classification, graphical models, and others. Within the Bayesian setting, inference typically proceeds by updating the Gaussian priors for the coefficients with the likelihood induced by probit or logit regressions for the binary responses. In this updating, the apparent absence of a tractable posterior has motivated a variety of computational methods, including Markov Chain Monte Carlo routines and algorithms which approximate the posterior. Despite being routinely implemented, current Markov Chain Monte Carlo methodologies face mixing or time-efficiency issues in large p and small n studies, whereas approximate routines fail to capture the skewness typically observed in the posterior. This article shows that the posterior distribution for the probit coefficients has a unified skew-normal kernel, under Gaussian priors. This result allows fast and accurate Bayesian inference for a wide class of applications, especially in large p and small-to-moderate n studies where state-of-the-art computational methods face substantial issues. These notable advances are quantitatively outlined in a genetic study, and further motivate the development of a wider class of conjugate priors for probit regression along with a novel independent sampler from the unified skew-normal posterior.