Despite the wide adoption of Bayesian methods in linear models, prior elicitation for linear models is still an open problem. The parameter prior distribution π γ can be expressed as π γ ( β γ , α , σ 2 ) = π γ ( β γ | α , σ 2 ) π γ ( α , σ 2 ) . A standard Jeffreys’ prior is often used for the intercept and error variance, which are often common to all models considered, i.e. π γ ( α , σ 2 ) ∝ σ − 2 . One of the most popular priors used for the model parameters β γ is Zellner’s g-prior [53], namely β γ | σ 2 , M γ ∼ N ( 0 , g σ 2 ( X γ T X γ ) − 1 ) .

This is popular because of its computational efficiency in evaluating marginal likelihoods and performing model search. It is also attractive because of its intuitive interpretation arising from analysis of a conceptual sample generated using the same design matrix X γ as in the data at hand. It is a special case of spike-and-lab family with the slab density given by the Normal density above and the spike being a point mass at zero. Also, g-priors are appealing in variable selection problems since they require the user to specify only value (or hyper-prior) for the scalar hyper-parameter g. This controls the prior variance of the parameters; the effective prior sample size is n / g. Several choices of g have been proposed [6, 13, 19, 23, 28, 32, 48, 54].

Based on an extensive simulation study, Porwal and Raftery [37] found that in comparing parameter priors for BMA, three adaptive g-priors performed the best among many popular choices across the statistical tasks of parameter estimation, interval estimation, model inference, point prediction and interval prediction. In what follows, we shall focus only on these three parameter prior choices, namely: •

g = n : First proposed by [13], it corresponds to a prior sample size equal to n and has been found to work well in high dimensional settings [52]. The complexity penalty for a model using this specification is effectively half that in the BIC [37].

In terms of theoretical properties, all three priors are model-selection consistent [32], except when the true model is the null model. This means that if the true model, denoted by M γ T belongs to the model space, then the posterior probability of the true model, P ( M γ = M γ T | y ) → 1 as the sample size n → ∞. None of the above priors suffers from Bartlett’s paradox [1]. BMA with g = n is subject to the “information paradox”, but it has been argued that information consistency is of little practical importance in real data applications [37]. The EB-local and Hyper-g BMA methods are not subject to the information paradox.

Model space priors require specification of the prior probabilities of all models M γ , indexed by the binary variable inclusion vector γ. A common approach is to consider the inclusion of each variable as an independent and exchangeable Bernoulli random variable with a common prior probability of inclusion θ, i.e. (2.1) p ( M γ | θ ) = ∏ i = 1 p θ γ i ( 1 − θ ) 1 − γ i = θ p γ ( 1 − θ ) p − p γ , where θ is the prior expected fraction of the β j ’s that are not zero and p γ = ∑ i = 1 p γ i is the total number of covariates included in the model M γ .

In the absence of prior information, a common choice is to set θ = 0.5. This induces a uniform prior over all models with p ( M γ ) = 2 − p for all γ ∈ { 0 , 1 } p , where p is the total number of covariates considered. The expected prior model size under the uniform model prior is p / 2. However, choosing θ = 0.5 does not provide any multiplicity control [19].

For a fixed value of θ, the above prior induces a binomial prior for model size S = ∑ i = 1 p γ i i.e. S ∼ B i n ( p , θ ), with prior mean p θ and prior variance p θ ( 1 − θ ). Another way to specify θ is by using the researcher’s prior belief about expected model size E [ S ]. Sala-i-Martin et al [46] (hereafter SDM) recommended a prior expected model size of 7 based on their growth regression analysis. Similar priors for expected model size have also been proposed elsewhere [45, 30]. Hence, we can define an SDM version of the above prior with θ S D M = 7 / p.

Any fixed choice of θ can lead to rather informative priors on model size p γ . One way to address this issue is by estimating θ from the data using an empirical Bayes (EB) approach [19]. The EB approach involves maximizing the marginal likelihood of the data given θ: (2.2) θ ˆ E B = arg max θ ∈ [ 0 , 1 ] ∑ γ P ( Y | M γ ) P ( M γ | θ ) .

However, maximization of (2.2) can be computationally challenging, especially when p is large since the sum is over all models. Moreover, when p is large, marginal likelihood evaluation for all models is not feasible and the sum is approximated based on the models explored by MCMC. To optimize the marginal likelihood in (2.2), we implement Algorithm 1, iterating between fitting the BMA approach to find likely models given θ and solving (2.2) using the fitted models to find a new θ:

An alternative way to reduce the sensitivity of the posterior distribution to prior assumptions is to use hierarchical modeling and specify a weak hyper-prior for θ. One choice of such a hyper-prior is a Beta distribution, θ ∼ Beta ( a , b ). Marginalizing out θ in (2.1), gives (2.3) P ( M γ | a , b ) = ∫ 0 1 p ( M γ | θ ) p ( θ ) d θ = B ( p γ + a , p − p γ + b ) B ( a , b ) , where B ( a ′ , b ′ ) = Γ ( a ′ ) Γ ( b ′ ) Γ ( a ′ + b ′ ) is the Beta function. It thus induces a Beta-Binomial ( a , b ) prior on the model size S with probability mass function P S ( s ) = p s B ( s + a , p − s + b ) B ( a , b ) .

Under a uniform prior on θ, i.e. when a = b = 1, (2.3) simplifies to p ( M γ ) = 1 p + 1 p p γ − 1 . This is a combination of a uniform prior over model size with a uniform prior over the models of same size given the model size.

Under a Beta-Binomial (BB) prior, the prior expected model size is E [ S ] = a a + b p. Similarly to a Bernoulli prior, we can elicit the prior in terms of the prior expected model size E [ S ]. To facilitate prior elicitation, we fix a = 1. We can then define an SDM version of the BB prior (BB-SDM) with an expected prior model size, such as E [ S ] = 7, by setting b S D M = p E [ S ] − 1. Note that SDM themselves [46] did not use a Beta-Binomial prior on models, but only a Bernouilli prior.

Alternatively, we can use an EB approach to learn b from the data. This can be done by maximizing the marginal likelihood given b, namely (2.4) b ˆ E B = arg max b ∈ ( 0 , ∞ ) ∑ γ P ( Y | M γ ) P ( M γ | a = 1 , b ) . We find the optimal value, b ˆ E B , using Algorithm 2.

For Zellner’s g-prior, we require the model size to be no larger than the number of regression coefficients that can be identified from the data, so that p γ < n − 2. Thus for higher dimensional datasets ( p > n ), we require that P ( M γ ) = 0 for all models with model size greater than n − 2. Hence, we use truncated versions of the priors defined in (2.1) and (2.3), namely (2.5) p ( M γ | θ ) ∝ θ p γ ( 1 − θ ) p − p γ 1 { p γ < n − 2 } , (2.6) P ( M γ | a , b ) ∝ B ( p γ + a , p − p γ + b ) B ( a , b ) 1 { p γ < n − 2 } .

Castillo et al [3] introduced complexity priors, also known in the literature as diffusing [35] or power priors [5]. Here the marginal probability of inclusion of any variable decays at the rate p − κ for some κ > 0, where p is the total number of possible covariates. This specifies a vanishing prior probability of large models and leads to a faster rate of rejection of spurious parameters, at the cost of slower rates of detection of active parameters [44]. Similar priors have also been used elsewhere [51, 35, 42].

The complexity prior is defined as p ( M γ ) ∝ p − κ | γ | 1 { | γ | ≤ s 0 } , where s 0 is a pre-specified integer specifying the maximum number of important covariates and | γ | is the model size. In the absence of external information, we set s 0 = min { n − 2 , p }. This prior is implemented in the BAS package as tr.power.prior(kappa,trunc). We implement Complexity priors with κ = 1 [43, 44], and with κ = 2 which is the default choice in the BAS package.

To illustrate the effect of different model space priors, we use two datasets from our analysis: Boston Housing ( n = 506 , p = 103 ) and Nutrimouse ( n = 40 , p = 120 ) (Figure 1). The solid lines show the independent Bernoulli prior from (2.1), while the dashed lines represent the Beta-Binomial prior in (2.3) and the dash-dotted lines illustrate Complexity priors. For the Nutrimouse dataset ( p > n ), we use the truncated versions as discussed above. The colors group different flavors of methods: (i) Uniform versions with θ = 0.5 or b = 1 (blue), (ii) SDM versions with expected prior model size 7 (red), (iii) EB versions with θ or b learned from the data (green), (iv) Complexity prior with κ = 1 (orange), and (v) Complexity prior with κ = 2 (brown).

