Binary regression models are essential tools in statistical modeling, particularly for understanding relationships between binary outcomes and predictor variables. Traditional approaches often rely on the logit or probit link functions due to their theoretical appeal and interpretability. However, these standard link functions may fail when the sigmoid function of the underlying data exhibits heavy tails or extreme responses [10, 21, 29], as often encountered in genetics and biomedical research. For example, binary outcomes in genetic studies, such as the presence or absence of a specific disease, are influenced by complex biological processes that standard links may inadequately capture. Moreover, inference results from these models may become unreliable in the presence of outlying covariates or incoherent responses [22, 24, 40]. Addressing these limitations requires exploring alternative approaches that better account for the complexities of such data.

In this context, the t-distribution emerges as an attractive alternative due to its flexibility in modeling heavy-tailed data, controlled by the degrees of freedom parameter v. This study aims to further investigate the t-link model’s potential to enhance robustness and mitigate sensitivity to outliers. Specifically, we explore its application to high-dimensional variable selection, a scenario that has received little, if any, attention in the current literature.

Bayesian variable selection is a theoretically intricate and computationally demanding problem. These challenges stem from the balance required between numerical computation, prior specification, and analytical evaluation, as highlighted in earlier works [18, 12, 35]. High-dimensional variable selection introduces additional computational burdens due to the increasing number of covariates. To address this, [32] proposed the skinny Gibbs sampler using the continuous spike-and-slab prior introduced by [31], demonstrating its consistency in selecting the true model in logistic regression. [38] extended this algorithm to the probit model, further proving its theoretical properties.

Building on these earlier work, we propose the hierarchical skinny Gibbs sampler with t-link (HSGT). The HSGT method extends the skinny Gibbs sampler framework by incorporating the t-distribution in the link function. Furthermore, we enhance the algorithms of [32] and [38] by introducing a hierarchical prior on the variance component of active covariates. This update eliminates the time-consuming and problematic tuning process required in previous approaches. The advantages of our method are twofold: (1) the t-link provides a flexible and robust family of link models, with probit and logit as special cases, enabling a more accurate representation of underlying biological and medical phenomena. Theoretically, the probit link can be recovered as v → ∞, while the logit link can be approximated by a t-link with a specific value of v [1, 32]. And its robustness in estimation and inference under outliers - both in covariates and responses - is well established in the literature [40]. and (2) the skinny Gibbs sampler ensures computational efficiency and scalability to large datasets, a critical requirement in genomics and other high-dimensional fields.

The paper is organized as follows. In Section 2, we introduce the model for variable selection in t-link binary regression and discuss the hierarchical Gibbs algorithms, including the exact Gibbs (HEGT) and the skinny Gibbs (HSGT) methods. The theoretical foundations of the HSGT method, along with consistency results, are detailed in Section 3, with comprehensive proofs provided in Appendix A. Simulation studies evaluating model selection under various scenarios are presented in Section 4, with a particular focus on comparing the performance of the t-link and probit link models, both with and without outliers. In Section 5, we apply the proposed method to a genomic dataset, providing a thorough comparison in prediction accuracy among models. Section 6 concludes with a discussion.

We consider a scenario of binary classification in high-dimensional analysis, specifically focusing on a t-link function to model binary response in medicine or genomics studies where outliers and heavy-tailed distributions are observed. For instance, identifying key genetic variants (SNPs) associated with the presence or absence of a disease (e.g., cancer). Let E i ∈ { 0 , 1 } represent the binary response variable, where E i = 1 indicates that the ith subject has the disease, and E i = 0 otherwise, i = 1 , … , n. The covariate vector x i = ( x i 1 , x i 2 , … , x i p ) ′ ∈ R p contains the expression levels of p genes for the ith subject. The probability of E i = 1 is modeled through the generalized liner regression model framework as (2.1) P ( E i = 1 ∣ x i , β , v ) = Ψ ( x i ′ β , v ) , where Ψ ( · ) is the cumulative distribution function of a distribution, such as the standard normal distribution for the probit link and the standard t-distribution for the t-link. The parameter v may be a scalar or a vector of parameters that is fixed or estimated in the regression model. In the t-link regression, v represents the degrees of freedom. The regression coefficients β is a p × 1 vector, which satisfies ‖ β ‖ 0 < < p, meaning that only a small subset of the components of β are non-zero. This sparsity assumption is especially relevant in high-dimensional contexts where the number of predictors p is much larger than the sample size n.

