Conjugate Bayes for probit regression via unified skew-normal distributions
Daniele Durante
















Daniele DuranteConjugate Bayes for probit regression via unified skew-normals





Department of Decision Sciences and Bocconi Institute for Data Science and Analytics, Bocconi University, Via Röntgen 1, 20136, Milan, Italy daniele.durante@unibocconi.it 


Conjugate Bayes for probit regression 
 via unified skew-normals
    Daniele Durante
    Received: date / Accepted: date
================================================================




Regression models for dichotomous data are ubiquitous in statistics. Besides being useful for inference on binary responses, these methods serve also as  building blocks in more complex formulations, such as density regression,  nonparametric classification and graphical models. Within the Bayesian framework, inference proceeds by updating the priors for the coefficients, typically set to be Gaussians,  with the likelihood induced by probit or logit regressions for the 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,  Markov Chain Monte Carlo strategies face mixing or time-inefficiency issues in large p and small n studies, whereas approximate routines fail to capture the skewness typically observed in the posterior. This article proves that the posterior  distribution for the probit coefficients has a unified skew-normal kernel, under Gaussian priors. Such a novel result allows efficient 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 notable issues. These  advances are  outlined in a genetic study, and further motivate the development of a wider class of conjugate  priors for probit models along with methods to obtain independent and identically distributed samples from the unified skew-normal posterior.



Bayesian inference; Binary data; Conjugacy; Probit regression; Unified skew-normal.




§ INTRODUCTION


There is a relevant interest in several fields towards learning how the probability mass function of a binary response y ∈{0;1} varies with a set of observed predictors x=(x_1, …, x_p)^∈^p  <cit.>. To address this goal, common formulations assume y is a Bernoulli variable whose  probability parameter changes with a linear combination of the predictors under a probit or logit mapping. In the first case (y=1 | x,β)=Φ(x^β), whereas in the second  (y=1 | x,β)={1+exp(-x^β) }^-1, and the goal is to provide inference on β=(β_1, …, β_p)^∈^p. 

Although frequentist inference for the above class of models is well-established <cit.>, the Bayesian approach has attracted an increasing interest since it allows borrowing information, uncertainty quantification, shrinkage and tractable inference via the posterior distribution for the regression coefficients <cit.>. Besides this, predictor-dependent models for binary data are also useful building blocks in more complex Bayesian formulations, such as density regression models <cit.>, additive trees <cit.>,  nonparametric classification <cit.>, graphical models <cit.>,    and others. Although these methods provide popular learning procedures,  there are still computational barriers. Indeed, unlike Bayesian regression with Gaussian data, there are no results on the availability of tractable posteriors for β in regression models for Bernoulli data, under the common Gaussian priors on the coefficients  <cit.>. 


Motivated by the above issue, several computational methods have been proposed for Bayesian regression with binary response data. Popular routines consider data augmentation strategies relying on hierarchical representations which provide conjugate full conditionals, within a Markov Chain Monte Carlo <cit.>. Although being routinely implemented, these methods face poor convergence and mixing in practice, especially for imbalanced datasets <cit.>. A solution is to consider alternative strategies, including carefully tuned or adaptive Metropolis–Hastings <cit.>, as well as more recent generalizations of Hamiltonian Monte Carlo, such as the no u-turn sampler by <cit.>. Both  procedures guarantee computational advantages in large n and moderate p settings compared to data augmentation Markov Chain Monte Carlo. However, when p is large, Metropolis–Hastings has difficulties in exploring the parametric space, whereas Hamiltonian Monte Carlo tends to be expensive  <cit.>. Laplace approximations, variational Bayes and expectation propagation  scale-up computations, but commonly provide Gaussian approximations which affect quality of inference  when the posterior is skewed. This is a common situation in regression for binary data <cit.> and, as will be discussed in the rest of this article, a property which is inherent to the posterior. Refer to <cit.>  for a thorough discussion and comparison among these methods with a specific focus on probit regression. 

Although providing state-of-the-art methods in Bayesian regression for binary response data, the aforementioned strategies are still  sub-optimal compared to situations in which the posterior belongs to a known and tractable class of random variables. Indeed, this result could facilitate the calculation of several quantities relevant to posterior inference, without relying on Monte Carlo methods. This article proves that when the focus is on probit regression models under Gaussian priors for the coefficients, the posterior belongs to the  class of unified skew-normal distributions <cit.>. These variables already appeared in probit models to obtain  flexible link functions via skewed  latent data, instead of Gaussian ones <cit.>. However, this is a different contribution compared to the one in the present article.

To the best of the author's knowledge, the above result is not available in the literature, but can provide important advances. In fact, as discussed in this article, the unified skew-normal posterior guarantees closure properties  <cit.> in addition to explicit formulas for the marginal, joint and conditional posteriors along with predictive distributions and marginal likelihoods for model selection. These quantities involve cumulative distribution  functions Φ_n(·) of n-variate Gaussians, and hence can be efficiently evaluated  for small-to-moderate n studies via recent minimax tilting methods <cit.>. The contribution by <cit.> is additionally useful to obtain independent samples from the posterior exploiting a representation of the unified skew-normal via a linear combination of p-variate Gaussians and n-variate truncated Gaussians. Such results are valid for any dimension and provide, per se, key methodological advances which motivate future theoretical studies and facilitate formal understanding of the skewness typically observed in the posterior. However, the associated inference strategies require evaluation of Φ_n(·) or sampling from n-variate truncated Gaussians, and hence are of practical usefulness in small-to-moderate n studies with, typically, few hundreds of units and any, even huge, p. This scenario is the most challenging for current Markov Chain Monte Carlo routines  <cit.> and, in fact, the methods for posterior inference arising from the new unified skew-normal result significantly improve state-of-the-art algorithms in these situations, thus covering also an important computational gap; refer to <ref> for a more detailed discussion on these aspects. 



