Revisiting High Dimensional Bayesian Model Selection for Gaussian Regression

TL;DR

This paper analyzes Bayesian model selection in high-dimensional Gaussian regression, focusing on the interplay between prior choices, model complexity, and the asymptotic behavior of Bayes factors, with implications for consistency and false discovery control.

Contribution

It introduces a novel analysis of the Zellner-Siow prior combined with a Poisson prior, revealing conditions for model selection consistency and false discovery control in high dimensions.

Findings

Model selection consistency restricts true model dimension growth.

Large covariate effects are needed to overcome complexity penalties.

Asymptotic behavior of Bayes factors involves stable laws and heavy-tailed variables.

Abstract

Model selection for regression problems with an increasing number of covariates continues to be an important problem both theoretically and in applications. Model selection consistency and mean structure reconstruction depend on the interplay between the Bayes factor learning rate and the penalization on model complexity. In this work, we present results for the Zellner-Siow prior for regression coefficients paired with a Poisson prior for model complexity. We show that model selection consistency restricts the dimension of the true model from increasing too quickly. Further, we show that the additional contribution to the mean structure from new covariates must be large enough to overcome the complexity penalty. The average Bayes factors for different sets of models involves random variables over the choices of columns from the design matrix. We show that a large class these random variables have no moments asymptotically and need to be analyzed using stable laws. We derive the domain of attraction for these random variables and obtain conditions on the design matrix that provide for the control of false discoveries.

Equations

$$\begin{array}[]{rcl}{y}=\alpha\mathbb{1}+{X\beta}+{\epsilon}\text{, where ${\epsilon}\sim N(0,\sigma^{2}{I}_{n\times n})$}.\end{array}$$

$$\begin{array}[]{rcl}T&:=&\left\{i\in\left\{1,\ldots,p\right\}|\beta_{i}\neq 0\right\},\end{array}$$

$$\begin{array}[]{rcl}\beta|\sigma^{2},M_{A}&\sim&\text{N}\left(0,g\sigma^{2}\left(X_{A}^{\prime}X_{A}\right)^{-1}\right),\end{array}$$

$$\begin{array}[]{rcl}{y}&=&\alpha\mathbb{1}+{X\beta}+{\epsilon}\\ \epsilon&\sim&N(0,\sigma^{2}\mathit{I}_{n\times n})\\ \alpha,\sigma^{2}&\sim&\frac{1}{\sigma^{2}}\\ \beta_{A}&\sim&Normal(0,\frac{\sigma^{2}}{\omega}\left(\tilde{X}^{\prime}_{A}\tilde{X}_{A}\right)^{-1})\\ \omega&\sim&\pi(\omega),\end{array}$$

$$\begin{array}[]{rcl}\stackbin[n\to\infty]{}{\text{lim}}\text{Pr}\left(M_{T}|y\right)&=&1,\end{array}$$

$$\begin{array}[]{rcl}\stackbin[n\to\infty]{}{\text{lim}}\stackbin[M_{A}\neq M_{T}]{}{\sum}\frac{m\left(y|M_{A}\right)}{m\left(y|M_{T}\right)}\frac{\pi\left(M_{A}\right)}{\pi\left(M_{T}\right)}&=&0.\end{array}$$

$$\begin{array}[]{rcl}BF_{A:T}=\frac{m\left(y|M_{A}\right)}{m\left(y|M_{T}\right)},\end{array}$$

$$\begin{array}[]{rcl}BF_{A:0}&=&\int_{0}^{\infty}\left(1+\omega\right)^{\frac{n-|A|}{2}}\omega^{\frac{|A|-1}{2}}\left(1-R_{A}^{2}+\omega\right)^{-\frac{n-1}{2}}\pi(\omega)\text{d}\omega,\end{array}$$

$$\begin{array}[]{rcl}BF_{A:T}&=&\frac{\left(1-R_{T}^{2}\right)^{{\frac{n-|T|-1}{2}}}}{\left(1-R_{A}^{2}\right)^{\frac{n-|A|-1}{2}}}\frac{\Gamma\left(\frac{|A|}{2}\right)\Gamma\left(\frac{n-|A|}{2}\right)}{\Gamma\left(\frac{|T|}{2}\right)\Gamma\left(\frac{n-|T|}{2}\right)},\end{array}$$

$$\begin{array}[]{rcl}\pi(M_{A})&=&\pi(|A|){p\choose|A|}^{-1},\end{array}$$

$$\begin{array}[]{rcl}c_{1}p^{-c_{2}}\leq&\frac{\pi(|A|)}{\pi(|A|-1)}&\leq c_{3}p^{-c_{4}},\end{array}$$

