Evolutionary Stochastic Search for Bayesian model exploration

Introduction

This paper is a contribution to the methodology of Bayesian variable selection for linear Gaussian regression models, an important problem which has been much discussed both from a theoretical and a practical perspective (see Chipman et al., 2001 and Clyde and George, 2004 for literature reviews). Recent advances have been made in two directions, unravelling the theoretical properties of different choices of prior structure for the regression coefficients (Fernández et al., 2001; Liang et al., 2008) and proposing algorithms that can explore the huge model space consisting of all the possible subsets when there are a large number of covariates, using either MCMC or other search algorithms (Kohn et al., 2001; Dellaportas et al., 2002; Hans et al., 2007).

In this paper, we propose a new sampling algorithm for implementing the variable selection model, based on tailoring ideas from Evolutionary Monte Carlo (Liang and Wong, 2000; Jasra et al., 2007; Wilson et al., 2009) in order to overcome the known difficulties that MCMC samplers face in a high dimension multimodal model space: enumerating the model space becomes rapidly unfeasible even for a moderate number of covariates. For a Bayesian approach to be operational, it needs to be accompanied by an algorithm that samples the indicators of the selected subsets of covariates, together with any other parameters that have not been integrated out. Our new algorithm for searching through the model space has many generic features that are of interest per se and can be easily coupled with any prior formulation for the variance-covariance of the regression coefficients. We illustrate this by implementing $g$-priors for the regression coefficients as well as independent priors: in both cases the formulation we adopt is general and allows the specification of a further level of hierarchy on the priors for the regression coefficients, if so desired.

The paper is structured as follows. In Section 2, we present the background of Bayesian variable selection, reviewing briefly alternative prior specifications for the regression coefficients, namely $g$-priors and independent priors. Section 3 is devoted to the description of our MCMC sampler which uses a wide portfolio of moves, including two proposed new ones. Section 4 demonstrates the good performance of our new MCMC algorithm in a variety of real and simulated examples with different structures on the predictors. In Section 4.2 we complement the results of the simulation study by comparing our algorithm with the recent Shotgun Stochastic Search algorithm of Hans et al. (2007). Finally Section 5 contains some concluding remarks.

Background

Let $y=\left(y_{1},\ldots,y_{n}\right)^{T}$ be a sequence of $n$ observed responses and $x_{i}=\left(x_{i1},\ldots,x_{ip}\right)^{T}$ a vector of predictors for $y_{i}$, $i=1,\ldots,n$, of dimension $p\times 1$. Moreover let $X$ be the $n\times p$ design matrix with $i$th row $x_{i}^{T}$. A Gaussian linear model can be described by the equation

where $\alpha$ is an unknown constant, $1_{n}$ is a column vector of ones, $\beta=\left(\beta_{1},\ldots,\beta_{p}\right)^{T}$ is a $p\times 1$ vector of unknown parameters and $\varepsilon\sim N\left(0,\sigma^{2}I_{n}\right)$.

Suppose one wants to model the relationship between $y$ and a subset of $x_{1},\ldots,x_{p}$, but there is uncertainty about which subset to use. Following the usual convention of only considering models that have the intercept $\alpha$, this problem, known as variable selection or subset selection, is particularly interesting when $p$ is large and parsimonious models containing only a few predictors are sought to gain interpretability. From a Bayesian perspective the problem is tackled by placing a constant prior density on $\alpha$ and a prior on $\beta$ which depends on a latent binary vector $\gamma=\left(\gamma_{1},\ldots,\gamma_{p}\right)^{T}$, where $\gamma_{j}=1$ if $\beta_{j}\neq 0$ and $\gamma_{j}=0$ if $\beta_{j}=0$, $j=1,\ldots,p$. The overall number of possible models defined through $\gamma$ grows exponentially with $p$ and selecting the best model that predicts $y$ is equivalent to find one over the $2^{p}$ subsets that form the model space.

Given the latent variable $\gamma$, a Gaussian linear model can therefore be written as

where $\beta_{\gamma}$ is the non-zero vector of coefficients extracted from $\beta$, $X_{\gamma}$ is the design matrix of dimension $n\times p_{\gamma}$, $p_{\gamma}\equiv\gamma^{T}1_{p}$, with columns corresponding to $\gamma_{j}=1$. We will assume that, apart from the intercept $\alpha$, $x_{1},\ldots,x_{p}$ contains no variables that would be included in every possible model and that the columns of the design matrix have all been centred with mean .

It is recommended to treat the intercept separately and assign it a constant prior: $p\left(\alpha\right)\propto 1$, Fernández et al. (2001). When coupled with the latent variable $\gamma$, the conjugate prior structure of $\left(\beta_{\gamma},\sigma^{2}\right)$ follows a normal-inverse-gamma distribution

with $a_{\sigma},b_{\sigma}>0$. Some guidelines on how to fix the value of the hyperparameters $a_{\sigma}$ and $b_{\sigma}$ are provided in Kohn et al. (2001), while the case $a_{\sigma}=b_{\sigma}=0$ corresponds to the Jeffreys’ prior for the error variance, $p\left(\sigma^{2}\right)\propto\sigma^{-2}$. Taking into account (1), (2), (3) and the prior specification for $\alpha$, the joint distribution of all the variables (based on further conditional independence conditions) can be written as

The main advantage of the conjugate structure (2) and (3) is the analytical tractability of the marginal likelihood whatever the specification of the prior covariance matrix $\Sigma_{\gamma}$:

where $S\left(\gamma\right)=C-M^{T}K_{\gamma}^{-1}M$, with $C=\left(y-\bar{y}_{n}\right)^{T}\left(y-\bar{y}_{n}\right)+m_{\gamma}^{T}\Sigma_{\gamma}^{-1}m_{\gamma}$, $M=X_{\gamma}^{T}\left(y-\bar{y}_{n}\right)+\Sigma_{\gamma}^{-1}m_{\gamma}$ and $K_{\gamma}=X_{\gamma}^{T}X_{\gamma}+\Sigma_{\gamma}^{-1}$ (Brown et al., 1998).

While the mean of the prior (2) is usually set equal to zero, $m_{\gamma}=0$, a neutral choice (Chipman et al., 2001; Clyde and George, 2004), the specification of the prior covariance $\Sigma_{\gamma}$ matrix leads to at least two different classes of priors:

When $\Sigma_{\gamma}=gV_{\gamma}$, where $g$ is a scalar and $V_{\gamma}=\left(X_{\gamma}^{T}X_{\gamma}\right)^{-1}$, it replicates the covariance structure of the likelihood giving rise to so called $g$-priors first proposed by Zellner (1986).

When $\Sigma_{\gamma}=cV_{\gamma}$, but $V_{\gamma}=I_{p_{\gamma}}$ the components of $\beta_{\gamma}$ are conditionally independent and the posterior covariance matrix is driven towards the independence case.

We will adopt the notation $\Sigma_{\gamma}=\tau V_{\gamma}$ as we want to cover both prior specification in a unified manner. Thus in the $g$-prior case, $\Sigma_{\gamma}=\tau\left(X_{\gamma}^{T}X_{\gamma}\right)^{-1}$ while in the independent case, $\Sigma_{\gamma}=\tau I_{p_{\gamma}}$. We will refer to $\tau$ as the variable selection coefficient for reasons that will become clear in the next Section.

To complete the prior specification in (4), $p\left(\gamma\right)$ must be defined. A complete discussion about alternative priors on the model space can be found in Chipman (1996) and Chipman et al. (2001). Here we adopt the beta-binomial prior illustrated in Kohn et al. (2001)

with $p_{\gamma}\equiv\gamma^{T}1_{p}$, where the choice $p\left(\gamma\left|\omega\right.\right)=\omega^{p_{\gamma}}\left(1-\omega\right)^{p-p_{\gamma}}$ implicitly induces a binomial prior distribution over the model size and $p\left(\omega\right)=\omega^{a_{\omega}-1}\left(1-\omega\right)^{b_{\omega}-1}/B\left(a_{\omega},b_{\omega}\right)$. The hypercoefficients $a_{\omega}$ and $b_{\omega}$ can be chosen once $E\left(p_{\gamma}\right)$ and $V\left(p_{\gamma}\right)$ have been elicited. In the “large $p$, small $n$” framework, to ensure sparse regression models where $p_{\gamma}\ll p$, it is recommended to centre the prior for the model size away from the number of observations.