The observed-data likelihood function of all n observations is given as (2.2) L n β ∣ E , X , v = ∏ i = 1 n Ψ ( x i ′ β , v ) E i 1 − Ψ ( x i ′ β , v ) ( 1 − E i ) . Using the latent variable approach [1] and representing the t-distribution as a scale mixture of Gaussian, the complete-data likelihood function is expressed as: (2.3) L ( β , Z , ω ∣ E , X , v ) = ∏ i = 1 n N ( Z i ∣ x i ′ β , ω i 2 ) IG ( ω i 2 ∣ v / 2 , v / 2 ) × ( E i I ( Z i ≥ 0 ) + ( 1 − E i ) I ( Z i < 0 ) ) , where Z i | ω i 2 is a normal with mean x i ′ β and variance ω i 2 and ω i 2 has an inverse-gamma distribution with the scale and the shape as v / 2. The aim is to efficiently infer the sparsity pattern encoded in the regression coefficients β through a hierarchical Bayesian model with consistency guarantees. Throughout this methodological exposition, we assume that the columns of X are standardized to have zero mean and unit variance.

In this section, we present the asymptotic properties of the hierarchical prior setup of the algorithm on binary regression with the student t-link. Let t = { t 1 , t 2 , … , t t } ⊆ [ p ] represent the true model, indicating that the non-zero locations of the true coefficient vector correspond to t = ( j , j ∈ t ). We assume | t | is fixed. Let β 0 ∈ R p be the true coefficient vector, and β 0 , t ∈ R | t | represent the vector of true non-zero coefficients.

For a given model k with active components indexed by k ⊆ [ p ] with size | k | at most m n , we define ℓ n ( β k ) as the log-likelihood and S n ( β k ) = ∂ ℓ n ( β k ) / ∂ β k as the score function. The restriction on the largest model size is a common practice in Bayesian high-dimensional asymptotic theory to effectively control the posterior ratio and achieve desirable consistency results. The variable selection consistency by the skinny Gibbs has been shown under the logit link [32] and the probit link [38], respectively. Here, we extend the consistency results to a t-link model with a known degrees of freedom.

The Hessian matrix of ℓ n ( β k ) is defined as H n ( β k ) = − ∂ 2 ℓ n ∂ β k ∂ β k ′ = ∑ i = 1 n x i , k x i , k ′ w i ( β k ) = X k ′ W ( β k ) X k , where W ( β k ) = Diag { w i ( β k ) } for i = 1 , … , n with w i ( β k ) as E i ( v + 1 ) v 1 + μ i , k 2 v − 1 μ i , k r i , k + r i , k 2 + ( 1 − E i ) · − ( v + 1 ) v 1 + μ i , k 2 v − 1 μ i , k R i , k + R i , k 2 with r i , k = ψ ( μ i , k , v ) / Ψ ( μ i , k , v ) and R i , k = ψ ( μ i , k , v ) / [ 1 − Ψ ( μ i , k , v ) ], ψ ( μ i , k , v ) and Ψ ( μ i , k , v ) as the density and the cumulative probability functions evaluated at μ i , k of a standard t distribution with v degrees of freedom, and μ i , k = x i , k β i , k .

