BSGS: Bayesian Sparse Group Selection

An R package BSGS is provided for the integration of Bayesian variable and sparse group selection separately proposed by Chen et al. (2011) and Chen et al. (in press) for variable selection problems, even in the cases of large p and small n. This package is designed for variable selection problems including the identification of the important groups of variables and the active variables within the important groups. This article introduces the functions in the BSGS package that can be used to perform sparse group selection as well as variable selection through simulation studies and real data.


Introduction
Variable selection is a fundamental problem in regression analysis, and one that has become even more relevant in current applications where the number of variables can be very large, but it is commonly assumed that only a small number of variables are important for explaining the response variable.This sparsity assumption enables us to select the important variables, even in situations where the number of candidate variables is much greater than the number of observations.BSGS is an R package designed to carry out a variety of Markov chain Monte Carlo (MCMC) sampling approaches for variable selection problems in regression models based on a Bayesian framework.In this package, we consider two structures of variables and create functions for the corresponding MCMC sampling procedures.In the first case where the variables are treated individually without grouping structure, two functions, CompWiseGibbsSimple and CompWiseGibbsSMP, are provided to generate the samples from the corresponding posterior distribution.In the second case, it is assumed that the variables form certain group structures or patterns, and thus the variables can be partitioned into different disjoint groups.However, only a small number of groups are assumed to be important for explaining the response variable, i.e. the condition of the group sparsity, and we also assume that sparse assumption is held for the variables within the groups.This problem is thus termed a sparse group selection problem Simon et al. (2013); Chen et al. (in press), and the goal is to select the important groups and also identify the active variables within these important groups simultaneously.There are two functions to handle the sparse group selection problems, BSGS.Simple and BSGS.Sample, which are used to generate the corresponding posterior samples.Once the posterior samples are available, we then can determine the active groups and variables, estimate the parameters of interest and make other statistical inferences.This paper is organized as follows.We first briefly introduce statistical models that are used to deal with the problems of variable selection in the BSGS package.We then describe the tuning parameters in the functions in the BSGS package.Two simulations are used to illustrate the details of the implementations of the functions.Finally we present a real economic example to demonstrate the BSGS package.

Framework of BSGS
We start with the introduction of individual variable selection problems, and then turn our attention to sparse group selection.For completeness, we describe the model and priors so that one may easily change the inputs of functions in the BSGS package for any purpose.

Variable selection
Consider a linear regression model given by where Y = (Y 1 , . . ., Y n ) is the response vector of length n, X = [X 1 , . . ., X p ] is an n × p design matrix, with X i as the corresponding i-th variable (regressor) as a potential cause of the variation in the response, β is the corresponding unknown p × 1 coefficient vector, and ε = (ε 1 , . . ., ε n ) is the error term, which is assumed to have a normal distribution with a zero mean, and the covariance matrix σ 2 I n , I n is the n × n identity matrix.To achieve variable selection, we select a subset of X 1 , . . ., X p to explain the response Y.For this purpose, following the stochastic search variable selection method George and McCulloch (1993), a latent binary variable γ i taking the value of 0 and 1 is introduced The R Journal Vol.7/2, December 2015 ISSN 2073-4859 to indicate whether or not the corresponding variable is selected.That is, if γ i = 1, the variable X i is selected, and otherwise it is not selected.
In this Bayesian framework we basically follow the prior assumption in Chen et al. (2011).We assume the prior distribution of γ i is a Bernoulli distribution with probability 1 − ρ i , and then given γ i , we assume the prior for β i is a mixture of point mass at 0 denoted by δ 0 and a normal distribution as follows where τ i is a pre-specified constant.Moreover, (γ i , β i ) are assumed to be independent for i = 1, . . ., p. Lastly, σ 2 is set to follow an inverse Gamma distribution, Based on the model setting given above, there are two Gibbs samplers for this variable selection.The first procedure is the componentwise Gibbs sampler (CGS) which was introduced by Geweke (1996) and also mentioned in Chen et al. (2011).The other is the stochastic matching pursuit (SMP) algorithm Chen et al. (2011) in which the variables compete to be selected.