$$\begin{array}[]{rcl}\pi(|A|)&=&\frac{\lambda^{|A|}\textrm{exp}\left\{-\lambda\right\}}{|A|!},\end{array}$$

$$||X_i|{|}_2^2=n\text{, for all }i\in \left\{1,\dots ,|F|\right\}.$$

$$\begin{array}[]{c}\stackbin[A\in\left\{1,\ldots,p\right\}]{}{\text{min}}\nu_{\text{min}}\left(\frac{1}{n}X_{A}^{T}X_{A}\right)\geq\zeta_{min}\end{array}$$

$$\begin{array}[]{rcl}\boldsymbol{E}\left[\stackbin[k\notin B]{}{\text{max}}\stackbin[B\subset F;B\supset T]{}{\text{max}}\left(\frac{|<(I-H_{B})x_{k},Z>|}{\sqrt{n}}\right)\right]\leq\frac{\sqrt{\frac{\zeta_{min}}{2}*\textrm{log}\left(\textrm{log}\left(n\right)\right)}}{2},\end{array}$$

$$\text{stackbin}[n\to \infty ]\text{lim}\frac{||X_T{\beta }_T|{|}^2}{n{\sigma }_T^2}=C_1$$

$$\frac{C_2}{|T|}\le {\tau }_T{\beta }_{min}^2\le \frac{C_3}{|T|}$$

$$C_2>\frac{10}{{\zeta }_{min}},$$

$$\text{stackbin}[n\to \infty ]\text{lim}\text{Pr}\left(M_T|y\right)=1$$

$$\begin{array}[]{rcl}\mathcal{M}_{c}&:=&\left\{M_{A};A\supset T,~{}|A|-|T|=c\right\}.\end{array}$$

$$\begin{array}[]{rcl}\stackbin[M_{A}\in\mathcal{M}_{c}]{}{\sum}\frac{m\left(y|M_{A}\right)}{m\left(y|M_{T}\right)}\frac{m\left(M_{A}\right)}{m\left(M_{T}\right)}=\frac{\overline{BF_{A:T}}}{c!}\lambda^{c},\end{array}$$

$$\begin{array}[]{rcl}\text{\lx@cref{creftype~refnum}{eq::bayes_factor}}&\approx&\textrm{exp}\left\{\frac{\eta_{i}}{2}\right\}\frac{\Gamma\left(\frac{|A|}{2}\right)\Gamma\left(\frac{n-|A|}{2}\right)}{\Gamma\left(\frac{|T|}{2}\right)\Gamma\left(\frac{n-|T|}{2}\right)},\end{array}$$

$$\begin{array}[]{rcl}\stackbin[n\to\infty]{}{\text{lim}}\frac{\sum_{i=1}^{m}\delta_{i}-mb_{m}}{a_{m}}&=&S_{\alpha}(1,1,0),\end{array}$$

$$\begin{array}[]{rcl}\stackbin[n\to\infty]{}{\text{lim}}\stackbin[M_{A}\in\mathcal{M}_{c}]{}{\sum}\frac{m\left(y|M_{A}\right)}{m\left(y|M_{T}\right)}\frac{m\left(M_{A}\right)}{m\left(M_{T}\right)}&\neq&0.\end{array}$$

$$\begin{array}[]{rcl}\stackbin[M_{A}\in\mathcal{M}_{c}]{}{\sum}\frac{m\left(y|M_{A}\right)}{m\left(y|M_{T}\right)}\frac{m\left(M_{A}\right)}{m\left(M_{T}\right)}&\approx&\overline{\delta_{i}}*\textrm{log}\left(n\right)^{-\frac{c}{2}}\\ &\to&\frac{\textrm{log}\left(m\right)^{\frac{c}{2}-1}}{\Gamma\left(\frac{c}{2}\right)}\left[S_{\alpha}(1,1,0)+\frac{\textrm{log}\left(m\right)}{\frac{c}{2}}\right]\textrm{log}\left(n\right)^{-\frac{c}{2}}\\ &\to&h(c)\neq 0,\end{array}$$

$$\begin{array}[]{ccc}||H_{A}x_{i}||&\leq&||x_{i}||\end{array}$$

$$\begin{array}[]{ccc}||H_{A}x_{i}||^{2}&=&\langle H_{A}x_{i},H_{A}x_{i}\rangle\\ &=&x_{i}^{\prime}H_{A}x_{i}\\ &=&\langle x_{i},H_{A}x_{i}\rangle\\ &\stackbin[]{(i)}{\leq}&||x_{i}|||H_{A}x_{i}||\end{array}$$