Before presenting the main results, we introduce the following notation: (i) for any x , y ∈ R, x ∨ y and x ∧ y represent the maximum and minimum of x and y, respectively. For any positive real sequences x n and y n , (ii) x n ≲ y n or x n = O ( y n ) implies there exists C > 0 a constant such that | x n | ≤ C | y n | for all large n. (iii) x n ≪ y n or x n = o ( y n ) if x n / y n → 0 as n → ∞. (iv) If x n ∼ y n , there exists K 1 > K 2 > 0 such that K 2 < y n / x n ≤ x n / y n < K 1 .

Let a = ( a 1 , a 2 , … , a p ) ′ ∈ R p then ‖ a ‖ 2 and ‖ a ‖ 1 is the ℓ 2 and ℓ 1 norm of a respectively and ‖ a ‖ max = max 1 ≤ j ≤ p | a j |. For any real symmetric matrix M, λ max ( M ) and λ min ( M ) are the maximum and minimum eigenvalues of M, respectively. The following conditions are assumed to obtain the asymptotic results:

Condition (C1) (Condition on dimension (n and p)) For some constant 0 < d ′ < 1, let p ≫ n, log p = o ( n ) and m n = O ( ( n / log p ) 1 − d ′ 2 ∧ p ), as n → ∞.Condition 3 guarantees that our proposed method is capable of handling high-dimensional setting analysis where the number of predictors increases at a sub-exponential rate relative to n. The parameter m n restricts the analysis to a set of sufficiently large models. Also, to avoid over-fitting, it is required that m n ≪ n. Similar assumptions regarding model size limitations are frequently encountered in the sparse estimation literature [14, 34, 28].

Condition (C2) (Condition on Design Matrix X) For x i ∈ R p , i = 1 , 2 , … , n. i.

(Boundedness) P ( | | x i | | m a x ≤ C 0 ) = 1 for some C 0 > 0 and P ( · ) is the probability measure.

Condition (C3) (Conditions on True Model) For some constant c 0 > 0, min j ∈ t ‖ β 0 , j ‖ 2 2 ≥ c 0 | t | Λ | t | ( log p n ∨ 1 log p ) , and Λ | t | = max t : | t | ≤ c 0 λ max ( n − 1 X k ′ X k ) for any integer c 0 > 0.This condition (also known as beta-min condition) guarantees true regression coefficients β 0 have a finite number of nonzero entries and a bounded ℓ 1 -norm. It holds if we assume that the true inactive beta components ( β t c ) can be nonzero but ‖ X t c β t c ‖ 2 = O ( 1 ) in line with similar assumption in [33].

Condition (C4) By adapting the line of reasoning(Conditions on prior hyperparameters) n τ 0 2 = O ( 1 ) and for some constants a 1 , a 2 > 0, the hyperparameters satisfy a 1 < r , s < a 2 and q 2 ∼ p − 1 , where r and s are the shape and scale parameters of the inverse gamma prior on τ 1 2 defined in Section 2.1.Condition 3 suggests appropriate conditions for the hyperparameters. Similar assumptions have also been considered in [43], [26] and [9]. We impose a prior on τ 1 2 with its value estimated within the MCMC iterations. This eliminates the need for the tuning procedure required in [32].

Lemma 1 is used to prove Theorem 1 which demonstrates the proposed method guarantees posterior ratio consistency.

We assess the performance of the t-links, the logit link, and the probit link in variable selection, particularly in the presence of outlying observations. In addition, we evaluate the performance under both the exact and the skinny Gibbs algorithms. For each data replicate, ten models are fitted and compared, including the hierarchal skinny Gibbs t-link model(HSGT), the exact Gibbs t-link model (HEGT), the skinny Gibbs logit model (SLogit)[37], the exact Gibbs logit model (ELogit), the skinny Gibbs probit model (SProbit) and the exact Gibbs probit model (EProbit)[38]. For the t-link models, we consider ν = 1 ( Cauchy ) , 3, and 7, respectively, where ν = 1 represents the heaviest tail case among integer values of ν, and ν = 7 is sometimes used as an approximation to the logit link. In addition, three LASSO models with different links are fitted to provide a frequentist benchmark and to further illustrate the robustness of the link functions under alternative estimation methods.