2 Priors for the variable selection coefficient τ𝜏\tau

It is a known fact that $g$-priors have two attractive properties. Firstly they possess an automatic scaling feature (Chipman et al., 2001; Kohn et al., 2001). In contrast, for independent priors, the effect of $V_{\gamma}=I_{p_{\gamma}}$ on the posterior distribution depends on the relative scale of $X$ and standardisation of the design matrix to units of standard deviation is recommended. However, this is not always the best procedure when $X$ is possibly skewed, or when the columns of $X$ are not defined on a common scale of measurement. The second feature that makes $g$-priors particularly appealing is the rather simple structure of the marginal likelihood (5) with respect to the constant $\tau$ which becomes

where, if $m_{\gamma}=0$, $S\left(\gamma\right)=\left(y-\bar{y}_{n}\right)^{T}\left(y-\bar{y}_{n}\right)-\frac{\tau}{1+\tau}\left(y-\bar{y}_{n}\right)^{T}X_{\gamma}\left(X_{\gamma}^{T}X_{\gamma}\right)^{-1}X_{\gamma}^{T}\left(y-\bar{y}_{n}\right)$. For computational reasons explained in the next Section, we assume that (7) is always defined: since we calculate $S\left(\gamma\right)$ using the QR-decomposition of the regression $\left(X_{\gamma},y-\bar{y}_{n}\right)$ (Brown et al., 1998), when $n\leq p_{\gamma}$, $S\left(\gamma\right)=\left(y-\bar{y}_{n}\right)^{T}\left(y-\bar{y}_{n}\right)/\left(1+\tau\right)$. Despite the simplicity of (7), the choice of the constant $\tau$ for $g$-priors is complex, see Fernández et al. (2001), Cui and George (2008) and Liang et al. (2008).

Historically the first attempt to build a comprehensive Bayesian analysis placing a prior distribution on $\tau$ dates back to Zellner and Siow (1980), where the data adaptivity of the degree of shrinkage adapts to different scenarios better than assuming standard fixed values. Zellner-Siow priors, Z-S hereafter, can be thought as a mixture of $g$-priors and an inverse-gamma prior on $\tau$, $\tau\sim InvGa(1/2,n/2)$, leading to

Liang et al. (2008) analyse in details Z-S priors pointing out a variety of theoretical properties. From a computational point of view, with Z-S priors, the marginal likelihood $p\left(y\left|\gamma\right.\right)=\int p\left(y\left|\gamma,\tau\right.\right)p\left(\tau\right)d\tau$ is no more available in closed form, something which is advantageous in order to quickly perform a stochastic search (Chipman et al., 2001). Even though Z-S priors need no calibration and the Laplace approximation can be derived (Tierney and Kadane, 1986), see Appendix A.2, never became as popular as $g$-priors with a suitable constant value for $\tau$. For alternative priors, see also Cui and George (2008) and Liang et al. (2008).

2.2 Independent priors

When all the variables are defined on the same scale, independent priors represent an attractive alternative to $g$-priors. The likelihood marginalised over $\alpha$, $\beta_{\gamma}$ and $\sigma^{2}$ becomes

where, if $m_{\gamma}=0$, $S\left(\gamma\right)=\left(y-\bar{y}_{n}\right)^{T}\left(y-\bar{y}_{n}\right)-\left(y-\bar{y}_{n}\right)^{T}X_{\gamma}\left(X_{\gamma}^{T}X_{\gamma}+\tau I_{p_{\gamma}}\right)^{-1}X_{\gamma}^{T}\left(y-\bar{y}_{n}\right)$. Note that (9) is computationally more demanding than (7) due to the extra determinant operator.

Geweke (1996) suggests to fix a different value of $\tau_{j}$, $j=1,\ldots,p$, based on the idea of “substantially significant determinant” of $\Delta X_{j}$ with respect to $\Delta y$. However it is common practice to standardise the predictor variables, taking $\tau=1$ in order to place appropriate prior mass on reasonable values of the regression coefficients (Hans et al., 2007). Another approach, illustrated in Bae and Mallick (2004), places a prior distribution on $\tau_{j}$ without standardising the predictors.

Regardless of the prior specification for $\tau$, using the QR-decomposition on a suitable transformation of $X_{\gamma}$ and $y-\bar{y}_{n}$, the marginal likelihood (9) is always defined.

MCMC sampler

In this Section we propose a new sampling algorithm that overcomes the known difficulties faced by MCMC schemes when attempting to sample a high dimension multimodal space. We discuss in a unified manner the general case where a hyperprior on the variable selection coefficient $\tau$ is specified. This encompasses the $g$-prior and independent prior structure as well as the case of fixed $\tau$ if a point mass prior is used.

The multimodality of the model space is a known issue in variable selection and several ways to tackle this problem have been proposed in the past few years. Liang and Wong (2000) suggest an extension of parallel tempering called Evolutionary Monte Carlo, EMC hereafter, Nott and Green, N&G hereafter, (2004) introduce a sampling scheme inspired by the Swendsen-Wang algorithm while Jasra et al. (2007) extend EMC methods to varying dimension algorithms. Finally Hans et al. (2007) propose when $p>n$ a new stochastic search algorithm, SSS, to explore models that are in the same neighbourhood in order to quickly find the best combination of predictors.

We propose to solve the issue related to the multimodality of model space (and the dependence between $\gamma$ and $\tau$) along the lines of EMC, applying some suitable parallel tempering strategies directly on $p\left(y\left|\gamma,\tau\right.\right)$. The basic idea of parallel tempering, PT hereafter, is to weaken the dependence of a function from its parameters by adding an extra one called “temperature”. Multiple Markov chains, called “population” of chains, are run in parallel, where a different temperature is attached to each chain, their state is tentatively swap at every sweep by a probabilistic mechanism and the latent binary vector $\gamma$ of the non-heated chain is recorded. The different temperatures have the effect of flatting the likelihood. This ensures that the posterior distribution is not trapped in any local mode and that the algorithm mixes efficiently, since every chain constantly tries to transmit information about its state to the others. EMC extents this idea, encompassing the positive features of PT and genetic algorithms inside a MCMC scheme.

Since $\beta$ and $\sigma^{2}$ are integrated out, only two parameters need to be sampled, namely the latent binary vector and the variable selection coefficient. In this set-up the full conditionals to be considered are

where $L$ is the number of chains in the the population and $t_{l}$, $1=t_{1}<t_{2}<\cdots<t_{L}$, is the temperature attached to the $l$th chain while the population $\bm{\gamma}$ corresponds to a set of chains that are retained simultaneously. Conditions for convergence of EMC algorithms are well understood and illustrated for instance in Jasra et al. (2007).

At each sweep of our algorithm, first the population $\bm{\gamma}$ in (10) is updated using a variety of moves inspired by genetic algorithms: “local moves”, the ordinary Metropolis-Hastings or Gibbs update on every chain; and “global moves” that include: i) selection of the chains to swap, based on some probabilistic measures of distance between them; ii) crossover operator, i.e. partial swap of the current state between different chains; iii) exchange operator, full state swap between chains. Both local and global moves are important although global moves are crucial because they allow the algorithm to jump from one local mode to another. At the end of the update of $\bm{\gamma}$, $\tau$ is then sampled using (11).

The implementation of EMC that we propose in this paper includes several novel aspects: the use of a wide range of moves including two new ones, a local move, based on the Fast Scan Metropolis-Hastings sampler, particularly suitable when $p$ is large and a bold global move that exploits the pattern of correlation of the predictors. Moreover, we developed an efficient scheme for tuning the temperature placement that capitalises the effective interchange between the chains. Another new feature is to use a Metropolis-within-Gibbs with adaptive proposal for updating $\tau$, as the full conditional (11) is not available in closed form.

In what follows, we will only sketch the rationale behind all the moves that we found useful to implement and discuss further the benefits of the new specific moves in Section 4.1. For the “large $p$, small $n$” paradigm and complex predictor spaces, we believe that using a wide portfolio of moves is needed and offers better guarantee of mixing.