$$\frac{1}{\sqrt{n}}||\left(\mathbf{I}-H_A\right)x_i||\le 1$$

$$\begin{array}[]{c}\nu_{\text{min}}\left(\frac{1}{n}X^{\prime}_{A^{\prime}/A}\left({\mathbf{I}}-H_{A}\right)X_{A^{\prime}/A}\right)\geq\zeta_{min}\end{array}$$

$${\zeta }_{max}\ge \frac{a^{\prime }Aa}{a^{\prime }a}\ge {\zeta }_{min}$$

Full text

1 Introduction

We are considering a generic Gaussian regression problem with a response vector ${y}\in\mathbb{R}^{n}$ and a design matrix $X\in\mathbb{R}^{n\times p}$, such that

[TABLE]

In eq. 1, $\mathbb{1}$ is a $n$-dimensional vector of 1’s, $\alpha$ is an intercept term included in all the models, and $\beta\in\mathbb{R}^{p}$ is a vector of regression coefficients. We assume that the data are generated by the true model $M_{T}$. The dimensions of the true model and the full model are assumed to be growing with the sample size, hence the problem is in the high-dimensional regime. To address the problem, we also make a sparse assumption such that a large portion of the coefficients in the full model are actually zero, and the goal here is to recover the exact support set,

[TABLE]

of the true model.

The high-dimensional sparse solutions for the Gaussian regression has been extensively investigated, and one of the most famous and inspiring methods is the Least Absolute Shrinkage and Selection Operator, so-called the Lasso [17]. It is well-known that the estimation of the Lasso achieves sparsity while suffers from a consistent bias, and there has been tremendous efforts for fixing the inadequacies of the Lasso, such as, the Elastic net [25], the SCAD [5], etc. Meanwhile, Bayesian variants based on the idea of shrinking coefficients have also received great attentions, e.g., the Bayesian Lasso [14], the Horseshoe prior [2] etc. All of these methods attain the sparse solution by shrinking certain coefficients to zero, and are assessed by some measurements of predictions. Rather than focusing on the criteria based on various $\mathit{L}_{q}$ norms $\mathbb{E}||\hat{\beta}-\beta_{T}||_{q}^{q}$, typically $L_{1}$ and $L_{2}$ norms, the current paper pays attentions to variable-selection consistency under 0-1 loss through Bayesian methodology.

The consistency of Bayesian model selection has been established in many papers. The general procedure of Bayesian model selection is usually based on the posterior evidence provided by the Bayes factor, which makes a direct comparison between the marginals of the null and the alternative models [6]. The mathematical nature of the Bayes factor prevents the deployment of the improper prior, since the nuisance constant in the priors would cause the identifiable problem for any statistical inference. To fix this problem, the intrinsic prior was deliberately designed [1] to cancel out the nuisance constants, and the consistency of model selection for pair-wise comparison between the null and alternative model has been established [3], and for the situations where the dimensions of the models are growing with the sample size [12]. However, the pair-wise consistency can not guarantee the posterior probability of the true model going to one due to the massive number of the models.

Another choice of prior distribution for model selection is famous Zellner’s $g$-prior given by,

[TABLE]

where $A\in\left\{1,\ldots,p\right\}$ and $X_{A}$ is the corresponding design matrix composed with the predictors whose indexes are in the set $A$. The popularity of Zellner’s $g$-prior is primarily due to its feasibility of the computation of the corresponding Bayes factor with a closed form. However, this relatively simple form of Bayes factor arises several paradoxes of model selection , and many empirical choices for ‘g’ parameter have been proposed to remedy the situation [10]. Liang et al., [10] proposed the mixtures of $g$-prior by randomizing ‘g’ parameter to overcome the paradoxes, and showed the consistency of model selection with the mixtures of $g$-prior, especially with Zellner&Siow prior [24]. It is also interested to notice that the intrinsic prior can be represented as the mixture of $g$-prior with a proven consistency of model selection [22].

One vital assumption in these papers is that the dimension of the full model is not growing with the sample size, hence the total number of the models is well controlled. Whether or not the size of the true model is growing also has important consequence on the assumptions of the design matrix and the consistency of model selection. With a fixed true model, the procedures, i.e., [13] [21], can tolerate more crucial conditions, such as the size of the full model can be growing at an exponential rate of the sample size or the eigenvalues of the gram matrix can converge to zero. On the other hand, the conditions or assumptions would be harsher, if the dimension of the true model is growing, especially when the situation of the sparsity is close to linear sparsity [20]. Also, when the dimensions are growing, the learning rate of the Bayes factor may break down, e.g., the models nesting the true model become indistinguishable to the true model. Hence, the prior on the model space is often required to offer extra penalties of the dimensionality. These priors, e.g., the sparsity prior [4] [23] or the truncated Poisson prior [21], are all designed to make penalization on the overfitted models and control the model size for models that are not identifiable from the data [9].