§ POSTERIOR INFERENCE IN PROBIT REGRESSION VIA UNIFIED SKEW-NORMALS




 §.§ The unified skew-normal distribution


Before deriving the unified skew-normal posterior induced by Gaussian priors for the β coefficients in a probit regression, let us first introduce the unified skew-normal distribution. Recalling  <cit.>, this random variable unifies different generalizations of the multivariate skew-normal  z ∼sn_p(ξ,Ω,α) <cit.> whose density 2 ϕ_p(z-ξ; Ω) Φ{α^ω^-1(z-ξ) } is obtained by modifying the one of a p-variate Gaussian N_p(ξ, Ω) with the  cumulative distribution function of a standard normal evaluated at α^ω^-1(z-ξ), where ω is a p × p diagonal matrix containing the square root of the diagonal elements in Ω. This strategy  introduces skewness in N_p(ξ, Ω) controlled by  α=(α_1, …, α_p )^∈^p, with ξ=(ξ_1, …, ξ_p)^∈^p and Ω driving location and variability, respectively <cit.>. Indeed, when α=0_p the multivariate skew-normal coincides with N_p(ξ, Ω), whereas, setting p=1, leads to a univariate skew-normal sn(ξ,ω^2,α)  <cit.>. 

Motivated by the success of the above formulation in different studies <cit.>, several extensions have been proposed to incorporate additional properties. Two important generalizations are obtained by adding an additional parameter γ in Φ{α^ω^-1(z-ξ) } to develop the multivariate extended skew-normal <cit.>, and by allowing the skewness-inducing mechanism to be multivariate, thus providing the closed skew-normal family <cit.> which includes a skewness matrix Δ∈^p × n and an n × n full-rank scale Γ in Φ_n(·). Besides increasing flexibility, these extensions allow closure properties for marginals, conditionals and joint distributions, thus providing a general class. <cit.> unify the above generalizations within a single and tractable unified skew-normal representation, obtaining the following density function

    ϕ_p(z -ξ; Ω) Φ_n{γ+Δ^Ω̅^-1ω^-1(z-ξ); Γ-Δ^Ω̅^-1Δ}/Φ_n(γ;Γ),

for z ∼_p,n(ξ,Ω,Δ,γ,Γ). In (<ref>), ϕ_p(z -ξ; Ω) denotes the density of a p-variate Gaussian with expectation ξ=(ξ_1, …, ξ_p)^∈^p,  and p× p variance-covariance matrix Ω=ωΩ̅ω obtained via the quadratic combination between a correlation matrix Ω̅ and a diagonal matrix ω containing the square root of the diagonal elements in Ω. The quantities Φ_n{γ+Δ^Ω̅^-1ω^-1(z-ξ); Γ-Δ^Ω̅^-1Δ} and Φ_n(γ;Γ) denote instead the cumulative distribution functions of the multivariate Gaussians N_n(0_n, Γ-Δ^Ω̅^-1Δ) and N_n(0_n, Γ), evaluated at γ+Δ^Ω̅^-1ω^-1(z-ξ) and γ, respectively,  with the p × n matrix Δ having the main effect on skewness. In fact, when  Δ is zero, (<ref>) coincides with the density of a N_p(ξ, Ω). The vector γ∈^n adds additional flexibility in departures from normality, consistent with the  multivariate extended skew-normal. Refer to <cit.> and <cit.> for details. 

It shall be noticed that <cit.> add a further condition which restricts the (n+p) × (n+p) matrix Ω^*, having blocks Ω_[11]^*=Γ, Ω_[22]^*=Ω̅ and Ω_[21]^*=Ω_[12]^*=Δ, to be a full-rank correlation matrix. As will be clarified in <ref>, this identifiability restriction is not required in the Bayesian setting. In fact, the parameters of the unified skew-normal posterior for the coefficients β are functions of the observed data and the pre-specified hyperparameters of the Gaussian prior, thus avoiding identifiability issues. Nonetheless, such a parameterization will be maintained in the  article to inherit the classical results of the unified skew-normal distribution and to ensure identifiability of the prior, when the findings for the Gaussian case will be generalized to the entire class of  unified skew-normals priors. Sections <ref>–<ref> prove that the posterior for the β coefficients in a probit model with Gaussian priors is a unified skew-normal and study the consequences of this novel finding in posterior inference. 



 §.§ Unified skew-normal posterior for Bayesian probit regression with Gaussian priors


To introduce the reader to the general case consisting of n observations from a probit model with Gaussian prior π(β)=  ϕ_p(β-ξ; Ω), let us first consider a simple setting with a single data point y and one covariate x, such that (y | x, β) ∼{Φ(x β)} and β∼ N(0,1). Although this scenario is uncommon in practice, it provides key intuitions on the role of x ∈ and y ∈{0;1} in driving departures from normality in the posterior distribution.  Indeed, consistent with Lemma <ref>, (β| y, x) has a unified skew-normal posterior when π(β)=ϕ(β). See Appendix A for proofs.


Let (y | x, β) ∼{Φ(x β)} and set π(β)=ϕ(β;1)=ϕ(β), then (β| y, x) ∼_1,1{0,1, (2y-1)x(x^2+1)^-1/2,0,1}, for every x ∈ and y ∈{0;1}. 