From a notational point of view, we will use the double indexing $\gamma_{l,j}$, $l=1,\ldots,L$ and $j=1,\ldots,p$ to denote the $j$th latent binary indicator in the $l$th chain. Moreover we indicate by $\gamma_{l}=\left(\gamma_{l,1},\ldots,\gamma_{l,p}\right)^{T}$ the vector of binary indicators that characterise the state of the $l$th chain of the population $\bm{\gamma}=\left(\gamma_{1},\ldots,\gamma_{L}\right)$.

Given $\tau$, we first tried the simple MC3 idea of Madigan and York (1995), also used by Brown et al. (1998) where add/delete and swap moves are used to update the latent binary vector $\gamma_{l}$. For an add/delete move, one of the $p$ variables is selected at random and if the latent binary value is the proposed new value is $1$ or vice versa. However, when $p\gg p_{\gamma_{l}}$, where $p_{\gamma_{l}}$ is the size of the current model for the $l$th chain, the number of sweeps required to select by chance a binary indicator with a value of $1$ follows a geometric distribution with probability $p_{\gamma}/p$ which is much smaller than $1-p_{\gamma}/p$ to select a binary indicator with a value of . Hence, the algorithm spends most of the time trying to add rather than delete a variable. Note that this problem also affects RJ-type algorithms (Dellaportas et al., 2002). On the other hand, Gibbs sampling (George and McCulloch, G&McC hereafter, 1993) is not affected by this issue since the state of the $l$th chain is updated by sampling from

where $\gamma_{l,j^{-}}$ indicates for the $l$th chain all the variables, but the $j$th, $j=1,\ldots,p$ and $\gamma_{l,j}^{\left(1\right)}=\left(\gamma_{l,1},\ldots,\gamma_{l,j-1},\gamma_{l,j}=1,\gamma_{l,j+1},\ldots,\gamma_{l,p}\right)^{T}$. The main problem related to Gibbs sampling is the large number of models it evaluates if a full Gibbs cycle or any permutation of the indices is implemented at each sweep. Each model requires the direct evaluation, or at least the update, of the time consuming quantity $S\left(\gamma\right)$, equation (7) or (9), making practically impossible to rely solely on the Gibbs sampler when $p$ is very large. However, as sharply noticed by Kohn et al. (2001), it is wasteful to evaluate all the $p$ updates in a cycle because if $p_{\gamma_{l}}$ is much smaller than $p$ and given $\gamma_{l,j}=0$, it is likely that the sampled value of $\gamma_{l,j}$ is again .

Global move: crossover operator

The first step of this move consists of selecting the pair of chains $\left(l,r\right)$ to be operated on. We firstly compute a probability equal to the weight of the “Boltzmann probability”, $p_{t}\left(\gamma_{l}\left|\tau\right.\right)=\exp\left\{f\left(\gamma_{l}\left|\tau\right.\right)/t\right\}/F_{t}$, where $f\left(\gamma_{l}\left|\tau\right.\right)=\log p\left(\gamma_{l}\left|y,\tau\right.\right)+\log p\left(\gamma_{l}\right)$ is the log transformation of the full conditional (10) assuming $t_{l}=1$ $\forall l$, $l=1,\ldots,L$, and $F_{t}=\sum_{l=1}^{L}\exp\left\{f\left(\gamma_{l}\left|\tau\right.\right)/t\right\}$ for some specific temperature $t$, and then rank all the chains according to this. We use normalised Boltzmann weights to increase the chance that the two selected chains will give rise, after the crossover, to a new configuration of the population with higher posterior probability. We refer to this first step as “selection operator”.

Suppose that two new latent binary vectors are then generated from the selected chains according to some crossover operator described below. The new proposed population of chains $\bm{\gamma}^{\prime}=\left(\gamma_{1},\ldots,\gamma_{l}^{\prime},\ldots,\gamma_{r}^{\prime},\ldots,\gamma_{L}\right)$ is accepted with probability

where $Q_{t}\left(\bm{\gamma}\rightarrow\bm{\gamma}^{\prime}\left|\tau\right.\right)$ is the proposal probability, see Liang and Wong (2000).

In the following we will assume that four different crossover operators are selected at random at every EMC sweep: $1$-point crossover, uniform crossover, adaptive crossover (Liang and Wong, 2000) and a novel block crossover. Of these four moves, the uniform crossover which “shuffles” the binary indicators along all the selected chains is expected to have a low acceptance, but to be able to genuinely traverse regions of low posterior probability. The block crossover essentially tries to swap a group of variables that are highly correlated and can be seen as a multi-points crossover whose crossover points are not random but defined from the correlation structure of the covariates. In practice the block crossover is defined as follows: one variable is selected at random with probability $1/p$, then the pairwise correlation $\rho\left(X_{j},X_{j^{\prime}}\right)$ between the $j$th selected predictor and each of the remaining covariates, $j^{\prime}=1,\ldots,p$, $j^{\prime}\neq j$, is calculated. We then retain for the block crossover all the covariates with positive (negative) pairwise correlation with $X_{j}$ such that $\left|\rho\left(X_{j},X_{j^{\prime}}\right)\right|\geq\rho_{0}$. The threshold $\rho_{0}$ is chosen with consideration to the specific problem, but we fixed it at $0.25$. Evaluation of block crossover and comparisons with other crossover operators are presented on a real data example in Section 4.1.

Global move: exchange operator

The exchange operator can be seen as an extreme case of crossover operator, where the first proposed chain receives the whole second chain state $\gamma_{l}^{\prime}=\gamma_{r}$, and vice versa. In order to achieve a good acceptance rate, the exchange operator is usually applied on adjacent chains in the temperature ladder, which limits its capacity for mixing. To obtain better mixing, we implemented two different approaches: the first one is based on Jasra et al. (2007) and the related idea of delayed rejection (Green and Mira, 2001); the second, a bolder “all-exchange” move, is based on a precalculation of all the $L\left(L-1\right)/2$ exchange acceptance rates between all chains pairs (Calvo, 2005). Full relevant details are presented in Appendix A.1. Both of these bold moves perform well in the real data applications, see Section 4.1, and simulated examples, see Section 4.2, thus contributing to the efficiency of the algorithm.

Temperature placement

As noted by Goswami and Liu (2007), the placement of the temperature ladder is the most important ingredient in population based MCMC methods. We propose a procedure for the temperature placement which has the advantage of simplicity while preserving good accuracy. First of all, we fix the size $L$ of the population. In doing this, we are guided by several considerations: the complexity of the problem, i.e. $E\left(p_{\gamma}\right)$, the size of the data and computational limits. We have experimented and we recommend to fix $L\geq 3$. Even though some of the simulated examples had $p_{\gamma}\simeq 20$ (Section 4.2), we found that $L=5$ was sufficient to obtain good results. In our real data examples (Section 4.1), we used $L=4$ guided by some prior knowledge on $E\left(p_{\gamma}\right)$. Secondly, we fix at an initial stage, a temperature ladder according to a geometric scale such that $t_{l+1}/t_{l}=b$, $b>1$, $l=1,\ldots,L$ with $b$ relatively large, for instance $b=4$. To subsequently tune the temperature ladder, we then adopt a strategy based on monitoring only the acceptance rate of the delayed rejection exchange operator towards a target of $0.5$. Details of the implementation are left in Appendix A.1

2 Adaptive Metropolis-within-Gibbs for τ𝜏\tau

Various strategies can be used to avoid having to sample from the posterior distribution of the variable selection coefficient $\tau$. The easiest way is to integrate it out through a Laplace approximation (Tierney and Kadane, 1986) or using a numerical integration such as quadrature on an infinite interval. We do not pursue these strategies and the reasons can be summarised as follows. Integrating out $\tau$ in the population implicitly assumes that every chain has its own value of the variable selection coefficient $\tau_{l}$ (and of the latent binary vector $\gamma_{l}$). In this set-up, two unpleasant situations can arise: firstly, if a Laplace approximation is applied, equilibrium in the product space is difficult to reach because the posterior distribution of $\gamma_{l}$ depends, through the marginal likelihood obtained using the Laplace approximation, on the chain specific value of the posterior mode for $\tau_{l}$, $\hat{\tau}_{\gamma_{l}}$ (details in Appendix A.2). Since the strength of $X_{\gamma_{l}}$ to predict the response is weakened for chains attached to high temperatures, it turns out that for these chains, $\hat{\tau}_{\gamma_{l}}$ is likely to be close to zero. When the variable selection coefficient is very small, the marginal likelihood dependence on $X_{\gamma_{l}}$ decreases even further, see for instance (7), and chains attached to high temperatures will experience a very unstable behaviour, making the convergence in the product space hard to reach. In addition, if an automatic tuning of temperature ladder is applied, chains will increasingly be placed at a closer distance in the temperature ladder to balance the low acceptance rate of the global moves, negating the purpose of EMC.