Here we assume ν as known and fix it at a specific value. The proposed algorithm is flexible and can be applied to any value of ν. The “best” ν for a given dataset can be determined ad hoc using standard model selection criteria, such as cross validation, DIC or WAIC. All the theoretical properties are proved for any fixed value of ν.

In the simulation, the first | k | columns of the design matrix X correspond to active covariates with non-zero coefficients, while the remaining columns represent inactive covariates with zero regression coefficients. Each row x i of X is independently generated from a normal distribution with a p-dimensional covariance matrix Σ. The binary response E is generated from a Bernoulli distribution with the probability of success Ψ ( X β , v ), where Ψ ( · ) is the cumulative distribution function of a standard normal for the probit link ( v → ∞), a standard logistic distribution for the logit link, and a standard t distribution with the degrees of freedom v = 1 , 3 , 7. We set | k | = 4, n = 100, β = ( 3 , 1.5 , 1.0 , 0.5 ) ′ . We vary the true link models Ψ − 1 (the probit link, the logit link or the Cauchy link ( v = 1)), the dimension in X ( p = 100 or 500) and its correlation structure (independent or dependent). For independent X, the covariance matrix Σ = I; in the case of correlated covariates, Σ is structured into a two-block matrix corresponding to X k and X k c : with a constant correlation of 0.1 between blocks, and within each block, an AR(1) correlation structure with autocorrelation ρ values of 0.5 across the covariates. Thus, there are a total of 12 scenarios to simulate. A total of 50 data replicates are generated under each scenario.

For each simulated dataset, we also generate two modified datasets identical to the original except for one observation The first dataset is altered to create a “bad leverage point” as described in [22]. This observation is randomly selected from among those with E i = 1 and the value of the first covariate is changed to −10 (note the regression coefficient associated with the first covariate is 3). It is referred to as a “bad leverage point” because it exerts leverage while contradicting the association between the response and explanatory variable observed in the rest of the data. The second dataset creates a “non-leverage outlier” by selecting the observation with the largest magnitude of E i − Ψ ( x i β ) and set the value of binary observation as 1 − E i . A comparison of the fits from the 10 Bayesian models and three LASSO models, both with and without the outlier, highlights the robustness of the t-link model compared to the probit link.

For all Bayesian models, we set τ 0 2 to 1 / n and select q = P ( γ i = 1 ) such that P ∑ j = 1 p γ j = 1 > K γ = 0.1, with K γ = max ( 10 , log ( n ) ). As discussed, unlike [32] and [38], the slab variance parameter τ 1 2 is imposed an inverse gamma prior with the hyperparameters s and r set to 2 and 1, respectively, providing a prior mean of 1 and an infinite prior variance. A total of 6, 000 iterations are run for the posterior sampling with a burn-in period of 2,000 iterations.

The metrics used for evaluating the methods are divided into two categories: variable selection performance and prediction performance. The variable selection metrics include sensitivity (SEN), specificity (SPE), and Matthew’s correlation coefficient (MCC), defined as follows: SEN = Sensitivity = T P T P + F N , SPE = Specificity = T N T N + F P , MCC = T P × T N − F N × F P ( T P + F P ) ( T P + F N ) ( T N + F P ) ( T N + F N ) , where TP, TN, FP, and FN represent true positives, true negatives, false positives, and false negatives, respectively. The MCC score provides a more comprehensive measure of performance than other metrics, as it incorporates all elements of the confusion matrix into a single summary statistic [13]. In addition, we present the out-of-sample mean square prediction error (MSPE) defined as: MSPE = | | E − p ˆ | | 2 2 n , where p ˆ is the estimated Prob ( E = 1 ).

Table 1 presents the comparison results, reporting average metrics from simulations based on data generated from the probit link model under two covariate correlation structures with p = 500. Results on the remaining 11 scenarios are provided in Tables S1–S5 of the Supplementary Materials.