Figure <ref> provides the density function of the unified skew-normal posterior for β in the illustrative example, under different combinations of x and y. As expected, (2y-1)x(x^2+1)^-1/2 controls  skewness. Indeed, the higher |x| the more skewness is observed in the posterior. This skewness is either positive or negative depending on the sign of (2y-1)x.   To better understand this result, note that the unified skew-normal in Lemma <ref>  coincides with a basic sn{0,1,(2y-1)x}.


The above results apply more generally to independent response data y_1, …, y_n from a probit model (y_i | x_i, β) ∼{Φ(x_i^β) }, for i=1, …, n, where x_i=(x_i1, …, x_ip)^∈^p denotes the vector of covariates for unit i and β=(β_1, …, β_p)^∈^p the associated coefficients.  Indeed, based on Theorem <ref>, when β has a Gaussian prior β∼ N_p(ξ,Ω) with mean ξ∈^p and full-rank variance-covariance matrix Ω=ωΩ̅ω, the posterior coincides with a unified skew-normal. 


If y=(y_1, …, y_n)^ comprises conditionally independent binary response data from a probit model (y_i | x_i, β) ∼{Φ(x_i^β)}, for i=1, …, n, and β∼ N_p(ξ,Ω), then 

    (β| y, X) ∼_p,n(ξ_,Ω_,Δ_ ,γ_, Γ_),

with the posterior parameters defined as a function of the data and the prior parameters via 
    ξ_= ξ,   Ω_= Ω,   Δ_= Ω̅ω D^s^-1,   γ_= s^-1Dξ,   Γ_= s^-1(D Ω D^+_n)s^-1,
 for every n× p data matrix D=(2y_1-1, …, 2y_n-1)X and any n × n positive diagonal matrix s={(d_1^Ω d_1+1)^1/2, …, (d_n^Ω d_n+1)^1/2}. The generic vector d_i^ denotes instead the i-th row of D, whereas X is the design matrix and I_n the n × n identity matrix.



Adapting equation (<ref>) to the results in Theorem <ref>, it can be immediately noticed, after minor mathematical simplifications, that the density function of the unified skew-normal posterior is 

    π(β| y,X)=ϕ_p(β - ξ; Ω) Φ_n{s^-1Dβ;(s s^)^-1}/Φ_n{s^-1Dξ; s^-1(D Ω D^+ I_n )s^-1},

where Φ_n{s^-1Dξ; s^-1(D Ω D^ + I_n)s^-1} defines the normalizing constant. To clarify  the role of the prior parameters ξ and Ω, along with that of the data y and X, let us consider a constructive representation of the posterior. In particular, adapting known results  from   unified skew-normals  <cit.> to the specific posterior in Theorem <ref>, it can be shown that (β| y, X) has the stochastic representation in Corollary <ref>.


If  (β| y, X) has the unified skew-normal distribution from Theorem <ref>, then

    (β| y, X)=ξ+ω{V_0+ Ω̅ω D^(D Ω D^+I_n)^-1sV_1},      (V_0 ⊥ V_1),

with V_0 ∼ N_p{0_p,Ω̅-Ω̅ω D^(D Ω D^+I_n)^-1D ωΩ̅} and V_1 from a zero mean n-variate truncated normal with covariance matrix s^-1(D Ω D^+ I_n)s^-1 and truncation below - s^-1Dξ.  



Based on (<ref>), ξ has a main role on the location, but has also an effect in controlling departures from normality  since it appears in the truncation s^-1Dξ. Instead, the prior variance-covariance matrix Ω mainly affects scale, via ω, and posterior dependence among the β parameters, while contributing also to the weight matrix assigned to the multivariate truncated Gaussian V_1, along with its variability. Finally, the data in D play a major role in controlling departures from normality. Indeed, if D has elements close to 0, the multivariate truncated Gaussian V_1 has a negligible importance compared to the multivariate Gaussian V_0 in  (<ref>).



 §.§ Inference, prediction and variable selection under unified skew-normal posteriors


A general primary focus in Bayesian regression studies is on marginal posteriors (β_j | y, X), for j=1, …, p, their associated moments and more complex functionals including measures of posterior dependence along with  credible intervals or regions. A fundamental property of unified skew-normals, which potentially facilitates this type of inference, is that such a class of  variables is closed under marginalization, linear combinations and conditioning <cit.>. In particular, adapting the derivations in  <cit.> to  Theorem <ref>, each marginal posterior  is still from a unified skew-normal for every  β_j, j=1, …, p. More specifically, (β_j | y, X)∼_1,n(ξ_ j,Ω_ jj,Δ_ j ,γ_, Γ_) with Δ_ j denoting the j-th row of Ω̅ω D^s^-1,  ξ_ j the j-th element of the prior mean vector ξ, Ω_ jj the entry [jj] in Ω, whereas γ_ and Γ_ coincide with those already defined in Theorem <ref>. A similar result holds also for sub-vectors of coefficients (β_𝒥| y, X), with 𝒥⊂{1; …; p},  linear combinations (a+A^β|  y, X), and conditional posteriors (β_𝒥| y, X,β_𝒥^*), with 𝒥⊂{1; …; p}, 𝒥^* ⊂{1; …; p}, and 𝒥∩𝒥^*=∅. Refer to <cit.> for details to obtain the parameters of these unified skew-normals  from simple transformations of those in  Theorem <ref>. Note that some linear combinations of β, such as x^β, are of particular interest. 