In this paper the convergence is reached instead in the product space $\prod\nolimits_{l=1}^{L}\left[p\left(\gamma_{l}\left|y,\tau\right.\right)\right]^{1/t_{l}}p\left(\tau\right)$, i.e. the whole population is conditioned on a value of $\tau$ common to all chains. This strategy will alleviate the problems highlighted before allowing for faster convergence and better mixing among the chains. The procedure just described comes with an extra cost, i.e. sampling the value of $\tau$. However, this step is inexpensive in relation to the cost required to sample $\gamma_{l}$, $l=1,\ldots,L$. There are several strategies that can be used to sample $\tau$ from (11). We found useful to apply the idea of adaptive Metropolis-within-Gibbs described in Roberts and Rosenthal (2008). Conditions for the asymptotic convergence and ergodicity are guaranteed as we enforce the diminishing adaptive condition, i.e. the transition kernel stabilises as the number of sweeps goes to infinity and the bounded convergence condition, i.e. the convergence time of the kernel is bounded in probability. In our set-up using an adaptive proposal to sample $\tau$ has several benefits; amongst others it avoids the known problems faced by the Gibbs sampler when the prior is proper, but relatively flat (Natarajan and McCulloch, 1998) as can happen for Z-S priors when $n$ is large or for the independent case considered by Bae and Mallick (2004). Moreover, given an upper limit on the number of sweeps, the adaptation guarantees a better exploration of the tails of $p\left(\tau\left|y\right.\right)$ than with a fixed proposal. For details of the implementation and discussion of conditions for convergence, see Appendix A.2.

3 ESS algorithm

In the following, we refer to our proposed algorithm, Evolutionary Stochastic Search as ESS. If $g$-priors are chosen the algorithm is denoted as ESS$g$, while we use ESS$i$ if independent priors are selected (the same notation is used when $\tau$ is fixed or given a prior distribution). Without loss of generality, we assume that the response vector and the design matrix have both been centred and, in the case of independent priors, that the design matrix is also rescaled. Based on the two full conditionals (10) and (11) and the local and global moves introduced earlier, our ESS algorithm can be summarised as follows.

Given $\tau$, sample the population’s states $\bm{\gamma}$ from the two steps:

With probability $0.5$ perform local move and with probability $0.5$ apply at random one of the four crossover operators: $1$-point, uniform, block and adaptive crossover. If local move is selected, use FSMH sampling scheme independently for each chain (see Appendix A.1). Moreover every $100$ sweeps apply on the first chain a complete scan by a Gibbs sampler.

Perform the delayed rejection exchange operator or the all-exchange operator with equal probability. During the burn-in, only select the delayed rejection exchange operator.

When $\tau$ is not fixed but has a prior $p\left(\tau\right)$, given the latent binary configuration $\bm{\gamma}=\left(\gamma_{1},\ldots,\gamma_{L}\right)$, sample $\tau$ from an adaptive Metropolis-within-Gibbs sampling (Section 3.2).

From a computational point of view, we used the same fast form for updating $S\left(\gamma\right)$ as Brown et al. (1998), based on the QR-decomposition. Besides its numerical benefits, QR- decomposition can deal with the case $p_{\gamma}\geq n$. This avoids having to restrict the search to models with $p_{\gamma}<n$, and helps mixing during the burn-in phase.

Performance of ESS

The first real data example is an application of linear regression to investigate genetic regulation. To discover the genetic causes of variation in the expression (i.e. transcription) of genes, gene expression data are treated as a quantitative phenotype while genotype data (SNPs) are used as predictors, a type of analysis known as expression Quantitative Trait Loci (eQTL).

Here we focus on the ability of ESS to find a parsimonious set of predictors in an animal data set (Hubner et al., 2005), where the number of observations, $n=29$, is small with respect to the number of covariates $p=1,421$. This situation, where $n\ll p$, is quite common in animal experiments since environmental sources of variation are controlled as well as the biological diversity of the sample. For illustration, we report the analysis of one gene expression response, where we apply ESS$g$ with and without the hyperprior on $\tau$, see Table 1– eQTL. In the former case, thanks to the adaptive proposal, the Markov chain for $\tau$ mixes very well reaching an overall acceptance rate which is close to the target value $0.44$. Convergence issue is not a problem since the trace of the proposal’s standard deviation stabilises quickly and well inside the bounded conditions, see Figure 3.

In both cases a good mixing among the $L=4$ chains is obtained (Figure 1, top panels, ESS$g$ with $\tau=29$). Although in the case depicted in Figure 1 with fixed $\tau$, the convergence is reached in the product space $\prod\nolimits_{l=1}^{L}\left[p\left(\gamma_{l}\left|y\right.\right)\right]^{1/t_{l}}$, by visual inspection we see that each chain marginally reaches its equilibrium with respect to the others; moreover, thanks to the automatic tuning of the temperature placement during the burn-in, the distributions of the chains log posterior probabilities overlap nicely, allowing effective exchange of information between the chains. Table 1–eQTL, confirms that the automatic temperature selection works well (with and without the hyperprior on $\tau$) reaching an acceptance rate for the monitored exchange (delayed rejection) operator close to the selected target of $0.50$. The all-exchange operator shows a higher acceptance rate, while, in contrast to Jasra et al. (2007), the overall crossover acceptance rate is reasonable high: in our experience the good performance of the crossover operator is both related to the selection operator (Section 3.1) and the new block crossover which shows an acceptance rate far higher than the others. Finally the computational time on the same desktop computer (see details in Appendix B.3) is rather similar with or without the hyperprior $\tau$, $28$ and $30$ minutes respectively for $25,000$ sweeps with $5,000$ as burn-in.

The main difference among the two implementations of ESS$g$ is related to the posterior model size: when $\tau$ is fixed at $\tau=29$ (Unit Information Prior, Fernández et al., 2001), there is more uncertainty and support for larger models, see Figure 2 (a). In both cases we fix $E\left(p_{\gamma}\right)=4$ and $V\left(p_{\gamma}\right)=2$, following prior biological knowledge on the genetic regulation. The posterior mean of the variable selection coefficient is a little smaller than the Unit Information Prior, with ESS$g$ coupled with the Z-S prior favouring smaller models than when $\tau$ is set equal to $29$. The best model visited (and the corresponding $R_{\gamma}^{2}=1-S(\gamma)/y^{T}y$) is the same for both version of ESS$g$, while, when a hyperprior on $\tau$ is implemented, the “stability index” which indicates how much the algorithm persists on the first chain top $1,000$ (not unique) visited models ranked by the posterior probability (Appendix B.3), shows a higher stability, see Table 1– eQTL. In this case, having a data-driven level of shrinkage helps the search algorithm to better discriminate among competing models.

Our second example is related to the application of model (1) in another genomics example: $10,000$ SNPs, selected genome-wide from a candidate gene study, are used to predict the variation of Mass Spectography metabolomics data in a small human population, an example of a so-called mQTL experiment. A suitable dimension reduction of the data is performed to divide the spectra in regions or bins and $\log_{10}$-transformation is applied in order to normalise the signal.

We present the key findings related to a particular metabolite bin, but the same comments can be extended to the analysis of the whole data set, where we regressed every metabolites bin versus the genotype data ($n=50$ and $p=10,000$). In this very challenging case, we still found an efficient mixing of the chains (see Table 1–mQTL). Note that in this case the posterior mean of $\tau$, $63.577$, is a little larger than the Unit Information Prior, $\tau=n$, although the influence of the hyperprior is less important than in the previous real data example, see Figure 2 (b). In both examples, the posterior model size favours clearly polygenic control with significant support for up to four genetic control points (Figure 2) highlighting the advantage of performing multivariate analysis in genomics rather than the traditional univariate analysis.

As expected in view of the very large number of predictors, in the mQTL example the computational time is quite large, around $5$ hours for $20,000$ sweeps after a burn-in of $5,000$, but as shown in Table 1 by the “stability index” ($\approx 0$), we believe that the number of iterations chosen exceeds what is required in order to visit faithfully the model space. For such large data analysis tasks, parallelisation of the code could provide big gains of computer time and would be ideally suited to our multiple chains approach.