In this paper, we show that the model selection consistency through the truncated Poisson prior on the model space and the modified Zellner&Siow prior on the regression parameters. The paper is organized as follows. In Section 2, the main results is presented. We review the assumptions of the design matrix used in other papers, list the assumptions that we adopt, and emphasize the corresponding impacts on the consistency. The difference between the sparsity prior and the truncated Poisson prior is highlighted. In Section 3, we show an interesting utility of the stable law and the domain of attraction on the overfitted models. Section 4 contains the conclusion and possible directions of future works. The proofs of the theorems are in Section 5.

1.1 Notification

As stated above, let $A$, $T$, and $F$ denote the indices set of the predictors in the testing model, the true model, and the full model respectively. For a vector $\beta\in\mathbb{R}^{p}$ and a set $A\in\left\{1,\ldots,p\right\}$, $\beta_{A}$ is the vector of $\beta_{i}$s for $i\in A$, and $|A|$ is the cardinality of $A$. Let $M_{0}$ denote the null model, the model only contains the intercept term.

2 Main results

2.1 Mixture of $g$-prior & the Bayes factor

Now we specify the prior distributions of the parameters. For any model $M_{A}$, the model and the priors are

[TABLE]

where $\omega=\frac{1}{g}$, $\tilde{X}^{\prime}_{A}\tilde{X}_{A}=\frac{1}{n}X_{A}^{\prime}X_{A}$, which corresponds to the recommendation of unit information prior [6]. Zellner&Siow [24] place a multivariate Cauchy distribution on the coefficient vector, which can be represented as a mixture of normal distributions with the mixing parameter $\omega\sim\text{Gamma}\left(\frac{1}{2},\frac{1}{2}\right)$. The Zellner-Siow prior is a multivariate extension of Jeffreys’s idea on the normal mean hypothesis problem, where Jeffrey argued that Cauchy prior is the simplest form to fix the paradoxes mentioned above.

The goal of this paper is to show the selection consistency, such that

[TABLE]

which is equivalent to show that

[TABLE]

The Bayes factor is defined as the ratio of the marginal densities between two competing models in eq. 6, such as

[TABLE]

where $m\left(y|M_{A}\right)$ is the marginal density of the model $M_{A}$ after integrating out the regression parameters. For simple computation, it is worth noting that the marginal density $m\left(y|M_{A}\right)$ is invariant if switching $\tilde{X}^{\prime}_{A}\tilde{X}_{A}$ by $X_{A}^{\prime}X_{A}$ and assigning $\text{Gamma}\left(\frac{1}{2},\frac{n}{2}\right)$ to $\omega$. Integrating out the parameters $\alpha$, $\beta_{A}$, and $\sigma^{2}$ is straight forward, and the Bayes factor can be represented as an univariate integration of $\omega$, such as

[TABLE]

where $R_{A}^{2}=\frac{y^{\prime}\left(H_{A}-H_{1}\right)y}{y^{\prime}(I-H_{1})y}$ is the coefficient of determination and $H_{A}$ is the projection matrix generated by the design matrix $X_{A}$. In [10], the authors suggested to approximate the integral in eq. 8 by the Laplace method [18], which requires to solve a cubic function of $\omega$. Since the dimension of the true model is assumed to be fixed in [10], the solution is asymptotically stable as $n\to\infty$, which is not the case in the current set-up. Actually, the unsatisfied approximation is the main reason that the Z&S prior never got popular in the first place, even with appealing statistical properties.

In this paper, we choose the Beta-prime distribution with the shape parameters $\frac{1}{2}$ and $\frac{n-|A|-1}{2}$ as a modified version of the Z&S prior. The parameters are deliberately chosen to cancel out the term $\left(1+\omega\right)^{\frac{n-|A|}{2}}$ in eq. 8, also to maintain an asymptotically equivalent behavior of the density function to the Zellner&Siow prior at the tail area and the origin. The same set-up has been adopted in [11], and the resulting Bayes factor between the testing model $M_{A}$ and the true model $M_{T}$ can be shown as

[TABLE]

where $\Gamma(\bullet)$ is the Gamma function.

Remark 1