The Bernoulli model space priors are very concentrated around their mean, p θ. The complexity priors are concentrated around smaller model sizes with a mode at 0. The BB priors, on the other hand, are more diffuse, implying more prior uncertainty about model size. For the Nutrimouse dataset, all the model space priors assign zero probability to any model with size greater than ( n − 2 ) = 38. Among the Bernoulli versions, θ = 0.5 implies a prior mode around min { p / 2 , n − 2 } while θ S D M has a prior model size of 7 (the same as the prior mean). The prior model size distribution induced by the Ber ( θ E B ) prior adapts based on the data, with a prior mode between θ = 0.5 and θ S D M for the Boston Housing dataset, while having the lowest prior mode among the Bernoulli priors considered for the Nutrimouse dataset. The BB ( 1 , 1 ) prior corresponds to a uniform prior over model size. The BB ( 1 , θ S D M ) and BB ( 1 , θ E B ) priors both induce a model size distribution with prior mode at zero.

We based our analysis on 14 publicly available datasets, of which six are available from UCI machine learning repository and the others are examples available in the literature. The sample size and number of candidate variables along with the data sources for the different datasets are listed in Table 2. These include the classical statistical setting with n > p, high dimensional datasets with p > n, and intermediate settings where n ≈ p. For each dataset, continuous covariates are standardized to have zero mean and variance 1 and the response variable is centered to have zero mean. Details of the choice of datasets and additional pre-processing can be found in [37].

For each dataset, we selected a data generating model that closely approximates the real dataset. We carry out all-subsets regression for datasets with p < 30 using the leaps package [33] in R. We then selected the largest model with all variables significant at the 0.05 level. For datasets with p > 30, all subsets regression is not computationally feasible. For these datasets, we obtained a filtered list of variables using iterative sure independence screening [12]. If the filtered list contained more than 30 variables, we selected the top 30 variables with the highest R 2 values under univariate regression. All-subsets regression was then applied to the filtered list of covariates to obtain the data generating model for our study. A summary of the data generating model used for each dataset can be found in the supplementary materials.

We used the data generating model and parametric bootstrapping to generate 100 bootstrapped datasets with the same design matrix X and different simulated response vectors Y. Each of the resulting simulated datasets had the same design matrix and error distribution as the real dataset on it was based.

We compared the performance of different parameter and model space prior combinations on these simulated datasets using the following metrics:

We also compared methods based on their out-of-sample predictive performance. We divided each dataset into 100 random 75%–25% train-test splits. We trained the methods on the training data and used the test data to assess the predictive performance using the metrics described below:

We also recorded the average size of the sampled models for each dataset and the average CPU time (in seconds) to carry out BMA for one bootstrapped dataset.

The results are shown in Table 3. We used the combination of the g-prior with g = n as the parameter prior and the Beta-Binomial ( 1 , 1 ) model space prior as the reference. Note that the g-prior with g = n was found to the best parameter prior by [37]. Metrics for all other combinations were calculated relative to the reference metric, and averaged across datasets. Detailed results of performance metrics for the simulation studies based on each of the 14 datasets can be found in the Supplementary materials. The “Score” column contains the average of the scores for PointEst, IntEst, Inference, Prediction and IntPred under each method. We used the Score column to rank the methods.

For each metric, we color the methods based on their performance relative to the reference metric. A method is colored green if it performed similarly or better than the reference method, yellow if it performed somewhat worse, and orange if it performed substantially worse.

For all choices of parameter prior, Beta-Binomial ( 1 , 1 ) was the top scoring model space prior. The three Beta-Binomial versions with g = n were the top three methods across statistical tasks. The Complexity priors with κ = 1 and κ = 2 were the worst performing model space priors. The uniform prior denoted by Ber ( θ = 0.5 ) also performed less well than the Beta-Binomial priors. This ranking of methods was consistent across different performance metrics.

Most parameter and model prior combinations selected sparser models than the g = n and BB ( 1 , 1 ) prior combination, with the exception of methods involving the uniform model space prior. The complexity priors selected very sparse models compared to our baseline, which may be explained by the strong sparsity induced by the prior. This may also explain the poor performance of the complexity priors across statistical tasks. Notably, the rankings of the prior combinations are similar for the different tasks. In particular, the rankings for prediction are consistent with those for point estimation and parameter inference, with a correlation of 0.77 between scores for point estimation and point prediction.

We also note that the EB model space priors tended to outperform the corresponding SDM model space priors when combined with the Hyper-g and EB-local parameter priors. However, the results with the EB model space priors took longer to computer on average because of the optimisation procedure. The hyper-g parameter priors are the slowest due to the integral calculations required in the posterior computation. In general, the Beta-Binomial priors performed better than the Bernoulli and complexity priors.