[Table 1 about here – Figure 1 about here – Figure 2 about here – Figure 3 about here]

We also evaluate the superiority of our ESS algorithm, and in particular the FSMH scheme and the block crossover, with respect to more traditional EMC implementations illustrated for instance in Liang and Wong (2000). Albeit we believe that using a wide portfolio of different moves enables any searching algorithm to better explore complicated model spaces, we reanalysed the first real data example, eQTL analysis, comparing: (i) ESS$g$ with only FSMH as local move vs ESS with only MC3 as local move; (ii) ESS$g$ with only block crossover vs ESS$g$ with only 1-point, only uniform and only adaptive crossover respectively. To avoid dependency of the results on the initialisation of the algorithm, we replicated the analysis $25$ times. Moreover, to make the comparison fair, in experiment (i) we run the two versions of ESS$g$ for a different number of sweeps ($25,000$ and $350,000$ with $5,000$ and $70,000$ as burn-in respectively), but matching the number of models evaluated. Results are presented in Table 2. We report here the main findings:

over the $25$ runs, ESS$g$ with FSMH reaches the same top visited model $68$% (17/25) while ESS$g$ with MC3 the same top model only $28$%, with a fixed $\tau$, and $88$% and $40$% respectively with Z-S prior. This ability is extended to the top models ranked by the posterior probability, data not shown, providing indirect evidence that the proposed new move helps the algorithm to increase its predictive power. The great superiority when FSMH scheme are implemented can be explained by comparing subplot (a) and (c) in Figure 1: the exchange of information between chains for ESS$g$ with MC3 as local move when $p>n$ (and $p\gg p_{\gamma}$) is rather poor, negating the purpose of EMC. ESS$g$ with MC3 has more difficulties to reach convergence in the product space and, in contrast to ESS$g$ with FSMH, the retained chain does not easily escape from local modes. This later point can be seen looking at Figure 1 (d) which magnifies the right hand tail of the kernel density of $\log p\left(\gamma\left|y\right.\right)$ for the recorded chain, pulling together the $25$ runs: interestingly ESS$g$ with FSMH is less “bumpy”, showing a better ability to escape from local modes and to explore more efficiently the right tail.

Regarding the second comparison when $\tau$ is fixed, ESS$g$ with only block crossover beats constantly the other crossover operators, with $80$% vs about $60$%, in terms of best model visited (Table 2) and models with higher posterior probability (data not shown), has higher acceptance rate (Table 3), showing also a great capacity to accumulate posterior mass as illustrated in Figure 4. The specific benefit of the block crossover is less pronounced when a prior on $\tau$ is specified, but we have already noticed that in this case having a hyperprior on $\tau$ greatly improves the efficiency of the search.

[Table 2 about here – Table 3 about here – Figure 4 about here]

2 Simulation study

We briefly report on a comprehensive study of the performance of ESS in a variety of simulated examples as well as a comparison with SSS. To make comparison with SSS fair, we use ESS${i}$, the version of our algorithm which assumes independent priors, $\Sigma_{\gamma}=\tau I_{p_{\gamma}}$,with $\tau$ fixed at $1$. Details of the simulated examples (6 set-ups) and how we conducted the simulation experiment (25 replication of each set-up) are given in Appendix B. The rationale behind the construction of the examples was to benchmark our algorithm against both $n>p$ and $p>n$ cases, to use as building blocks intricate correlation structures that had been used in previous comparisons by G&McC (1993, 1997) and N&G (2004), as well as a realistic correlation structure derived from genetic data, and to include elements of model uncertainty in some of the examples by using a range of values of regression coefficients.

In our example we observe an effective exchange of information between the chains (reported in Table 4) which shows good overall acceptance rates for the collection of moves that we have implemented. The dimension of the problem does not seem to affect the acceptance rates in Table 4, remarkably since values of $p$ range from $60$ to $1,000$ between the examples. We also studied specifically the performance of the global moves (Table 5) to scrutinise our temperature tuning and confirmed the good performance of ESS$i$ with good frequencies of swapping (not far from the case where adjacent chains are selected to swap at random with equal probability) and good measures of overlap between chains.

All the examples were run in parallel with ESS${i}$ and SSS 2.0 (Hans et al., 2007) for the same number of sweeps (22,000) and matching hyperparameters on the model size. Comparison were made with respect to the marginal probability of inclusion as well as the ability to reach models with high posterior probability and to persist in this region. For a detailed discussion of all comparison, see Appendix B.3.

Overall the covariates with non-zero effects have high marginal posterior probability of inclusion for ESS$i$ in all the examples, see Figure 6. There is good agreement between the two algorithms in general, with additional evidence on some examples (Figure 6 (c) and (d)) that ESS$i$ is able to explore more fully the model space and in particular to find small effects, leading to a posterior model size that is close to the true one. Measures of goodness of fit and stability, Table 6, are in good agreement between ESS$i$ and SSS. The comparison highlight that a key feature of SSS, its ability to move quickly towards the right model and to persist on it, is accompanied by a drawback in having difficulty to explore far apart models with competing explanatory power, in contrast to ESS$i$ (contaminated example set-up). Altogether ESS$i$ shows a small improvement of $R_{\gamma}^{2}$, related to its ability to pick up some of the small effects that are missed by SSS. Finally ESS$i$ shows a remarkable superiority in terms of computational time, especially when the simulated (and estimated) $p_{\gamma}$ is large. Altogether our comparisons show that we have designed a fully Bayesian MCMC-EMC sampler which is competitive with the effective search provided by SSS$i$.

In the same spirit of the real data example analysis, we also evaluate the superiority of the FSMH scheme with respect to more traditional EMC implementations, i.e when a MC3 local move is selected. While both versions of the search algorithm visit almost the same top models ranked by the posterior probability, ESS persists more on the top models.

[Table 4 about here – Table 5 about here – Table 6 about here

Figure 5 about here – Figure 6 about here]

Discussion

The key idea in constructing an effective MCMC sampler for $\gamma$ and $\tau$ is to add an extra parameter, the temperature, that weakens the likelihood contribution and enables escaping from local modes. Running parallel chains at different temperature is, on the other hand, expensive and the added computational cost has to be balanced against the gains arising from the various “exchanges” between the chains. This is why we focussed on developing a good strategy for selecting the pairs of chains, using both marginal and joint information between the chains, attempting bold and more conservative exchanges. Combining this with an automatic choice of the temperature ladder during burn-in is one of the key element of our ESS algorithm. Using PT in this way has the potential to be effective in a wide range of situations where the posterior space is multimodal.

To tackle the case where $p$ is large with respect to $p_{\gamma}$, the second important element in our algorithm is the use of a Metropolised Gibbs sampling-like step performed on a subset of indices in the local updating of the latent binary vector, rather than an MC3 or RJ-like updating move. The new Fast Scan Metropolis Hastings sampler that we propose to perform these local moves achieves an effective compromise between full Gibbs sampling that is not feasible at every sweep when $p$ is large and vanilla add/delete moves. Comparison of FSMH vs MC3 scheme on a real data example and simulation study shows the superiority of our new local move.

When a model with a prior on the variable selection coefficient $\tau$ is preferred, the updating of $\tau$ itself present no particular difficulties and is computationally inexpensive. Moreover, using an adaptive sampler makes the algorithm self contained without any time consuming tuning of the proposal variance. This latter strategy works perfectly well both in the $g$-prior and independent prior case as illustrated in Sections 4.1 and 4.2. Our current implementation does not make use of the output of the heated chains for posterior inference. Whether gains in variance reduction could be achieved in the spirit of Gramacy et al. (2007) is an area for further exploration, which is beyond the scope of the present work.

Our approach has been applied so far to linear regression with univariate response $y$. An interesting generalisation is that of a multidimensional $n\times q$ response $Y$ and the identification of regressors that jointly predict the $Y$ (Brown et al., 1998). Much of our set-up and algorithm carries through without difficulties and we have already implemented our algorithm in this framework in a challenging case study in genomics with multidimensional outcomes.

Acknowledgements