Sparse group selection
In traditional variable selection problems, each variable X i is treated individually; however, in some real applications, a group of variables that behave similarly may be more meaningful in explaining the response.In other words, a group of variables potentially plays a more important role in explaining the response than a single variable can.The variable selection problem thus becomes a group selection one.In group selection problems, the variables within the selected groups are all treated as important.However, this assumption might not be held in practice such as the climate application in Chatterjee et al. (2012).Instead of the group selection problem, we thus consider approaches to sparse group selection, in which the sparse assumption is held for groups and the variables within groups.
Here the goal is not only to select the influential groups, but also to identify the important variables within these.To this end, a Bayesian group variable selection approach, the Group-wise Gibbs sampler (GWGS) Chen et al. (in press), is applied.Suppose that, in terms of expert or prior knowledge, the explanatory variables X i 's are partitioned into g non-overlapping groups in a regression model.Each group l contains j = 1, . . ., p l variables with ∑ g l=1 p l = p.Now the model is rewritten as where is the p l × 1 coefficient vector of group l.The GWGS works by introducing two nested layers of binary variables to indicate whether a group or a variable is selected or not.At the group-selection level, we define a binary variable η l = 1 if group l is selected, and η l = 0 otherwise.At the variable-selection level, we define another binary variable γ li = 1 if the variable i within group l, X li , is selected, and γ li = 0 otherwise.We assume the group indicator, η l , has the Bernoulli distribution with the probability 1 − θ l .Within group l, the prior distribution of γ li conditional on the indicator η l is defined as where δ 0 is a point mass at 0 and Ber(1 − ρ li ) is Bernoulli distributed with the probability 1 − ρ li .Equation ( 5) implies that if the l-th group is not selected, it turns out that γ li = 0 for all i.The prior distribution of the coefficient β li given η l and γ li is given by where τ li is a pre-specified value Chen et al. (2011).Finally, the variance σ 2 is assumed to have an inverse Gamma distribution, that is, σ 2 ∼ IG(ν/2, νλ/2).We also assume (η l , γ l1 , β l1 , . . ., γ l p l , β l p l ), l = 1, • • • , g, are a priori independent.
Two sampling procedures are proposed in Chen et al. (in press) for Bayesian sparse group selection.The first is the GWGS.In the GWGS, simulating the indicator variable η l from the posterior distribution would be computationally intensive, especially when the number of variables within the group is large.To address this issue, Chen et al. (in press) proposed a modified and approximation approach, a sample version of GWGS.In this a Metropolis-Hastings algorithm is adopted to replace the Gibbs sampling method in GWGS.

Implementation
In this section, we describe the default tuning parameters and some details in the implementation of the functions in BSGS.

Hyperparameter set-up
• The tuning parameters, ν and λ: The parameters in the prior distribution of σ 2 are suggested by George and McCulloch (1993) setting the default values of ν = 0 and λ being any positive value.A data-driven choice for this is proposed by Chipman et al. (1997), which sets the value of ν around 2 and the value of λ is set up to be the 99% quantile of the prior of σ 2 that is close to Var(Y).In addition, George and McCulloch (1997) simply set ν = 10 and λ = Var(Y).In our experiences, a larger value of ν tends to result in a larger estimate of σ.
• The parameter τ: Now we consider the assignment of the value of τ, the prior variance of the regression coefficient for the active variable.It was found that the larger the value of τ, the smaller the conditional probability of γ = 1 is.As a result, a large value of τ favors a more parsimonious mode.In contrast, a small value of τ would yield more complex models.
• The parameters ρ and θ: The default of the prior inclusion probabilities of groups and variables is set equal to 0.5.In the case when p is much greater than n and only a small number of variables are considered active, we would assign a larger values to ρ and θ to reflect the prior belief of sparsity.

Stopping rule
The posterior distribution is not available in explicit form so we use the MCMC method, and specifically Gibbs sampling to simulate the parameters from this distribution Brooks et al. (2011).To implement the Gibbs sampler, the full conditional distributions of all parameters must be determined.A derivation of the full conditional distributions is provided in Chen et al. (in press).When these have been obtained, the parameters are then updated individually using a Gibbs sampler (where available), or a Metropolis-Hastings sampling algorithm.An MCMC sample will converge to the stationary distribution, i.e. the posterior distribution.We use the batch mean method to estimate the Monte Carlo standard error (MCSE) of the estimate of σ 2 and then decide to stop the simulation once the MCSE is less than a specified value, cf.Flegal et al. (2008).The default minimum number of iterations is 1000.
If the MCSE does not achieve the prespecified value, an extra 100 iterations are run until the MCSE is less than the prespecified value.The sample can then be used for statistical inference.