The major difference between the current paper and [23] is the choice of the original $g$-prior or the mixture of $g$-prior. Certainly, the original $g$-prior with an empirical choice on $g$, e.g., $g=\text{max}\left\{n,p^{2}\right\}$ or $g=p^{2\alpha}$ [23] [7], does exhibit certain flexibility for the model selection. However, the information paradox associated with the original $g$-prior demands more assumptions on the signals to prove the consisitency [10] [16] [23]. Also, when the null model is the true model, consistency only holds true for the Zellner-Siow prior, but does not hold for the empirical $g$-prior [10]. The advantages of adopting the mixture of $g$-prior are easy to observe. First, it avoids the tuning step of $g$, since it has been integrate out for the marginals. Second, it can be shown that either when $M_{T}=M_{N}$ or $M_{T}\neq M_{N}$, the Bayes factor associated with the mixture of $g$ holds the consistency of the selection, where it is not the case for the empirical $g$-prior. ∎

2.2 Model prior

In the fixed dimension regime, the prior of models is usually given equal prior probability and would be canceled out during the procedure. Also, the convergence of the pair-wise comparison is enough for establishing the consistency, since the summation in eq. 6 is only over a finite number of models. However, things are different in the high-dimensional regime, and the convergence of the pair-wise comparison doesn’t ensure the posterior selection consistency in eq. 5.

For a generic model $M_{A}$, the prior probability of choosing $M_{A}$ can be represented as

[TABLE]

where $\pi(|A|)$ is the prior probability of the size $|A|$, and all the models with same size equally share the same probability mass $\pi(|A|)$. In [4] [13] [23], the authors chose the sparse prior to introduce strong penalty on the dimension with the form

[TABLE]

for positive $c_{i}$. If we assumed $p=O(n)$, then eq. 11 is a very strong penalization, and the consistencies have been proved in the papers adopted the sparse prior. Equation 11 is also equivalent to $\pi(\gamma_{j}=1)=O(p^{-c_{i}})$, which is the prior inclusion probability of adding one generic covariate. From this perspective, it is obvious that the prior inclusion probability only depends on the size of the full model without any consideration of the size of the current testing model, which is less adjustable or flexible for a variety of situations. Also, eq. 11 lacks a properly realistic explanation to the practitioners at face value. Last but not the least, when comparing the growing true model to a finite model nested in the true model, $\frac{\pi\left(M_{A}\right)}{\pi\left(M_{T}\right)}$ in eq. 6 increases in an exponential rate of the sample size, which imposes harsh and unrealistic conditions on design matrix and signal noise ratio to achieve model selection consistency.

We propose to use a truncated Poisson prior on the size of models, such that

[TABLE]

where $\lambda$ is the rate parameter, which can be tuned by practitioners. The derivation of the truncated Poisson prior comes from using the self-similarity property for model spaces with finite $p$ and letting $p\to\infty$ with an easy proof in [21]. The prior inclusion probability of the truncated Poisson distribution is only associated with the testing model size $|A|$, which is more adjustable than the sparsity prior. When dealing with the overfitted models, the truncated Poisson prior provides a well-designed penalty to convert the summation of the Bayes factors generated by the models with same size to its arithmetic average, hence the convergence is readdressed by the probabilistic property of the average of the Bayes factor instead of the overly-strong dimensional penalty such as the sparsity prior. Also, the average of the Bayes factor is an average of the random variables generated by randomly choosing the extraneous predictors, which leads to an interesting application of the stable law and the domain of attraction. We will give a more detailed explanation in Section 3.

2.3 Assumptions

The assumptions that we will make in this section are crucial to the proof of the consistency. These assumptions reflect the integrity of the problem, i.e., the model with too strict assumptions would lost practical usage to real-world data yet the model without any assumptions could not be proved to attain the consistency.

Assumption 1** (Conditions on the design matrix).**

• (i)

Assume all the columns of the design matrix are standardized, such as

[TABLE]

• (ii)

There exists a positive constant $\zeta_{min}$, such that

[TABLE]

where $\nu(H)$ is the eigenvalue of a symmetric matrix $H$.

• (iii)

We assume a standardized errors constraining assumption, such as,

[TABLE]

where $Z\sim N(0,\mathit{I}_{n\times n})$, and $H_{B}$ is the projection matrix generated by the design matrix $X_{B}$.

∎

Assumption 2** (Condition on the mean structure).**

Given $X_{T}$, $\beta_{T}$, and $\sigma_{T}^{2}$, assume

[TABLE]

where $C_{1}$ is a finite positive constant. ∎

Assumption 3** (Conditions on the growing rates).**

The growing rate of the full model can be only as fast as a fraction of the sample size, i.e., $\stackbin[n\to\infty]{}{\text{lim}}\frac{|F|}{n}=f<1$. ∎

Assumption 4** (Condition on the minimum signal).**