The authors are thankful to Norbert Hubner and Timothy Aitman for providing the data of the eQTL example, Gareth Roberts and Jeffrey Rosenthal for helpful discussions about adaptation and Michail Papathomas for his detailed comments. Sylvia Richardson acknowledges support from the MRC grant GO.600609.

Appendix

Appendix A Technical details of EMC implementation

In this Section we will describe some technical details omitted from the paper and related to the sampling schemes we used for the population of binary latent vectors $\bm{\gamma}$ and the selection coefficient $\tau$.

Let $\gamma_{l,j}$, $l=1,\ldots,L$ and $j=1,\ldots,p$ to denote the $j$th latent binary indicator in the $l$th chain. As in Kohn et al. (2001), let $\gamma_{l,j}^{\left(1\right)}=\left(\gamma_{l,1},\ldots,\gamma_{l,j-1},\gamma_{l,j}=1,\gamma_{l,j+1},\ldots,\gamma_{l,p}\right)^{T}$ and $\gamma_{l,j}^{\left(0\right)}=\left(\gamma_{l,1},\ldots,\gamma_{l,j-1},\gamma_{l,j}=0,\gamma_{l,j+1},\ldots,\gamma_{l,p}\right)^{T}$. Furthermore let $L_{l,j}^{\left(1\right)}\propto p\left(y\left|\gamma_{l,j}^{\left(1\right)},\tau\right.\right)$ and $L_{l,j}^{\left(0\right)}\propto p\left(y\left|\gamma_{l,j}^{\left(0\right)},\tau\right.\right)$ and finally $\theta_{l,j}^{\left(1\right)}=p\left(\gamma_{l,j}=1\left|\gamma_{l,j^{-}}\right.\right)$ and $\theta_{l,j}^{\left(0\right)}=1-\theta_{l,j}^{\left(1\right)}$. From (6) it is easy to prove that

where $p_{\gamma_{l}}$ is the current model size for the $l$th chain. Using the above equation, for $\gamma_{l,j}=1$ the normalised version of (12) can be written as

where $S\left(1/t_{l}\right)=\left.\theta_{l,j}^{\left(1\right)}\right.^{1/t_{l}}\left.L_{l,j}^{\left(1\right)}\right.^{1/t_{l}}+\left.\theta_{l,j}^{\left(0\right)}\right.^{1/t_{l}}\left.L_{l,j}^{\left(0\right)}\right.^{1/t_{l}}$ with $\left[p\left(\gamma_{l,j}=1\left|y,\gamma_{l,j^{-}},\tau\right.\right)\right]^{1/t_{l}}$ defined similarly. Hence if $\left.\theta_{l,j}^{\left(1\right)}\right.^{1/t_{l}}$ is very small, then $\left[p\left(\gamma_{l,j}=1\left|y,\gamma_{l,j^{-}},\tau\right.\right)\right]^{1/t_{l}}$ is small as well. Therefore for the Gibbs sampler with a beta-binomial prior on the model space, the posterior probability of $\gamma_{l,j}=1$ depends crucially on $\left.\theta_{l,j}^{\left(1\right)}\right.^{1/t_{l}}$.

In the following we derive a Fast Scan Metropolis-Hastings scheme specialised for Evolutionary Monte Carlo or parallel tempering. We define $Q\left(1\rightarrow 0\right)=Q\left(\gamma_{l,j}^{\left(1\right)}\rightarrow\gamma_{l,j}^{\left(0\right)}\right)$ as the proposal probability to go from $1$ to and $Q\left(0\rightarrow 1\right)$ the proposal probability to go from to $1$ for the $j$th variable and $l$th chain. Moreover using the notation introduced before, the Metropolis-within-Gibbs version of (12) to go from to $1$ in the EMC local move is

with a similar expression for $\alpha_{l}^{\text{MwG}}\left(1\rightarrow 0\right)$. The proof of the Propositions are omitted since they are easy to check. We first introduce the following Proposition which is useful for the calculation of the acceptance probability in the FSMH scheme.

The FSMH scheme can be seen as a random scan Metropolis-within-Gibbs algorithm where the number of evaluations is linked to the prior/current model size and the temperature attached to the chain. The computation requirement for the additional acceptance/rejection step is very modest since the normalised tempered version of (A.1) is used.

Finally it can be proved that the Gibbs sampler is more efficient than the FSMH scheme, i.e. for a fixed number of iterations, Gibbs sampling MCMC standard error is lower than for FSMH sampler. However the Gibbs sampler is computationally more expensive so that, if $p$ is very large, as described in Kohn et al. (2001), FSMH scheme becomes more efficient per floating point operation.

Global move: exchange operator

The exchange operator can be seen as an extreme case of crossover operator, where the first proposed chain receives the whole second chain state $\gamma_{l}^{\prime}=\gamma_{r}$, and the second proposed chain receives the whole first state chain $\gamma_{r}^{\prime}=\gamma_{l}$, respectively.

In order to achieve a good acceptance rate, the exchange operator is usually applied on adjacent chains in the temperature ladder, which limits its capacity for mixing. To obtain better mixing, we implemented two different approaches: the first one is based on Jasra et al. (2007) and the related idea of delayed rejection (Green and Mira, 2001); the second one on Gibbs distribution over all possible chains pairs (Calvo, 2005).

The delayed rejection exchange operator tries first to swap the state of the chains that are usually far apart in the temperature ladder, but, once the proposed move has been rejected, it performs a more traditional (uniform) adjacent pair selection, increasing the overall mixing between chains on one hand without drastically reducing the acceptance rate on the other. However its flexibility comes at some extra computational costs and in particular the additional evaluation of the pseudo move necessary to maintain detailed balance (Green and Mira, 2001). Details are reported below.

Suppose two chains are selected at random, $l$ and $r$ with $l\neq r$, in order to swap their binary latent vector. Then, given that $\gamma_{l}^{\prime}=\gamma_{r}$, $\gamma_{r}^{\prime}=\gamma_{l}$ and $Q_{t}\left(\bm{\gamma}\rightarrow\bm{\gamma}^{\prime}\right)=Q_{t}\left(\bm{\gamma}^{\prime}\rightarrow\bm{\gamma}\right)$, (13) reduces to

Since the two chains are selected at random, the above acceptance probability decreases exponentially with the difference $\left|1/t_{l}-1/t_{r}\right|$ and therefore most of the proposed moves are rejected. If rejected, a delayed rejection-type move is applied between two random adjacent chains, with $l$ the first one and $s$, $\left|l-s\right|=1$, the second one, giving rise to the new acceptance probability

where the pseudo move $\bm{\gamma}^{\ast}$ is necessary in order to maintain the detailed balance condition (Green and Mira, 2001).

Temperature placement

A.2 Adaptive Metropolis-within-Gibbs for τ𝜏\tau

where $\hat{\lambda}$ is the posterior mode after the transformation $\lambda=\log\left(\tau\right)$, which is necessary to avoid problems on the boundary, $\sigma_{\hat{\lambda}}$ is the approximate squared root of the variance calculated in $\hat{\lambda}$ and $J\left(\cdot\right)$ is the Jacobian of the transformation. Details about Laplace approximation can be found in Tierney and Kadane (1986). Similar derivations when $p\left(\sigma^{2}\right)\propto\sigma^{-2}$ are presented in Liang et al. (2008). Finally throughout the presentation we will assume that $n>p_{\gamma}$ and that $a_{g}$ and $b_{g}$ are fixed small as in Kohn et al. (2001).

Cubic equation for Zellner-Siow priors If $p\left(\tau\right)=InvGa\left(a_{\tau},b_{\tau}\right)$ the posterior $\hat{\lambda}$ mode is the only positive root of the integrand function

where the last factor in the above equation $e^{\lambda}=\left|de^{\lambda}/d\lambda\right|$ is the Jacobian of the transformation. After the calculus of the first derivative of the log transformation and some algebra manipulations, it can be shown that $e^{\hat{\lambda}}$ is the solution of the cubic equation