The aforementioned results facilitate graphical representation of marginal or joint posteriors, along with calculation of posterior moments and credible intervals for the probit coefficients via one-dimensional integrals involving marginal posterior densities. This can be done via numerical integration <cit.> whenever it is possible to evaluate Φ_n(·) with efficiency and accuracy. When the focus is on posterior moments, another solution is to obtain such quantities via direct derivation of the moment generating function. Indeed, adapting the result in <cit.>  to   (<ref>) a similar strategy can be considered when studying the functionals of the unified skew-normal posterior, provided that (β| y, X) has moment generating function

    M(t)=exp(ξ^t+0·5 t^Ω t)Φ_n{ s^-1Dξ +s^-1DΩ t;s^-1(D Ω D^+I_n)s^-1}/Φ_n{s^-1Dξ  ; s^-1(D Ω D^+I_n)s^-1},    (t ∈^p).

Exploiting (<ref>) and adapting the derivations in <cit.> to the unified skew-normal in Theorem <ref>, the posterior expectation of β can be explicitly calculated as

    E(β| y, X)=ξ+1/Φ_n{s^-1Dξ; s^-1(D Ω D^+ I_n )s^-1}Ω D^s^-1η,

where η represents an n × 1 vector whose generic i-th component is equal to ϕ(γ̅_i )Φ_n-1(γ̅_-i-Γ̅_-iγ̅_i, Γ̅_-i,-i-Γ̅_-iΓ̅_-i^), with γ̅_̅i̅ and γ̅_-i denoting the i-th element of s^-1Dξ=γ_ and the (n-1) × 1 vector obtained by removing the i-th entry in γ_, respectively. Similarly, Γ̅_-i,-i defines the sub-matrix of s^-1(D Ω D^+ I_n )s^-1=Γ_ without the i-th row and column, whereas Γ̅_-i is the i-th column of Γ_ with the i-th row element removed. Computing the expectation via (<ref>) is more efficient than numerical integration since it requires calculation of n+1 cumulative distribution functions, which is typically much less than the number of evaluations of Φ_n(·) required in numerical integration of marginal posteriors. However, as noticed in  <cit.>, obtaining expressions for higher-order marginal and joint moments via direct derivation of (<ref>) requires tedious  calculations, thus motivating Monte Carlo methods based on samples from the posterior, as discussed in  <ref>. Refer also to  <cit.> and <cit.> for an expression of the  variance-covariance matrix and the cumulative distribution function of a generic unified skew-normal. Both quantities, appropriately computed under the parameters in Theorem <ref>, are useful in posterior inference, especially for credible intervals or regions.

Although inference on the posterior distribution of β is often of interest,  prediction of a future response y_∈{0;1} given the associated covariates x_∈^p and the current data (y, X) is a primary goal in applications of probit models to classification. Within the Bayesian framework, this task requires the derivation of the posterior predictive distribution (y_| y, X, x_), which is simply a Bernoulli having parameter (y_=1 | y, X, x_)=∫Φ(x_^β)π(β| y, X) β, in the binary case. According to Corollary <ref>, this probability parameter is available in explicit form.


If (y_i | x_i, β) ∼{Φ(x_i^β)}, for i=1, …, n, and β∼ N_p(ξ,Ω), then

    (y_=1 | y, X, x_)   =   Φ_n+1{s_^-1 D_ξ ; s_^-1(D_Ω D_^+I_n+1) s_^-1}/Φ_n{s^-1 Dξ ; s^-1(D Ω D^+I_n) s^-1},

with D_ representing the (n+1) × p matrix obtained by adding a last row d_^=x_^ to D, whereas s_={(d_1^Ω d_1+1)^1/2, …, (d_n^Ω d_n+1)^1/2, (d_^Ω d_+1)^1/2}. 



An advantage of (<ref>), compared  to Markov Chain Monte Carlo strategies  <cit.>, is that prediction does not require Monte Carlo integration for ∫Φ(x_^β)π(β| y, X) β via sampling of β from the posterior, and hence the computational burden does not depend on p.  As will be discussed in <ref>, this result is especially useful in large p and small-to-moderate n studies. 

The above derivations  are further helpful in obtaining explicit methods to perform Bayesian selection among models ℳ_1, …, ℳ_K characterizing, in general, different subsets 𝒥_1, …, 𝒥_K of covariates  entering the linear predictor. Although there are different strategies for model selection <cit.>, the general  approach formally defines prior probabilities (ℳ_1), …, (ℳ_K) for the set of models, and subsequently ranks them via the posterior probabilities (ℳ_k | y, X) ∝(ℳ_k)∫(y |ℳ_k,X, β_𝒥_k)π( β_𝒥_k|ℳ_k) β_𝒥_k, k=1, …, K <cit.>. Clearly, the major issue in this task is the calculation of ∫(y |ℳ_k, X, β_𝒥_k)π(β_𝒥_k|ℳ_k)β_𝒥_k which may be intractable in the absence of conjugacy, thus requiring Monte Carlo integration or  approximations <cit.>. This procedure can be implemented  also in probit models leveraging the methods  in <ref>, but inference and computational performance face the issues previously discussed. Corollary <ref> provides instead an explicit formula for the marginal likelihood in probit models with Gaussian priors, which can be easily evaluated, especially in large p and small-to-moderate n settings of interest in such studies.