Across all scenarios without data contamination (“No outlier”), the different link models exhibit similar performance. As expected, the true model generally performs slightly better than the others in terms of variable selection (higher MCC) and prediction accuracy (lower MSPE). Notably, the skinny Gibbs algorithm often outperforms the exact Gibbs algorithm for both the probit and t-link models, consistent with the findings of [32] and [38].

A comparison between the “No outlier” and “Bad Leverage Outlier” scenarios shows that the performance of the t-link models with ν = 1 or 3 experiences only minimal changes. In contrast, the results for the probit link show noticeable deterioration, both in variable selection (evidenced by lower MCC and higher FP and FN) and prediction (reflected in higher MSPE). Due to this one outlying observation, the probit model performs much worse than the t-link even when the datasets are generated under the probit link. These findings highlight the sensitivity of the probit link and the robustness of the t-link model. The t-link with ν = 7 and the logit link show generally comparable performance, both more robust than the probit link but less robust than the t-links with ν = 1 or ν = 3.

When comparing the “Non-leverage outlier” scenario with the other two cases, the t-links with ν = 1 and ν = 3 again demonstrate robustness in both variable selection and prediction. The other link functions show some deterioration, though to a lesser degree than under the “Bad Leverage Outlier” scenario, with the probit link remaining the most vulnerable.

Overall, the LASSO models perform worse than their Bayesian counterparts fitted with either the skinny Gibbs or exact Gibbs algorithms using the same link functions. Regarding the link functions, similar patterns are observed for the Cauchy, logit, and probit links under the LASSO method.

Figure 1 displays the runtime (in seconds) per iteration for the four methods, highlighting the computational efficiency of the skinny methods compared to the exact Gibbs methods. The computation time for HEGT and EProbit appears to increase quadratically with p, whereas HSGT and SProbit show a linear increase, consistent with their respective computational complexities. The t-link requires a slightly more time than the probit under either the exact or the skinny algorithms (by a margin of microseconds). This difference stems from the additional time spent sampling extra inverse-gamma latent variables, which is necessary when employing the t-link method.

The Fast Sampling Algorithm in [3] has a typical computational complexity of n 2 p, which is slightly higher than the n ( p ∨ | k | 2 ) complexity of Skinny Gibbs, where | k | denotes the active model size and typically satisfies | k | ≪ p. The Skinny Gibbs method is specifically designed to ensure scalability under sparsity, making it a more suitable choice for sparse Bayesian variable selection.

