The focus of this article centers around the linear model where the interest is to explore the linear association between a response variable and the covariates via Y ∼ N ( X β , σ 2 I ) where Y = ( Y 1 , … , Y n ) T is the n × 1 response variable, X = [ x 1 , … , x p ], is the n × p design matrix, and β = ( β 1 , … , β p ) T , is p × 1 vector of coefficients. The problem of variable selection can be treated as a model selection problem letting M ⊆ { all subsets of 1, …, p } as the model space under consideration. An additional notation of γ ⊂ { 0 , 1 } p is introduced to denote M γ , an individual member of M, indexed by the binary vector γ; while the null model which has no independent variable in the model is denoted by M 0 . The stochastic law of representation of Y then depends on ( x 1 , … , x p ) via X γ β γ where γ is working as a subscript of X ( β ) such that x j ( β j ) is present in the model whenever γ j = 1, j = 1 , … , p. It follows that there are 2 p models in the model space M, by virtue of which the model space easily becomes large even for moderate p thus precluding to visit all models in the model space.

As discussed before, our focus in this article is on the criterion based variable selection techniques as in any such criteria one must visit all the models in M to compare and conclude in the favor of a good model. Alternatively, one can find the model having a good, lowest in particular, value of a criterion by performing an optimization on the model space where the optimization can be carried out with respect to the criterion values of the candidate models. More generally, we consider a real valued function C on M, where C is the objective function which we want to minimize over M. We recall that, in variable selection setting, any model in the model space can be represented by the binary representation of γ. So when maximizing any objective function over the model space, the solution must belong to the set of binary numbers 0 , 1. This unique structure of the model space severely limits the choice of optimization methods.

Because of the special features of the maximization problem, we propose to conduct the maximization process stochastically using the widely known simulated annealing (SA), a stochastic optimization method. In what follows, we provide a brief review of an SA approach; for details, see [7]. We consider the model space M and define M ∗ ⊂ M to be the set of global minima of the function C, assumed to be a proper subset of M. For each i ∈ M, there exists a set M ( i ) ⊂ M ∖ { i }, called the set of neighbors of i. In addition, for every i, there exists a collection of positive coefficients q i j , j ∈ M ( i ), such that, ∑ j ∈ M ( i ) q i j = 1; so { q i j } = Q form a transition matrix, elements of which provide the transition probabilities of moving from i to j. It is assumed that j ∈ M ( i ) if and only if i ∈ M ( j ). We also define a nonincreasing function T : N → ( 0 , ∞ ) which is called the cooling schedule. Here N is the set of positive integers, and T ( t ) is called the temperature at time t.

Let ψ ( t ) be a discrete time inhomogeneous Markov chain on the model space M. The search process starts at an initial state ψ ( 0 ) ∈ M. Suppose at time t we arrive at the point i. We then choose a neighbor j of i at random according to probability q i j . Once j is chosen, and if C ( j ) ≤ C ( i ), then ψ ( t + 1 ) = j with probability 1; however if C ( j ) > C ( i ), then ψ ( t + 1 ) = j with probability exp [ − { C ( j ) − C ( i ) } / T ( t ) ], otherwise set ψ ( t + 1 ) = i; this gives raise to the so called Gibbs acceptance probability function. In supplementary material, we provide some technical clarity toward the performance of the proposed method.

[16] established that under regularity conditions, repeating the above steps with gradually reducing the temperature schedule guarantees that ψ ( t ) converges to the optimal set M ∗ , that is, for k ∈ N and all j ∈ S (2.1) lim n → ∞ Pr ( ψ ( n + k ) ∈ M ∗ | ψ ( k ) = i ) = lim n → ∞ Pr ( C ( ψ ( n + k ) ) (2.2) = C ∗ | ψ ( k ) = i ) = 1 where C ∗ = min j ∈ M C ( j ) and M ∗ = { i : i ∈ M , C ( i ) = C ∗ }.

The conditions for this result to hold are: ( i ) the probability of moving to j th model from i th model in p steps is positive, that is, q i j ( p ) > 0, q i i > 0 for all i, and ( i i i ) Q is irreducible.

The Bayesian approach to the variable selection problem is relatively straightforward. We express uncertainty about models by putting a prior distribution on the model M γ , The Bayesian linear model is thus defined as Y ∼ Pr ( y | θ , M γ ) , θ ∼ π ( θ | M γ ) , M γ ∼ Pr ( M γ ) , where π ( θ | M ) is the prior distribution of parameter θ = ( β , σ 2 ) under model M γ and Pr ( M γ ) is the prior on the model M γ . Then posterior distribution of θ is given by π ( θ | y , M γ ) = Pr ( y | θ , M γ ) π ( θ | M γ ) / Pr ( y | M γ ), where (2.3) Pr ( y | M γ ) = ∫ Pr ( y | θ , M γ ) π ( θ | M γ ) d θ is called the marginal likelihood or integrated likelihood of data y under model M γ . Then the posterior probability of model M γ can be expressed by (2.4) Pr ( M γ | y ) = Pr ( y | M γ ) Pr ( M γ ) ∑ M i ∈ M Pr ( y | M i ) Pr ( M i ) .

The Bayes factor [28, 3] for model M γ 1 against model M γ 2 is the ratio of their marginal likelihoods, BF 12 = Pr ( y | M γ 1 ) / Pr ( y | M γ 2 ). [28] stated that the Bayes factor is a summary of evidence for model M γ 1 against model M γ 2 and provided a table of cutoffs for interpreting log BF 12 . In general, the model with higher log-marginal likelihood is preferable in this model selection criterion.