Let $\beta_{min}$ denote $\stackbin[i\in T]{}{\text{min}}|\beta_{i}|$, then we assume that $\tau_{T}||\beta_{T}||_{2}^{2}=O(1)$, which implies $\tau_{T}\beta_{min}^{2}=\Theta(\frac{1}{|T|})$. Specifically, the minimum of the true coefficients is bounded above and below as

[TABLE]

where $C_{2}$ is a finite constant bounded above by zero, and $C_{3}$ is a finite constant implied by 2. ∎

Remark 2

In 1, $(i)$ and $(ii)$ are typically mild conditions that appear in many literatures considering the consistency, e.g., [13] [21] [23] [16]. $(ii)$ in 1 guarantees that each potentially testing model is not too close to be distinguish to each other. $(iii)$ is to constrain the projection of the standardized errors on the space generated by any extraneous predictors excluding the space generated by the design matrix nesting the true model. This assumption is also a mild condition just to prevent any extreme behavior from the overfitted model and their extraneous predictors.

2 is a reasonable condition putting onto the true mean structure, which surprisedly was not considered by many literatures explicitly. Without 2, under some circumstances, the true mean would grow to infinity with the sample size, which doesn’t make any sense. 3 assumes the growing rate of the full model is linear to the sample size. 4 is related to the information-theoretical capacity of the model to recovery the exact support [20]. As the growing rate of the true model hitting the limit, i.e., $\frac{n}{\textrm{log}\left(n\right)}$, the minimum of the signal has to correspond to such harsh condition to separate itself from the noise. The reason for this rate is shown in proposition 1. ∎

2.4 Result

Proposition 1**.**

Suppose that 1, 2, 3, and 4 hold, the fastest growing rate of the true model that the procedure can tolerate is $\mathcal{O}\left(\frac{n}{\textrm{log}\left(n\right)}\right)$.

Proof.

See Section 5. ∎

Here are the main theorem of this paper.

Theorem 1**.**

Suppose that 1, 2, 3, and 4 hold. Suppose $|T|=\frac{n}{\textrm{log}\left(n\right)}$. If

[TABLE]

then, we have

[TABLE]

with probability at least $1-\left[\textrm{log}\left(n\right)\right]^{-\frac{1}{4}}$.

Proof.

See Section 5. ∎

The constant $C_{2}$ characterizes the relationship between the minimums of the eigenvalues and the signals, which all hinges on the growing rate of the true model. The probability in theorem 1 is determined by the worst scenario, where the testing models are the overfitted models with only one extraneous predictor. The reason behind such a slowly convergent rate is due the fast growing rate of the true model.

It is reasonable to consider the scenario when the growing rate of the true model is $\mathcal{O}(n^{d})$ for $d<1$. It can be shown that the relatively strict assumptions can be loosen, and the probability associated with the consistency would be faster on a different order.

3 Stable law

One interesting encounter from the proof is an application of the stable law and the domain of attraction, when dealing with the overfitted models. Let $\mathcal{M}_{c}$ denote the model set of overfitted models with number of $c$ extraneous predictors, such that

[TABLE]

Then, it is straight forward to show

[TABLE]

where $\overline{BF_{A:T}}$ is the arithmetic average of the Bayes factor over the model set $\mathcal{M}_{c}$. As $|F|-|T|\to\infty$, each Bayes factor can be seen as a random variable which is randomly chosen from the set $\mathcal{M}_{c}$. Furthermore, as shown in eq. 9, the Bayes factor can be well approximated as

[TABLE]

where $\eta_{i}\sim\chi^{2}(c)$ for $i=1,\ldots,{|F|-|T|\choose c}$, and the covergent rate of the fractions associated with the terms of the Gamma functions is $\mathcal{O}\left(\textrm{log}\left(n\right)^{-\frac{c}{2}}\right)$. Hence, the average over the Bayes factor is equivalent to the average over the random variable $\textrm{exp}\left\{\frac{\eta_{i}}{2}\right\}$, which is generated by models $M_{A}\in\mathcal{M}_{c}$. However, a simple integration can show that the r.v. $\delta_{i}=\textrm{exp}\left\{\frac{\eta_{i}}{2}\right\}$ doesn’t have any finite moment. Hence, to show the convergence of the average, the stable law is required. The next theorem concludes the convergence of the average for i.i.d. situation.

Theorem 2**.**