Let ℳ_k denote the probit regression model for y_1, …, y_n including only the covariates with indices in the subset 𝒥_k ⊂{1;…;p}  and assume (β_𝒥_k|ℳ_k) ∼ N_p_k(ξ_k,Ω_k), with p_k=|𝒥_k| and  β_𝒥_k∈^p_k the probit coefficients for the covariates in model ℳ_k, then

    ∫(y |ℳ_k,X, β_𝒥_k)π( β_𝒥_k|ℳ_k) β_𝒥_k=Φ_n{s_k^-1D_kξ_k;s_k^-1(D_k Ω_k D_k^+I_n) s_k^-1},

for every model ℳ_k, k=1, …, K, where D_k= (2y_1-1, …, 2y_n-1)X_k ∈^n × p_k, s_k={(d_1k^Ω_k d_1k+1)^1/2, …, (d^_nkΩ_k d_nk+1)^1/2}∈^n × n_+, and X_k ∈^n × p_k denoting the n × p_k  design matrix of covariates with indices in 𝒥_k.



Equation (<ref>) is additionally useful to compute Bayes factors <cit.> and to perform Bayesian model averaging <cit.> without sampling from the posterior. 



 §.§ Computational considerations and sampling procedures




All the inference methods outlined in <ref> can, in principle, proceed via direct strategies without sampling from the posterior, thus improving the available procedures in large p applications. The only  barrier, which is relevant for  a large n, is evaluation of Φ_n(·). Quasi-randomized Monte Carlo <cit.> allows, in fact, accurate calculation of Φ_n(·) for small n, and have been recently improved via minimax tilting  <cit.> to ensure effective evaluation of Φ_n(·)  in moderate n studies. This procedure, available in the  library , has a rare vanishing asymptotic relative error, thus allowing tractable inference without sampling from the posterior in studies having, typically, few hundreds of units. This strategy is also useful in larger n applications when few evaluations of Φ_n(·) are required, as in prediction of not many outcomes and in selection among few models. However, for general inferential tasks requiring a plenty of evaluations of  Φ_n(·), such as in numerical integration, moments calculation and high-dimensional prediction or model selection, inference without sampling from the posterior might face non-negligible increments in computational time when n is large; refer to <cit.> for details on scalability in the evaluation of Φ_n(·). In this situation, sampling from the posterior provides a tractable and common strategy to obtain numerical evaluations of generic functionals via Monte Carlo integration approximating E{g(β) | y, X}=∫ g(β)π(β| y, X) β. Indeed, the availability of a large number R of  replicates from the unified skew-normal posterior, allows fast and accurate approximation of E{g(β) | y, X} via ∑_r=1^ g(β^(r))/R. 

Popular routines addressing the above goal require data augmentation Markov Chain Monte Carlo <cit.>, which provide poor performance, especially in imbalanced high-dimensional studies <cit.>. This issue can be addressed via   Algorithm <ref>, which combines the stochastic representation of the unified skew-normal posterior in Corollary <ref> with a new scheme proposed by <cit.> to obtain independent samples from  multivariate truncated Gaussians. This routine relies on minimax tilting  and accept-reject methods to improve the acceptance rate of classical rejection sampling, while  avoiding convergence and mixing issues of Markov Chain Monte Carlo methods. By combining this sampler with classical routines  for multivariate Gaussians, Algorithm <ref> inherits these properties, thus improving the computational methods discussed in <ref>, especially in large p  and small-to-moderate n applications. Clearly, when n increases and p decreases, sampling from the n-variate truncated Gaussian progressively affects computational time in favor of  more efficient Markov Chain Monte Carlo strategies which directly explore the p-dimensional parametric space <cit.>. In this situation, a possibility to scale-up the computations is to exploit the  structure of Algorithm <ref> to perform parallel computing.  Another alternative is to leverage the closure properties of multivariate truncated Gaussians under conditioning <cit.> and iteratively block-update sub-vectors of V_1 whose dimension still allows efficient sampling via <cit.>. Although this hybrid strategy could induce some auto-correlation in the posterior samples of β, the blocking approach typically guarantees  improvements in mixing and convergence <cit.>.

It is also worth noticing that <cit.> applied his accept-reject method to Bayesian probit regression. However, unlike  Algorithm <ref>, the proposed strategy  requires sampling from (n +p)-variate truncated Gaussians. Separating these two blocks, as in Algorithm <ref>, reduces computational complexity and allows parallel computing. A more similar representation can be found in <cit.> and  in the documentation of the  library  by <cit.>. In fact, the resulting routines are closely related to Algorithm <ref>. However, <cit.> and <cit.> base their derivations on different arguments without noticing that the posterior is indeed a unified skew-normal. This last result and its broader implications are arguably the most important contribution of the present article.





 

Finally, Algorithm <ref> can be also adapted to sample from a generic unified skew-normal. This can be broadly useful much beyond Bayesian  inference. An example is parametric bootstrap <cit.> for frequentist inference on  the unified skew-normal parameters. 




 §.§ A  class of conjugate unified skew-normal priors for Bayesian probit regression


The derivations in <ref> suggest the more general result outlined in Corollary <ref>, thereby allowing tractable inference in  Bayesian probit regression under more flexible priors for β.

If (y_i | x_i, β) ∼{Φ(x_i^β)} independently for i=1, …, n, and β is assigned a _p,m(ξ,Ω,Δ,γ,Γ) prior <cit.>, then

    (β| y, X) ∼_p,m+n(ξ_,Ω_,Δ_ ,γ_, Γ_),

with posterior parameters ξ_=ξ, Ω_=Ω, Δ_=(Δ  Ω̅ω D^s^-1), γ_=(γ^  ξ^ D^s^-1)^ and Γ_ characterizing an (m+n) × (m+n) full-rank correlation matrix having block entries Γ_[11]=Γ, Γ_[22]=s^-1(D Ω D^+I_n) s^-1, Γ_[21]=Γ_[12]^=s^-1DωΔ.