Statistical inference
• Variable and group selection criteria: In this package, we adopt the median probability criterion proposed by Barbieri and Berger (2004) for group and variable selections.Specifically, for the variable selection problem, we estimate the posterior inclusion probability P(γ i = 1|Y) from the posterior samples and then the i-th variable is selected into the model if the estimated posterior probability is larger than or equal to 1/2.Here instead of 1/2, this cut-off value is treated as a tuning parameter, α, which can be specified by users.For the sparse group selection problem, the estimated posterior probability of the l-th group is greater than or equal to α g , i.e.P(η l = 1|Y) ≥ α g , we then include X l into the model.Suppose the l-th group is selected, then the i-th variable within this group is selected if Here α g and α i are two pre-specified values between 0 and 1.

• Posterior estimates of regression coefficients:
We use the Rao-Blackwell method to estimate β by where β m is the sample in mth iteration, M is the number of iterations, and N M is the number of nonzero β m .

Evaluation of model estimation
Regarding the stability of the estimation, we compare the accuracy of selection of the variables by the following measures in the simulation studies: the True Classification Rate (TCR), the True Positive Rate (TPR), and the False Positive Rate (FPR).These are defined as follows TCR = number of correctly selected variables number of variables ; TPR = number of correctly selected variables number of active variables ; FPR = number of falsely selected variables number of inactive variables .
TCR is an overall evaluation of accuracy in the identification of the active and inactive variables.TPR is the average rate of active variables correctly identified, and is used to measure the power of the method.FPR is the average rate of inactive variables that are included in the regression, and it can be considered as the type I error rate of the selection approach.In these three criteria, it is preferred to have a larger value of TCR or TPR, or a smaller value of FPR.
We also report the mean squared error (MSE), where ŷi is the prediction of y i .This is used to evaluate whether the overall model estimation has a good fit with regard to the data set.

Examples
Two simulations and a real example are provided to demonstrate the use of functions in the BSGS package.

Simulation II
We provide another stimulation to illustrate the use of functions for sparse group selection.In the following simulation, the data Y of length n = 50 is generated from a normal distribution with a mean Xβ and σ 2 = 1, where X is from a multivariate distribution with mean 0 and covariance Σ with the pair correlation between variables equal to zero.There are 10 groups of variables.Each group contains 10 variables and some of them are active.More specifically, X l = X l1 , . . ., X l p l , k = 1, . . ., 10 and p l = 10 for all l.We assume the group l = 1, 2, 5, 8 are active.Variables p 1 = 7, 8, 9 in the group l = 1, p 2 = 1, 2 in the group l = 2, p 5 = 3 in the group l = 5, and p 8 = 7 in the group l = 8 are active.
## hyperparameters tau2.value<-rep(1, num.of.covariates) rho.value <-rep(0.5, num.of.covariates) theta.value<-rep(0.5, num.of.group.var)With the reasonable assignment of the initial values of parameters, we apply two functions, BSGS.Simple and BSGS.Sample, to estimate the posterior quantities of interest and in turn identify the important groups and variables.For illustration, we stop the simulation when the MCSE of estimate of σ 2 is less than 0.5.
## The posterior quantities estimated by two sampling approaches respectively.
One may be interested in how the different stopping points would affect the computational effort and model estimation.We thus compare the computational times and accuracy of parameter estimation in terms of TCR, TPR, FPR, and MSE for two different sampling approaches by using different MCSEs to stop the simulation and the results are shown in Tables 1 and 2. It has been found that for the case of n = 50 and p = 100 the simple version would converge faster than the sample version.However, the sample version would require less effort to produce the sample than the simple version.We further investigate the accuracy of parameter estimations and the selection of the variables in terms of MSE, TCR, TPR, and FPR for the two sampling approaches when the number of covariates increases but each group has 10 variables.We consider p = 300, 500, and 1000 and the simulation is terminated when the MCSE of the estimate of σ 2 is less than 0.5 for illustration.One may use different values of MCSE, but more computational time may be needed.Table 3 shows that the simulation stops with the same number of iterations and the parameter estimates and accuracy of variable selection show little difference.On the other hand, the results in Table 4 for the sample version GWGS show that more iterations are needed, under the same stopping rule.The two approaches thus perform equally well on the selection of variables, but by comparing MSEs it is evident that the simple version of GWGS outperforms with regard to the predictions.Although the sample version has less computational intensity, it needs more iterations to achieve the stopping point.If one is interested in selecting important variables, both approaches are effective.But if one is interested in choosing a model which fits the data well, it is thus suggested one uses the simple version approach.