Suppose $\delta_{i}=\textrm{exp}\left\{\frac{\eta_{i}}{2}\right\}$ and $\eta_{i}\stackbin[i.i.d.]{}{\sim}\chi^{2}(c)$, where $c$ is a finite positive integer and $i=1,\ldots,m$. Then given the constants $a_{m}=m*\frac{\textrm{log}\left(m\right)^{\frac{c}{2}-1}}{\Gamma\left(\frac{c}{2}\right)}$ and $b_{m}=\frac{\textrm{log}\left(a_{m}\right)^{\frac{c}{2}}}{\Gamma\left(\frac{c}{2}+1\right)}$, we have

[TABLE]

where $S_{\alpha}(1,1,0)$ is a Cauchy random variable.

This theorem reveals the fact that the average of the Bayes factor converges to a Cauchy r.v. and a slowing varying function associated with the total number of the models in the set $\mathcal{M}_{c}$. Since the convergent power provided by the Bayes factor is only $\mathcal{O}\left(\textrm{log}\left(n\right)^{-\frac{c}{2}}\right)$, the corollary next shows that it is impossible to show the convergence for the overfitted models without assumptions controlling the behavior of the errors.

Corollary 1**.**

Given a dimensional difference $c$ and the corresponding model set $\mathcal{M}_{c}$ defined in eq. 20, we have

[TABLE]

Proof.

Let $m={|F|-|T|\choose c}\asymp n^{c}$ denote the number of the models in $\mathcal{M}_{c}$. We have seen that the summation above can be represented the average of the r.v. $\delta_{i}$, and the Gamma functions in eq. 22 only has the rate of $\mathcal{O}\left(\textrm{log}\left(n\right)^{-\frac{c}{2}}\right)$. Suppose $\delta_{i}$s are i.i.d., hence

[TABLE]

where $h(c)$ is a function of the difference $c$, and not equal to 0 if $c$ is a finite positive integer. ∎

This corollary shows that the consistency requires extra assumptions, e.g., (iii) in 1, to control the projection of the errors on the subspace generated by the non-true predictors, when the growing rate of the true model has reached the limit, i.e., $\frac{n}{\textrm{log}\left(n\right)}$. If the growing rate of the true model is slower than the limit, then the convergent rate provided by the Bayes factor will overcome the slowly varying function in $a_{m}$ and $b_{m}$.

4 Conclusion

In this paper, the model selection problem of Gaussian regression in high-dimensional regime has been studied. The mixture of $g$-prior and the truncated Poisson prior have been proposed to show the consistency of model selection. We presented the consistency theorem under the situation of extreme growing rate. It showed that when the limit of the growing rate of the true model is approached, the consistency theorem requires more strict assumptions. The stable law has been applied to make the argument of the unidentifiable situation for the overfitted models, when there is no extra assumptions on the extraneous predictors and the growing rate of the true model reaches the limit.

For the next step, it is interesting to further investigate the stable law under the dependent situation due to the possible dependency structure between the extraneous covariates. Another direction would be to discuss the situation when the size of the full model is greater than the sample size. The authors in [23] simply truncated the model space by assigning zero probability to the models whose size is greater than the sample size. We consider to adopt the PCA method to reduce the dimension as the first step of the model selection.

5 Appendix

5.1 Lemmas

Let $H_{A}$ denote the projection matrix onto the span of $\left\{X_{A},A\in F\right\}$.

Lemma 1**.**

A projection is a non-expensive mapping[23]. In the words, for any column $x_{i}$, $i\in F$,

[TABLE]

Proof.

[TABLE]

(i) is by the Cauchy-Schwarz inequality. ∎

Lemma 2**.**

For any column $x_{i}$, $i\in F$,

[TABLE]

Proof.

Direct result from use of Lemma 1 and $(i)$ in 1. ∎

Lemma 3**.**

Under $(ii)$ 1, for any pair $A\in F$, $A^{\prime}\in F$, $A\subset A^{\prime}$, we have

[TABLE]

Proof.

W.t.l.g., assume $X_{A^{\prime}}=\left\{x_{A},x_{A^{\prime}/A}\right\}$, then by the formula of the blockwise inversion, one can show that the lower right corner of the matrix $\left(\frac{1}{n}{X_{A^{\prime}}}^{\prime}X_{A^{\prime}}\right)^{-1}$ is $\left(\frac{1}{n}X^{\prime}_{A^{\prime}/A}\left({\mathbf{I}}-H_{A}\right)X_{A^{\prime}/A}\right)^{-1}$. The rest follows. ∎

Lemma 4**.**

For any $n\times 1$ vector $a$ and a generic symmetric matrix $A$, we have

[TABLE]

where $\zeta_{min}$ and $\zeta_{max}$ are the smallest and largest eigenvalues of $A$.

Proof.

See [15]. ∎

Lemma 5**.**