According to Corollary <ref>, tractable inference  in Bayesian probit regression is possible under a broader class of priors. Indeed, all the methods in <ref>–<ref> also apply to  this more general case, since  the posterior in  (<ref>) is still a unified skew-normal. This ensures increased flexibility in prior specification, thus allowing departures from normality. Although the general unified skew-normal choice may be  uncommon in applied contexts, it shall be noticed that this class incorporates several  priors of interest, including multivariate Gaussians, independent skew-normals for each β_1, …, β_p, and multivariate skew-normals for β <cit.>.




§ EMPIRICAL STUDIES


To evaluate the methods developed in <ref> and compare performance with the popular strategies for Bayesian inference in probit regression discussed in <ref>, let us consider an online available dataset on the gene expression of n=74 normal and cancerous biological tissues at p-1=516 different tags <cit.>. An overarching focus in these applications is quantifying the effects of gene expression on the probability of a cancerous tissue and predicting the status of new tissues as a function of the gene expression <cit.>. Consistent with this goal, let us focus on studying the location of the posterior for β  and the predictive distribution in the Bayesian model (y_i | x_i, β) ∼{Φ(x_i^β)}, for i=1, …, n, with β∼ N_517( 0_517, 16· I_517) prior. In this probit regression, x_i denotes the vector having an intercept term and the gene expressions for tissue i, whereas y_i is either 1 or 0 if the tissue is cancerous or not, respectively. 

The choice of a weakly informative prior for the β coefficients is motivated by the guidelines in <cit.>  and by similar implementations from <cit.> and <cit.>.   In line with these contributions, the gene expressions at the  516 different tags have been also standardized to have mean 0 and standard deviation 0·5. To assess predictive performance, the prior for β is updated with the information of 50 randomly chosen observations, and out-of-sample classification via the posterior predictive distribution is made on the 24 held-out units.






Although other datasets could be considered, it shall be emphasized that state-of-the-art computational methods for probit regression provide valuable strategies in a variety of applications, but face  mixing and time-inefficiency issues in large p and small n studies  <cit.>. As shown in Figs. <ref>–<ref> and in Table <ref>, the novel results outlined in <ref>  allow notable improvements in these large p and small n studies, thus providing straightforward Bayesian inference in relevant applications where this task was previously impractical. To clarify these results, the strategies in  <ref> are compared with state-of-the-art procedures, covering the data augmentation  Gibbs sampler  by <cit.>, the Hamiltonian no u-turn sampler  in <cit.> and the adaptive  Metropolis–Hastings from <cit.>. To increase acceptance rate and efficiency, the starting Gaussian proposal for the Metropolis–Hastings routine has been initialized with the mean and the rescaled variance-covariance matrix provided by an expectation propagation approximation. Consistent with  <cit.> and <cit.>, the scaling factor has been set to  2·38^2/p. 






The above Markov Chain Monte Carlo routines were run for 20000 iterations after a burn-in of 5000, and can been easily implemented in , leveraging the libraries  ,  , and a combination of  and , respectively. Although certain routines  converged rapidly than others and with excellent mixing, the same  settings were considered for all the algorithms  to  facilitate comparison. The  sampling scheme proposed in Algorithm <ref> provides instead independent samples from the exact posterior and hence requires no burn-in or convergence checks. Refer to the Supplementary Materials for details on code and implementation. 




According to Table <ref>,  the Hamiltonian no u-turn sampler has the same mixing of Algorithm <ref> which displays, however, a significantly faster sampling speed. This could be due to the number of leap-frog steps required at each iteration of the no u-turn sampler  <cit.>. As expected, the data augmentation Gibbs sampler and Metropolis–Hastings display lower mixing, but provide similar or improved running time compared to Hamiltonian no u-turn samplers. However, as is clear from Figs. <ref>–<ref>,  this reduction in mixing has a direct effect on the accuracy of posterior inference and prediction. There is instead  an almost perfect match between Monte Carlo and direct estimates of posterior means and posterior predictive probabilities for the proposed Algorithm <ref> and the Hamiltonian no u-turn sampler. However, as already discussed, such a routine is significantly slower than Algorithm <ref> in this application. These computational gaps further increased when focusing on larger p studies, with the competing methods becoming rapidly impractical. Conversely, the inference and sampling methods relying on the unified skew-normal results have difficulties in scaling with n. This claim is confirmed by a voice rehabilitation study presented at the online repository in the Supplementary Materials. However,  also in this application having doubled n and almost halved p, Algorithm <ref> remains still competitive. 




§ FINAL CONSIDERATIONS AND FUTURE DIRECTIONS OF RESEARCH


This article shows that the posterior for the coefficients in a probit regression having Gaussian priors is a unified skew-normal  <cit.>, thus allowing key advances in  Bayesian modeling of binary response data, especially for large p and small-to-moderate n studies. Indeed, unified skew-normals have moment generating functions involving known quantities, tractable additive representations, and are closed under marginalization, conditioning and linear transformations, thus facilitating derivations of marginal likelihoods for model selection and posterior predictive distributions. As shown in the empirical assessments in <ref>, in small-to-moderate n settings with large or even huge p, posterior inference can proceed via direct methods or via an efficient sampler from the posterior, which notably improves available computational methods.