The PCR dataset originates from a study conducted by [27] to explore the genetics of two inbred mouse populations, B6 and BTBR. It includes gene expression data for 22,575 genes measured across 60 F2 mice, comprising 31 females and 29 males. The physiological phenotype, glycerol-3-phosphate acyltransferase (GPAT), was quantified using quantitative real-time Polymerase Chain Reaction (PCR). Both the phenotypic data (including GPAT measurements) and the gene expression data are publicly available on the GEO website (http://www.ncbi.nlm.nih.gov/geo) under accession number GSE3330.

GPAT is known to influence Hepatic Steatosis, a condition associated with obesity, with lower GPAT levels reducing its prevalence. Following the approach in [32], we define a binary response variable (E) based on GPAT values as E = I ( GPAT < Q ( 0.4 ) ), where Q ( 0.4 ) represents the 40th percentile of GPAT.

Given the high dimensionality of the gene expression data, an initial screening step was performed using the simple binary regression. This process identified the 99 most significant genes based on p-values, which were combined with the sex variable, resulting in a total of 100 predictors for further analysis. Figure 2 displays the boxplots of the standardized covariates in both datasets. The plot highlights the presence of numerous covariates with outlying values, which is a common occurrence in real-world data applications.

The six models, incorporating three link functions and two computational algorithms, are then fitted to the data. The initial number of active genes influencing GPAT is set to 10 ( | k | = 10), and the β vector is initialized to zero. Bayesian inference is performed using 4,000 MCMC samplings, following a burn-in period of 2,000 iterations. To assess the performance of these variable selection methods, we use 10-fold cross-validation and evaluate prediction accuracy based on mean squared prediction errors (MSPE), prediction accuracy and area under the curve (AUC). The dataset T is partitioned into 10 subsets, denoted T 1 , T 2 , … , T 10 . For each subset T r (where r = 1 , … , 10), variable selection is performed on the remaining data ( T ∖ T r ), and the model is used to compute predicted probabilities for the responses in T r . The MSPE for subset r is calculated as MSPE C V ( r ) = ∑ i ∈ S r ( E i − p ˆ i ) 2 / n r and Accuracy C V ( r ) = ( T P r + T N r ) / ( T P r + T N r + F P r + F N r ) where n r = | S r |, p ˆ i is the predicted probability for the ith observation, and E i is the observed response. The process is repeated for r = 1 , … , 10 and the overall cross-validated MSPE, accuracy and AUC are then obtained by averaging across all 10 subsets.

Figure 3 depicts the 10-fold cross-validation errors based on the selected number of covariates (model size). The MSPE C V is plotted against the different model sizes ranging from | k | = 1 to 6 with lower values indicating better predictive performance. It is evident that the HSGT and HEGT yield the smallest MSPE C V although the SProbit has the smallest MSPE when the number of active covariates is assumed to be one. The overall out-of-sample prediction accuracy in Figure 4 shows that the two t-link models demonstrate superior classification accuracy compared to the logit and probit methods except at model size 1.

Figure 5 illustrates the performance measured by AUC (Area Under the Curve) across varying model sizes. The plot shows that HSGT consistently achieves the highest AUC values as the model size increases, and has the highest AUC value from model size 3 to 6. The second best performed model is the HEGT. These results demonstrate the robustness of the t-link model and the practical performance of the scalable skinny Gibbs algorithm.

In this study, we establish the theoretical consistency results for the scalable algorithm to implement t-link (HSGT) in high-dimensional variable selection. The performance of different links and different algorithms is illustrated through simulation studies and an application to the real dataset. The t-link model is preferable, compared to the probit link, given its robustness in variable selection and prediction in the presence of outlying cases. HSGT methods offers competitive computational efficiency. As observed in [37], HSGT has a lower variability of the MCMC chains for active variables than HEGT and thus better performance in variable selection. Our application to real datasets further supports the effectiveness of HSGT. This outcome underscores its utility in handling high-dimensional sparse data, a common characteristic of modern datasets.

Our work opens the door to several promising research directions. First, it would be valuable to further investigate the interplay between link function choice and variable selection in high-dimensional data settings. While [35] explored this interplay using reversible jump Markov chain Monte Carlo (RJMCMC) to study the simultaneous selection of link functions and variables, we aim to address this question in high-dimensional settings by removing restrictions on the degrees of freedom v. Such an investigation would enhance the theoretical foundations and practical applications of Bayesian methods in modern data science.

Additionally, we are interested in establishing consistency results for more general link function models. We also plan to explore the use of Beta priors on q, as suggested by [41] and [11], and implement scalable spike-and-slab algorithms introduced by [4]. Another promising avenue is extending the methods to mixed-effects models, leveraging the normal scale mixture prior, as discussed by [45].

Robustness in generalized linear models can be pursued through the choice of link function or through robust estimation procedures such as M-estimators with adaptive weighting. In this study, we focus on the former. The weighting-based methods [23, 20] achieve robustness by downweighting observations with large residuals or leverage, thereby directly limiting the impact of contamination. Robustness by link modification stabilizes inference by the latent structures associated with the corresponding links. For example, [19] showed that the t-link’s influence function is bounded and its weight function down-weights extreme observations, thereby providing robust estimation under contamination and outperforming probit or logit links in such settings. While the weighting-based methods provides more flexibility, alternative links provide robustness while retaining standard estimation procedures.