Suppose $M_{A}$ is in the model set $\mathcal{M}_{c,k}:=\left\{M_{A}:|T|-|A|=c,A/T=k\right\}$, then it can be shown that

[TABLE]

Proof.

Notice that $T/A=c+k$.

[TABLE]

∎

Remark 3

there is a question. ∎

5.2 Proof of Theorem 1

5.2.1 $|A|<|T|$

We begin by proving the consistency under the situation $\mathcal{M}_{c,k,s}:=\left\{M_{A}:|T|-|A|=c,A/T=k\right\}$. First, notice that there are ${|T|\choose c+k}{|F|-|T|\choose k}$ models associated with $\mathcal{M}_{c,k}$ for specific $c$ and $k$. The Bayes factor for any $M_{A}\in\mathcal{M}_{c,k}$ is

[TABLE]

The terms associated with the coefficients of determination can be shown as

[TABLE]

where $(i)$ is due to Lemma 5. It can also be shown that

[TABLE]

Let $\stackbin[n\to\infty]{}{\text{lim}}\frac{n+nC_{1}}{n-|A|+(c+k)\beta_{min}^{2}\tau_{T}\zeta_{min}n}=C_{4}=O(1)$. Then,

[TABLE]

The summand then can be represented as

[TABLE]

If $|T|=\frac{n}{\textrm{log}\left(n\right)}$, then the consistency requires that $C_{2}>\frac{10}{\zeta_{min}}$. On the other hand, suppose that $|T|=n^{d}$ for $0\leq d<1$. Then $\beta_{min}^{2}\tau_{T}\zeta_{min}n=O(n^{1-d})$, which is growing faster than $\textrm{log}\left(n\right)$. This implies that the $\zeta_{min}$ can converge to zero with a mild speed under such circumstance.

5.2.2 $|A|>|T|$ and $A\cap T\neq T$

Under this circumstance, the dimensional penalty helps the procedure to pick up the true model comparing to the previous situation. Define the model set as

$\mathcal{M}_{c,k,l}:=\left\{M_{A}:|A|-|T|=c,T/A=k\right\}$, then the total number of models given any $c$ and $k$ is ${|T|\choose k}*{|F|-|T|\choose k+c}$. The terms in the summand excluding the terms associated with the coefficient of determination can be shown as

[TABLE]

It is straight forward to show that

[TABLE]

where $\stackbin[n\to\infty]{}{\text{lim}}\frac{n-|A|+n\zeta_{min}||\beta_{T/A}||_{2}^{2}\tau_{T}}{n+C_{1}n}=C_{5}\leq 1$. Combining Equation 38 and Equation 39, it gives us

[TABLE]

Given the condition $C_{2}>\frac{10}{\zeta_{m}in}$ in the LABEL:tm::under-fitted, the consistency holds for this situation.

5.3 Overfitted Model

Suppose that $|T|=\frac{tn}{\textrm{log}\left(n\right)}$, for $0<t<1$. Define the overfitted model set

$\mathcal{M}_{c}:=\left\{M_{A}:A\supset T,|A|-|T|=c\right\}$ associated with a given finite $c$. Notice that the summation over the model set $\mathcal{M}_{c}$ is

[TABLE]

The individual Bayes factor is asymptotically equivalent to

[TABLE]

where $\xi_{A}=\tau_{T}Y^{\prime}(H_{A}-H_{T})Y$ and $t=1$ for now. To show $\overline{BF}\to 0$, it suffices to show that

[TABLE]

This can be achieved by showing the probability of the event defined as Equation 43 going to 1, such as

[TABLE]

where $\epsilon_{n}=\left[\textrm{log}\left(n\right)\right]^{-\frac{c}{4}}\stackbin[n\to\infty]{}{\to}0$, and $w^{\prime}\stackbin[n\to\infty]{}{\to}\infty$ to guarantee the convergence. Notice that

[TABLE]

which means that it suffices to show

[TABLE]

Define a function

[TABLE]

For any $Z$, $Z_{1}\in\mathcal{R}^{n}\sim N(0,I_{n\times n})$, it can be shown

[TABLE]

where $(i)$ is due to the normalization assumption. Equation 48 shows that $V(Z)$ is a Lipschitz function, and $||V(\bullet)||_{Lip}=1$. By the Theorem 5.2.2 in [19], and also can be found in [8], we have

[TABLE]

Now we bound the expectation of $V(Z)$ by introducing the following assumptions

Assumption 5**.**

[TABLE]

∎

This assumption ensures that

[TABLE]

This completes the convergence as

[TABLE]

where $\epsilon_{n}=\left[\textrm{log}\left(n\right)\right]^{-\frac{c}{4}}$.