The above results could lead to computational gains also in more complex formulations relying on predictor-dependent observed or latent binary data, such as in mixture models for density regression <cit.>. For instance, leveraging results  in <ref>, binary classification via Gaussian processes <cit.> could avoid sampling or approximations by exploiting closure properties of unified skew-normals, especially under conditioning. Moreover, when binary regression serves as a latent dictionary function, sampling the binary data via (<ref>), instead of conditioning on β, could speed-up computations. Finally, the  novel conjugacy results in <ref>, open new avenues for  incorporation of skewness in prior specification.

There are also different directions for future advances. For instance, improved studies on the moment generating function of the unified skew-normal could facilitate direct calculation of relevant functionals without the need to sample from (β| y, X).  On the same line, improving the methods for efficient evaluation of Φ_n(·) in large n applications, either via data transformations, blocking methods <cit.> or recent algorithms <cit.>, could enlarge the range of applications which allow direct inference, prediction and model selection, without sampling from (β| y, X). Also approximations of the exact posterior, which preserve the skewness but allow analytical inference provide an interesting direction. Finally, more detailed studies on particular forms of  the prior variance-covariance matrix Ω, such as those associated with g-priors or generalized g-priors <cit.> and limiting cases arising from flat priors, could provide novel insights on the effect of these specific choices and, potentially, lead to simplifications in the unified skew-normal parameters which may ease posterior inference. It is also possible to consider hyper-priors for (ξ, Ω) such as the normal-inverse-Wishart. With this choice, Theorem <ref> holds only for the full conditional  (β| y, X, ξ, Ω). Moreover, it is straightforward to notice that (ξ, Ω|  y, X, β) has still a normal-inverse-Wishart kernel, since (ξ, Ω) only enter the Gaussian prior for β. Hence, although an hierarchical prior would not allow direct sampling from a unified skew-normal posterior, Gibbs sampling methods can be easily applied to this situation.

As discussed in <ref>, the availability of an exact posterior with tractable stochastic representations and closure properties can also motivate novel finite-sample and asymptotic theory. Finally, although the studies in <ref> and the discussion in <ref> provide the general guidelines on the practical usefulness of the unified-skew results in <ref>, additional quantitative assessments on scalability and its relations with specific prior settings or dataset structures, are certainly interesting.



§ ACKNOWLEDGEMENT

The author is grateful to Adelchi Azzalini, David Dunson and Giacomo Zanella for the stimulating discussions on this topic and the comments received on a first version of this article. The Editor, the Associate Editor and the referees are also acknowledged for their useful suggestions. This article was partially supported by the grant number 670337–INDIMACRO.




§ SUPPLEMENTARY MATERIALS

Codes and tutorial implementations for the methods associated with this article are available at .







§ APPENDIX A


Let Φ(x β)^y{1-Φ(xβ)}^1-y=Φ{(2y-1)x β} denote the probability mass function of y in Lemma <ref>. Direct application of the Bayes rule provides 

    π(β| y, x) ∝ϕ(β)Φ{(2y-1)x β}=  ϕ(β)Φ{(2y-1)x(x^2+1)^-1/2β; (x^2+1)^-1}.
 
Hence, letting ξ_=0, Ω_=1, Δ_=(2y-1)x(x^2+1)^-1/2, γ_=0,  and Γ_=(x^2+1)^-1+Δ_^Δ_=1, provides the kernel of the unified skew-normal in Lemma <ref>, with  correlation matrix Ω_^* having block entries Ω_[11]^*=Ω_[22]^*=1 and Ω_[21]^*=Ω_[12]^*=Δ_. 



Adapting the proof of Lemma <ref>, it is possible to write the joint probability mass function of the responses y  as ∏_i=1^nΦ{(2y_i-1)x_i^β}=Φ_n(Dβ; I_n)=Φ_n{s^-1 Dβ;(ss^)^-1}, with D and s defined as in Theorem <ref>.  Combining this likelihood for y with the Gaussian prior for β provides 

    π(β| y, X) ∝ϕ_p(β-ξ; Ω)Φ_n{s^-1 Dβ; (ss^)^-1}=ϕ_p(β-ξ; Ω)Φ_n{s^-1 Dξ +s^-1 D( β-ξ); (ss^)^-1}.

To establish the relation among the above kernel and  the unified skew-normal density in (<ref>), note that ( β-ξ)=ωΩ̅Ω̅^-1ω^-1( β-ξ). Therefore, letting ξ_=ξ, Ω_=Ω, Δ_=Ω̅ω D^ s^-1, γ_=s^-1 Dξ  and Γ_=(ss^)^-1+Δ_^Ω̅^-1Δ_=s^-1s^-1+s^-1 DωΩ̅Ω̅^-1Ω̅ω D^ s^-1=s^-1(D Ω D^+I_n) s^-1, provides the kernel of a unified skew-normal whose parameters coincide with those presented in Theorem  <ref>. To conclude the proof it is also necessary to guarantee that 

    Ω_^*=[ s^-1(D Ω D^+I_n) s^-1             s^-1 DωΩ̅;           Ω̅ω D^ s^-1                    Ω̅ ]=[ s^-1    0;    0 ω^-1 ]×[ D Ω D^+I_n         DΩ;        ΩD^          Ω ]×[ s^-1    0;    0 ω^-1 ]

is a full-rank correlation matrix. This last result can be easily proved by noticing that Ω_^* coincides with the correlation matrix of the random vector (z_1^, z_2^)^ where z_1=D z_2+ ϵ with E(ϵ)=0_n, E(ϵϵ^)=I_n, and z_2 is a p-variate random variable having zero mean and positive definite variance-covariance matrix E(z_2 z_2^)=Ω=ωΩ̅ω. Finally,  s={(d^_1Ω d_1+1)^1/2, …, (d^_nΩ d_n+1)^1/2} is the diagonal matrix with the square root of the diagonal elements of E(z_1z_1^)=D Ω D^+I_n.