MCSE # of iterations runtime (in
Next, we compare the computational time when the number of variables in the group increases.We perform an experiment in which there are seven groups, each containing 15 variables.The assignments of hyperparameters are the same as those in Simulation II.

A real economic example
This subsection further illustrates the functions in the BSGS based on an economic dataset from Rose and Spiegel Rose and Spiegel (2010, 2011, 2012) which is available at http://faculty.haas.berkeley.edu/arose The response variable is the 2008-2009 growth rate for the crisis measure.Rose and Spiegel originally consider 119 explanatory factors for the crisis for as many as 107 countries, but there are missing data for a number of these.

require(BSGS) ## the whole data set data(Crisis2008)
To maintain a balanced data set, we use 51 variables for a sample of 72 countries.For more information about the balanced data, please see the description of 'Crisis2008' in the BSGS.The balanced data is then analyzed to illustrate the main sampling function BSGS.Simple to simulate the sample from the posterior distribution.In the analysis, we demean the response so that it is not necessary to include the intercept into the design matrix.All variables are standardized except the dummy variables.As discussed in Ho (in press), these variables can be classified into the nine theoretical groups of the crisis' origin (the number in parentheses indicates the number of variables considered in the group): principal factors (10), financial policies (three), financial conditions (four), asset price appreciation (two), macroeconomic policies (four), institutions (11), geography (four), financial linkages (one), and trade linkages ( 12).Based on this information, we assign a group index to each variable.

set.seed(1) data
Group.Index <-rep(1:9, c(10, 3, 4, 2, 4, 11, 4, 1, 12)) Since the number of covariates within the group is moderate, it is recommended that the simple version GWGS be applied to generate the samples.In this example, we have tested different values of τ 2 and finally we set τ 2 = 10 due to the minimal MSE's.Here the stopping rule is when the MCSE of the estimate of σ 2 is less than or equal to 0.1 for the simple version.No prior information is provided to indicate which group or variable is more important, so we let θ = 0.5 and η = 0.5.We let ν = 0 and λ = 1 resulting in a non-informative prior for σ 2 .In each group, we will update parameters "Num.of.Iter.Inside.CompWise = 100" times within a group for stability.

Summary
This paper illustrated the usage of a new R package, BSGS, for identifying the important groups of variables and important variables in linear regression models.Furthermore, BSGS can be easily implemented with problems of a large p and small n.The grouping idea is also applicable to other regression and classification settings, for example, the multi-response regression and multi-class classification problems.We envision future additions to the package that will allow for extensions to these models.
We are confident that this package can be applied to many important real-world problems by keeping flexibility with regard to selecting variables within a group based on the hierarchical assignment of two layers of indicator variables.For instance, in the gene-set selection problem, a biological The R Journal Vol.7/2, December 2015 ISSN 2073-4859 Once the simulation stops, the posterior samples are used to estimate the posterior quantities of interest.One can then check the number of iterations and computational times for both approaches.One can then use the function CGS.SMP.PE to identify the important variables and to estimate the parameters.Due to the limitations of space, we do not include the estimates here.A variable i is considered important if the posterior probability of its indicator variable γ i = 1 is greater than or equal to α i = 1/2.Once the critical point is decided, two functions TCR.TPR.FPR.CGS.SMP and MSE.CGS.SMP are carried out to evaluate the performance on the model estimations in terms of TCP, TPR, FTP and MSE.

Table 1 :
Results based on simple version procedure for sparse group selection for different stopping points.

Table 2 :
Results based on sample version procedure for sparse group selection for different stopping points.
Table 5 shows the computational time plus

Table 3 :
Results based on the simple version procedure for sparse group selection for different numbers of covariates, with pair-correlation equal to 0.

Table 4 :
Results based on the sample version procedure for sparse group selection for different numbers of covariates, with pair-correlation equal to 0. the model estimations.It can be seen that the sample version is strongly recommended when the number of variables within a group is greater than 15.

Table 5 :
Comparison between simple and sample versions for sparse group selection when the number of covariates within a group is 15.