In modern era, Bayesian inference is typically done by Markov Chain sampling. The computation of Bayes factor from Markov Chain sampling, however, is generally difficult since the Markov Chain methods avoid the computation of the normalizing constant of the posterior and it is precisely this constant that is needed for the marginal likelihood.

The HPM has the highest posterior model probability among all models in the model space, that is, HPM = argmax γ ∈ M Pr ( M γ | y _ ). Under the notion of a data generating model (or the so-called true model) in the model space it can be shown that the data generating model is often asymptotically equivalent to the highest posterior model. For instance, this can be examined via consistency of posterior model probabilities [20] or via the Bayes factors [32]. [20] examined model consistency for g priors when g is fixed. [30] extended this for mixture of g priors and hyper g priors. [17] proved model consistency for spike and slab type priors. [12] and [32] proved consistency of objective Bayes procedures. On the other hand, [39] and [40] showed Bayes factor consistency for unbalanced ANOVA models and nested designs respectively. Moreover, [27] proved consistency for the true model when non local priors were specified on the parameters. However, they distinguished the true model consistency and pairwise Bayes factor consistency and argued that for large dimensional space pairwise consistency is misleading and hence not much useful.

We set C = negative posterior probability of model M γ for maximizing the posterior probabilities over the model space applying simulated annealing algorithm. In the SA approach, an appropriately chosen cooling schedule accelerates convergence. When T is very small, the time it takes for the Markov chain ψ ( t ) to reach equilibrium can be excessive. The main significance of cooling schedule is that, during the beginning of the search process it helps the algorithm to escape from the local modes and then when the search is actually in the neighborhood of the global optimum the algorithm tries to focus in that region by reducing the value of cooling schedule and thereby finding the actual optimum. There is a number of suggestions available in the literature to choose a functional form for cooling schedule.

A transition matrix definition is equally important in an SA algorithm. We define the ( i , j ) th element of the transition matrix Q as q i j = posterior probability of j th model sum of posterior probabilities of neighbors of i th model where j th model ∈ neighborhood of i th model.

For a given model M γ we define its neighborhood as { M γ 0 , M γ 00 }, where γ 0 is such that if | γ 0 − γ | = 1, that is, the model M γ 0 can be obtained from model M γ by either adding or deleting one predictor; γ 00 ′ 1 = γ ′ 1 and | γ 00 − γ | = 2, that is, model M γ 00 can be obtained from model M γ by swapping one predictor with another.

It is interesting to note that, our selection provides the advantage for getting different region of neighborhood at every step and thus eliminates the possibility of keeping old models in the search region which is the case in [5] and [24]. In this way our approach is different in the sense that the search procedure does not require a complicated and long Markov chain to converge. These ingredients give raise to our proposed stochastic search algorithm called SA-HPM the steps of which are described below. The approach is implemented in R package sahpm and is made available on R CRAN.

In practice, to make the computation stable, we suggest to calculate the log of posterior probabilities instead of posterior probabilities and use that as estimates of proposal distributions.

The ozone dataset has been considered in the literature frequently [5, 11] and consists of daily measurements of atmospheric ozone concentration (maximum one hour average) and eight meteorological quantities for 330 days of 1976 in the Los Angeles Basin. Among them one temperature predictor was dropped from the analysis due to the potential multicollinearity with another temperature variable. The Ozone35 data was then curated by considering the main effects, the second order effects, their first order interactions [21], which gives raise to a total of p = 35 covariates. Mainly for comparison purpose we make use of the n = 178 observations which were used in the analysis of [21]. The description of the predictor variables and the response variable is provided in Table 3. [21] illustrated that the posterior probability of the median probability model (MPM [2]) is 23 times lower than that of the highest posterior model. We considered a g-prior as in [21].

[21] considered a complete enumeration of this large model space using distributed computing over an extended time and reported the complete enumeration HPM to be the model (7, 10, 23, 26, 29). The Bayes factor of HPM and MPM, against M 0 are reported in Table 4. Our main contribution is that, the proposed SA-HPM method is able to recover the HPM 95 times out of 100 repetitions after a burn-in of 50 iterations in the stochastic chain. In particular, we note that, the SA-HPM method is extremely useful even for large model spaces.

In this example we consider gene expression data on 31 female mice and 29 male mice. A number of psychological phenotypes, including numbers of stearoyl-CoA desaturase 1 (SCD1), glycerol-3-phosphate acyltransferase (GPAT) and phos- phoenopyruvate carboxykinase (PEPCK), were measured by quantitative real-time RT-PCR, along with 22,575 gene expression values. The resulting data set is publicly available at http://www.ncbi.nlm.nih.gov/geo (accession number GSE3330). Following [35] we restrict ourselves into the consideration of the SCD1 response only.

Due to ultra large high-dimensional nature of this dataset it is beyond the reach of the ultra modern machinery to enumerate all the models in the model space. Hence, in order to find the highest posterior model it is necessary to make use of a model space search technique such as the SA-HPM method developed here. We employ our algorithm in this dataset to find the HPM model. When utilizing SA-HPM we omit the swapping step to minimize exploring the model space due to the ultra large size of that. We report our findings in Table 5 from which it can be noted that the resulting HPM is a sparse model with three variables when SA-HPM with g prior is fitted. Similarly, SA-HPM with piMOM prior discovers another sparse model with two predictors. It is interesting to note that one of them coincides with the model projected by the maximum aposteriori (MAP) estimate of the BayesS5 method. We notice that BayesS5 also results in a parsimonious model with two predictor variables.