The proof is a simple adaptation of equation (7.4) in <cit.> to the unified skew-normal posterior in Theorem <ref>. In particular,  according to  <cit.> the posterior in equation (<ref>) has the same distribution of the random variable

    ξ_+ω_(V_0+Δ_Γ^-1_V_1)

with V_0 ∼ N_p(0_p, Ω̅_-Δ_Γ^-1_Δ_^) and V_1 from an n-variate truncated normal with mean 0_n, covariance matrix Γ_ and truncation below - γ_. Substituting the posterior parameters in this stochastic representation with their expressions in Theorem  <ref> concludes the proof. To clarify this final claim, note that Δ_Γ^-1_= Ω̅ω D^s^-1s(D Ω D^+ I_n)^-1s=Ω̅ω D^(D Ω D^+ I_n)^-1s and that Δ_Γ^-1_Δ_^ coincides with the matrix Ω̅ω D^s^-1s(D Ω D^+ I_n)^-1ss^-1D ωΩ̅=Ω̅ω D^(D Ω D^+ I_n)^-1D ωΩ̅.
 


Recalling the expression for π(β| y, X) outlined in equation (<ref>), the posterior predictive probability (y_=1 | y, X, x_)=∫Φ(x_^β)π(β| y, X) β can be  expressed as 

    Φ_n{s^-1Dξ ; s^-1(D Ω D^+I_n) s^-1}^-1∫ϕ_p(β -ξ; Ω)Φ(x_^β)Φ_n{s^-1 Dβ; (ss^)^-1}β.

Exploiting the proof of Theorem <ref>, the quantity inside the above integral can be re-expressed as ϕ_p(β -ξ; Ω)Φ_n+1{s_^-1 D_β; (s_s_^)^-1}, with D_ and s_ defined as in Corollary <ref>. Comparing, now, this function with the density of the unified skew-normal posterior in  (<ref>) it can be immediately noticed that ϕ_p(β -ξ; Ω)Φ_n+1{s_^-1 D_β; (s_s_^)^-1}  is the kernel of a unified skew-normal with normalizing constant ∫ϕ_p(β -ξ; Ω)Φ(x_^β)Φ_n{s^-1 Dβ; (ss^)^-1}β=Φ_n+1{ s_^-1 D_ξ ; s_^-1( D_Ω D_^+I_n+1) s_^-1}. To conclude the proof, substitute this quantity in the above formula for the posterior predictive probability (y_=1 | y, X, x_). 



To prove Corollary <ref> simply notice that  ∫(y |ℳ_k, X, β_𝒥_k)π( β_𝒥_k|ℳ_k) β_𝒥_k is the normalizing constant of the posterior for β_𝒥_k in model ℳ_k. Hence, adapting equation (<ref>) to model ℳ_k, leads to ∫(y |ℳ_k, X, β_𝒥_k)π( β_𝒥_k|ℳ_k) β_𝒥_k=Φ_n{s_k^-1D_kξ_k;s_k^-1(D_k Ω_k D_k^+I_n) s_k^-1}.



To prove Corollary <ref> it suffices to generalize Theorem <ref>. In particular, adapting the proof of Theorem <ref> to the case in which β∼_p,m(ξ,Ω,Δ,γ,Γ), provides 

    π(β| y, X) ∝ϕ_p(β-ξ, Ω)Φ_n{s^-1 Dβ; (s s^)^-1}Φ_m{γ+Δ^Ω̅^-1ω^-1(β-ξ); Γ-Δ^Ω̅^-1Δ}.

To proceed with the proof, note that exploiting Theorem <ref>, it is possible to re-write Φ_n{s^-1 Dβ; (s s^)^-1} as Φ_n{s^-1 Dξ +(Ω̅ω D^ s^-1)^Ω̅^-1ω^-1( β-ξ); s^-1(D Ω D^+I_n) s^-1- s^-1 DωΩ̅Ω̅^-1Ω̅ω D^ s^-1}. 

Consistent with this result, let us define the appropriate parameters for which the above kernel coincides with the one of the unified skew-normal posterior in Corollary <ref>. This goal is easily accomplished by setting ξ_=ξ, Ω_=Ω,  Δ_=(Δ  Ω̅ω D^ s^-1), γ_=(γ^  ξ^ D^s^-1)^ and Γ_ a full-rank correlation matrix with blocks Γ_[11]=Γ, Γ_[22]=s^-1(D Ω D^+I_n) s^-1, Γ_[21]=Γ_[12]^=s^-1DωΔ. As in Theorem <ref>, it is also necessary to ensure that 

    [                     Γ           Δ^ω D^ s^-1                    Δ^                      ;              s^-1 DωΔ s^-1(D Ω D^+I_n) s^-1             s^-1 DωΩ̅;                     Δ           Ω̅ω D^ s^-1                    Ω̅ ],

is a full-rank correlation matrix. This result follows under minor modifications of the proof in Theorem <ref>, after noticing that by definition the (m+p) × (m+p) unified skew-normal prior matrix Ω^*  having block entries Ω_[11]^*=Γ, Ω_[22]^*=Ω̅=ω^-1Ωω^-1, Ω_[21]^*=Ω_[12]^*=Δ, is a full-rank correlation matrix.

biometrika