where $c_{1}=\left(2a_{\sigma}+n-1-p_{\gamma}\right)/2$, $c_{2}=\left(2a_{\sigma}+n-1\right)/2$, $c_{3}=2b_{\sigma}+y^{T}y$ and $c_{4}=2b_{\sigma}+y^{T}y\left(1-R_{\gamma}^{2}\right)$. Following Liang et al. (2008), since $\lim_{\lambda\rightarrow-\infty}\partial I_{\lambda}/\partial\lambda>0$, because $c_{3}b_{\tau}>0$, and $\lim_{\lambda\rightarrow\infty}\partial I_{\lambda}/\partial\lambda<0$, because $\left(c_{1}-c_{2}-a_{\tau}\right)c_{4}<0$, at least one real positive solution exists. Moreover since $-\left(c_{3}b_{\tau}\right)/\left(c_{1}-c_{2}-a_{\tau}\right)c_{4}>0$, the remaining two real solutions should have the same sign (Abramowitz and Stegun, 1970). A necessary condition for the existence of just one real positive solution is that the summation of all the pairs-products of the coefficients is negative

and this happens if $b_{\tau}/a_{\tau}>c_{3}/\left(c_{3}+c_{4}\right)$. When $R_{\gamma}^{2}\rightarrow 0$ and thus $c_{3}=c_{4}$, the above condition corresponds to $b_{\tau}>a_{\tau}/2$ and when $R_{\gamma}^{2}\rightarrow 1$, as $c_{3}/\left(c_{3}+c_{4}\right)\approx 1$ especially when $y^{T}y$ is large, which might be expected when $n$ becomes large, the condition is equivalent to $b_{\tau}>a_{\tau}$. Therefore it turns out that a sufficient condition for the existence of just one real positive solution in (A.1) is $b_{\tau}>a_{\tau}$.

The positive semidefiniteness of the approximate variance can be proved as follows. First of all it is worth noticing that all the terms in (A.8) are of the same order $O_{p}\left(e^{-\lambda}\right)$. Then, when $R_{\gamma}^{2}\rightarrow 0$, the positive semidefiniteness is always guaranteed, while when $R_{\gamma}^{2}\rightarrow 1$, provided that $y^{T}y$ is large, the middle term in (A.8) tends to zero and the condition is fulfilled if $b_{\tau}>c_{1}$.

Quadratic equation for Liang et al. (2008) prior If $p\left(\tau\right)\propto\left(1+\tau\right)^{-c_{\tau}}$, with $c_{\tau}>0$, $e^{\hat{\lambda}}$ is only the positive root of the integrand function

or, after the first derivative of the log transformation, the solution of the quadratic equation

with $c_{1}^{\ast}=\left[2a_{\sigma}+n-1-\left(p_{\gamma}+2c_{\tau}\right)\right]/2$ and $c_{2}$, $c_{3}$ and $c_{4}$ defined as above. The discriminant of the quadratic equation is $\Delta=\left(c_{1}^{\ast}c_{3}-c_{2}c_{4}c_{3}+c_{3}+c_{4}\right)^{2}-4\left(c_{1}^{\ast}-c_{2}+1\right)c_{4}c_{3}$ which is always greater than zero and therefore two real roots exist. Since one of them is positive in order to prove that (A.9) admits just one positive solution, it is necessary to show that

which is true provided that $\left(c_{1}^{\ast}-c_{2}+1\right)c_{4}c_{3}<0$. Moreover the approximate variance can be written as

which is positive semidefinite when $R_{\gamma}^{2}\rightarrow 0$ if $c_{2}>c_{1}^{\ast}$, which is always verified, while, if $R_{\gamma}^{2}\rightarrow 1$ and $y^{T}y$ is large, equation (A.10) is not positive unless $p_{\gamma}+2c_{\tau}>2a_{\sigma}+n-1$.

The explicit solution of the posterior mode is also available

which corresponds to MLE if $c_{\tau}=0$.

Diminishing adaptive and bounded conditions

Appendix B Performance of ESS: Simulation study

In this Section we report in details on the performance of ESS in a variety of simulated examples. Main conclusions are summarised in the Section 4.2.

Firstly we analyse the simulated examples with ESS${i}$ the version of our algorithm which assumes independent priors, $\Sigma_{\gamma}=\tau I_{p_{\gamma}}$, so as to enable comparisons with SSS which also implements an independent prior. Moreover, in order to make to comparison with SSS fair, in the simulation study only the first step of the algorithm described in Section 3.3 is performed, with $\tau$ fixed at $1$. As in SSS, standardisation of the covariates is done before running ESS${i}$. We run ESS${i}$ and SSS 2.0 (Hans et al., 2007) for the same number of sweeps (22,000) and with matching hyperparameters on the model size.

Secondly, to discuss the mixing properties of ESS when a prior $p\left(\tau\right)$ is defined on $\tau$, we implement both the $g$-prior and independent prior set-up for a particular simulated experiment. To be precise in the former case we will use the Zellner-Siow priors (8), and for the latter we will specify a proper but diffuse exponential distribution as suggested by Bae and Mallick (2004).

We apply ESS with independent priors to an extensive and challenging range of simulated examples with $\tau$ fixed at $1$: the first three examples (Ex1-Ex3) consider the case $n>p$ while the remaining three (Ex4-Ex6) have $p>n$. Moreover in all examples, except the last one, we simulate the design matrix, creating more and more intricated correlation structures between the covariates in order to test the proposed algorithm in different and increasingly more realistic scenarios. In the last example, we use, as design matrix, a genetic region spanning $500$-kb from the HapMap project (Altshuler et al., 2005).

Simulated experiments Ex1-Ex5 share in common the way we build $X$. In order to create moderate to strong correlation, we found useful referring to two simulated examples in George and McCulloch, G&McC hereafter, (1993) and in G&McC (1997): throughout we call $X_{1}$ ($n\times 60$) and $X_{2}$ $(n\times 15)$ the design matrix obtained from these two examples. In particular the $j$th column of $X_{1}$, indicated as $X_{\left(1\right)j}$, is simulated as $X_{\left(1\right)j}=X_{j}^{\ast}+Z$, where $X_{1}^{\ast},\ldots,X_{60}^{\ast}$ iid $\sim N_{n}\left(0,1\right)$ independently form $Z\sim N_{n}\left(0,1\right)$, inducing a pairwise correlation of $0.5$. $X_{2}$ is generated as follows: firstly we simulated $Z_{1},\ldots,Z_{15}$ iid $\sim N_{n}\left(0,1\right)$ and we set $X_{\left(2\right)j}=Z_{i}+2Z_{j}$ for $j=1,3,5,8,9,10,12,13,14,15$ only. To induce strong multicollinearity, we then set $X_{\left(2\right)2}=X_{\left(2\right)1}+0.15Z_{2}$, $X_{\left(2\right)4}=X_{\left(2\right)3}+0.15Z_{4}$, $X_{\left(2\right)6}=X_{\left(2\right)5}+0.15Z_{6}$, $X_{\left(2\right)7}=X_{\left(2\right)8}+X_{\left(2\right)9}-X_{\left(2\right)10}+0.15Z_{7}$ and $X_{\left(2\right)11}=X_{\left(2\right)14}+X_{\left(2\right)15}-X_{\left(2\right)12}-X_{\left(2\right)13}+0.15Z_{11}$. A pairwise correlation of about 0.998 between $X_{\left(2\right)j}$ and $X_{\left(2\right)j+1}$ for $j=1,3,5$ is introduced and similarly strong linear relationship is present within the sets $\left(X_{\left(2\right)7},X_{\left(2\right)8},X_{\left(2\right)9},X_{\left(2\right)10}\right)$ and $\left(X_{\left(2\right)11},X_{\left(2\right)12},X_{\left(2\right)13},X_{\left(2\right)14},X_{\left(2\right)15}\right)$.

Then, as in Nott and Green, N&G hereafter, (2004) Example 2, more complex structures are created by placing side by side combinations of $X_{1}$ and/or $X_{2}$, with different sample size. We will vary the number of samples $n$ in $X_{1}$ and $X_{2}$ as we construct our examples. The levels of $\beta$ are taken from the simulation study of Fernández et al. (2001), while the number of true effects, $p_{\gamma}$, with the exception of Ex3, varies from $5$ to $16$. Finally the simulated error variance ranges from $0.05^{2}$ to $2.5^{2}$ in order to vary the level of difficulty for the search algorithm. Throughout we only list the non-zero $\beta_{\gamma}$ and assume that $\beta_{\gamma^{-}}=0^{T}$. The six examples can be summarised as follows:

$X=X_{1}$ is a matrix of dimension $120\times 60$, where the responses are simulated from (1) using $\alpha=0$, $\gamma=\left(21,37,46,53,54\right)^{T}$, $\beta_{\gamma}=\left(2.5,0.5,-1,1.5,0.5\right)^{T}$, and $\varepsilon\sim N\left(0,2^{2}I_{120}\right)$. In the following we will not refer to the intercept $\alpha$ any more since, as described in Section 3.3 in the paper, we consider $y$ centred and hence there is no difference in the results if the intercept is simulated or not. This is the simplest of our example, although, as reported in G&McC (1993) the average pairwise correlation is about $0.5$, making it already hard to analyse by standard stepwise methods.

This example is taken directly from N&G (2004), Example 2, who first introduce the idea of combining simpler “building blocks” to create a new matrix $X$ : in their example $X=\left[X_{2}^{\left(1\right)}X_{2}^{\left(2\right)}\right]$ is a $300\times 30$ matrix, where $X_{2}^{\left(1\right)}$ and $X_{2}^{\left(2\right)}$ are of dimension $300\times 15$ and have each the same structure as $X_{2}$. Moreover $\gamma=\left(1,3,5,7,8,11,12,13\right)^{T}$, $\beta_{\gamma}=\left(1.5,1.5,1.5,1.5,-1.5,1.5,1.5,1.5\right)^{T}$ and $\varepsilon\sim N\left(0,2.5^{2}I_{300}\right)$. We chose this example for two reasons: firstly, since the correlation structure in $X_{2}$ is very involved, we test the proposed algorithm under strong and complicated correlations between the covariates; secondly, since $y$ is not simulated from the second “block”, we are interested to see if the proposed algorithm does not select any variable that belongs to the second group.

As in G&McC (1993), Example 2, $X=X_{1}$, is a $120\times 60$ matrix, $\beta=\left(\beta_{1},\ldots,\beta_{60}\right)^{T}$, $\left(\beta_{1},\ldots,\beta_{15}\right)=\left(0,\ldots,0\right)$, $\left(\beta_{16},\ldots,\beta_{30}\right)=\left(1,\ldots,1\right)$, $\left(\beta_{31},\ldots,\beta_{45}\right)=\left(2,\ldots,2\right)$, $\left(\beta_{46},\ldots,\beta_{60}\right)=\left(3,\ldots,3\right)$ and $\varepsilon\sim N\left(0,2^{2}I_{120}\right)$. The motivation behind this example is to test the strength of the proposed algorithm to select a subset of variables which is large with respect to $p$ while preserving the ability not to choose any of the first $15$ variables.

The design matrix $X$, $120\times 300$, is constructed as follows: firstly we create a new $120\times 60$ “building block”, $X_{3}$, combining $X_{2}$ and a smaller version of $X_{1}$, $X_{1}^{\ast}$, a $120\times 45$ matrix simulated as $X_{1}$, such that $X_{3}=\left[X_{2}X_{1}^{\ast}\right]$ (dimension $120\times 60$). Secondly we place side by side five copies of $X_{3}$, $X=\left[X_{3}^{\left(1\right)}X_{3}^{\left(2\right)}X_{3}^{\left(3\right)}X_{3}^{\left(4\right)}X_{3}^{\left(5\right)}\right]$: the new design matrix alternates blocks of covariates of high and complicated correlation, as in G&McC (1997), with regions where the correlation is moderate as in G&McC (1993). We simulate the response selecting $16$ variables from $X$, $\gamma=\left(1,11,30,45,61,71,90,105,121,131,150,165,181,191,210,225\right)^{T}$ such that every pair belongs alternatively to $X_{2}$ or $X_{1}$. We simulate $y$ using $\beta_{\gamma}=\left(2,-1,1.5,1,0.5,2,-1,1.5,1,0.5,2,-1,-1,1.5,1,0.5\right)^{T}$ with $\varepsilon\sim N\left(0,2.5^{2}I_{120}\right)$. This example is challenging in view of the correlation structure, the number of covariates $p>n$ and the different levels of the effects.

This is the most challenging example that we simulated and it is based on the idea of contaminated models. The matrix $X$, $200\times 1000$, is $X=\left[X_{3}^{\left(1\right)}X_{3}^{\left(2\right)}X_{3}^{\left(3\right)}X_{1}^{\ast\ast}X_{3}^{\left(4\right)}X_{3}^{\left(5\right)}X_{3}^{\left(6\right)}X_{3}^{\left(7\right)}X_{3}^{\left(8\right)}\right]$, with $X_{1}^{\ast\ast}$, a $200\times 520$ larger version of $X_{1}$. We partitioned the responses such that $y=[y_{1}y_{2}]^{T}$: $y_{1}$ is simulated from “model 1” ($\gamma^{1}=\left(701,730,745,763,790,805,825,850,865,887\right)$ and $\beta_{\gamma}^{1}=\left(2,-1,1.5,1,0.5,2,-1,1.5,2,-1\right)$) while $y_{2}$ is simulated from “model 2” ($\gamma^{2}=\left(1,38,63,98,125\right)$ and $\beta_{\gamma}^{2}=\left(2,-1,1.5,1,0.5\right)$). Finally, fixing $\varepsilon\sim N\left(0,0.05^{2}I_{200}\right)$ and the sample size in the two models such that $y_{1}$ and $y_{2}$ are vectors of dimension $1\times 160$ and $1\times 40$ respectively, $y$ is retained if, given the sampling variability, we find $R_{\gamma^{1}}^{2}\geq 0.6$ and $R_{\gamma^{1}}^{2}/8\leq R_{\gamma^{2}}^{2}\leq R_{\gamma^{1}}^{2}/10$: in this way we know that “model 1” accounts for most of the variability of $y$, but without a negligible effect for “model 2”. In this example, we measure the ability of the proposed algorithm to recognise the most promising model and therefore being robust to contaminations. However since ESS can easily jump between local modes we are also interested to see if “model 2” is selected.

The last simulated example is based on phased genotype data from HapMap project (Altshuler et al., 2005), region ENm014, Yoruba population: the data set originally contained 1,218 SNPs (Single Nucleotide Polymorphism) for 120 chromosomes, but after eliminating redundant variables, the design matrix reduced to $120\times 775$. While in the previous examples a “block structure” of correlated variables is artificially constructed, in this example blocks of linkage disequilibrium (LD) derive naturally from genetic forces, with a slow decay of the level of pairwise correlation between SNPs. Finally we chose $\gamma=\left(50,75,140,200,300,400,500,650,700,770\right)^{T}$ such that the effects are visually inside blocks of LD, with their size simulated from $\beta_{\gamma}\sim N\left(0,3^{2}I_{10}\right)$ with $\varepsilon\sim N\left(0,0.10^{2}I_{120}\right)$. Since the simulated effects can range roughly between $\left(-6,6\right)$, this will allow us to test also the ability of ESS$i$ to select small effects.

We conclude this Section by reporting how we conducted the simulation experiment: every example from Ex1 to Ex6 has been replicated $25$ times and the results presented for example Ex1 to Ex5 are averaged over the $25$ replicates. For Ex6 the effects size change so average across replicated is only done for the mixing properties. ESS$i$ with $\tau$ =1 was applied to each example/sample, recording the visited sequence of $\gamma_{1}$ for $20,000$ sweeps after a burn-in of $2,000$ required for the automatic tuning of the temperature placement, Section 3.1 With the exception of Ex2 and Ex3, where we used an indifferent prior, $p\left(\gamma\right)=\left(1/2\right)^{p}$, we analysed the remaining examples setting $E\left(p_{\gamma}\right)=5$ with $V\left(p_{\gamma}\right)=E\left(p_{\gamma}\right)\left(1-E\left(p_{\gamma}\right)/p\right)$ which corresponds to a binomial prior over $p_{\gamma}$. In order to establish the sensitivity of the proposed algorithm to the choice of $E\left(p_{\gamma}\right)$ we also analysed Ex1 fixing $E\left(p_{\gamma}\right)=10$ and $20$. Moreover in all the examples we chose $L=5$ with the starting value of $\bm{\gamma}$ chosen at random. The remaining two hyperparameters to be fixed, namely $a_{\sigma}$ and $b_{\sigma}$, are set equal to $a_{\sigma}=10^{-6}$ and $b_{\sigma}=10^{-3}$ as in Kohn et al. (2001) which corresponds to a relative uninformative prior.