\section{Introduction}
\label{sec:intro}
Suppose the linear regression model is used to relate $Y$ to the
$p$ potential predictors $x_1, \dots, x_p$,
\begin{equation} \label{full-model}
\bm{y} = \alpha \bm{1}_n+\bm{X}_F\bm{\beta}_F + \sigma_F\bm{\epsilon}_F ,
\end{equation}
 where the subscript $F$ refers to the full model $\mathcal{M}_F$.
In the model \eqref{full-model}, $\alpha$ is an unknown intercept parameter,
$\bm{1}_n$ is an $n\times 1$ vector of ones,
$\bm{X}_F=(\bm{x}_1,\dots, \bm{x}_p)$ is an $n \times p$ design matrix,
and $\bm{\beta}_F$ is a $p \times 1$ vector of unknown
regression coefficients.
In the error term of \eqref{full-model}, $\sigma_F$ is an unknown scalar
and $\bm{\epsilon}_F $ has a spherically symmetric (SS) distribution with the density
$f(\|\bm{\epsilon}_F\|^2)$, $E[\bm{\epsilon}_F]=\bm{0}_n$ and $\mathrm{Var}[\bm{\epsilon}_F]=\bm{I}_n$.
We assume that the columns of $\bm{X}_F$ have been standardized so that for 
$1 \leq i \leq p$, $\bm{x}'_i\bm{1}_n = 0$ and $\bm{x}'_i\bm{x}_i/n = 1$. 

We shall be particularly interested in the variable selection 
problem where we would like to select an unknown
subset of the effective predictors.
It will be convenient throughout to index each of these $2^p$ possible
subset choices by the vector
\begin{equation*}
\bm{\gamma} =(\gamma_1, \dots,\gamma_p)'
\end{equation*}
where $\gamma_i=0$ or $1$. 
We use $q_{\bm{\gamma}}=\bm{\gamma}'\bm{1}_p$ to denote the size of
the $\bm{\gamma}$th subset.
The problem then becomes that of selecting a submodel
of \eqref{full-model}
\begin{equation} \label{submodel-gamma}
\bm{y} = \alpha \bm{1}_n+\bm{X}_{\bm{\gamma}} \bm{\beta}_{\bm{\gamma}} + \sigma_{\bm{\gamma}}\bm{\epsilon}_{\bm{\gamma}}.
\end{equation}
In \eqref{submodel-gamma}, $\bm{X}_{\bm{\gamma}}$ is the $n \times q_{\bm{\gamma}}$ 
matrix whose columns 
correspond to the $\bm{\gamma}$th subset of $x_1,\dots,x_p$, $\bm{\beta}_{\bm{\gamma}}$
is a $q_{\bm{\gamma}} \times 1$ vector of unknown regression coefficients.
Let $\mathcal{M}_{\bm{\gamma}}$ denote the submodel given by \eqref{submodel-gamma}. 
We assume the error term $\bm{\epsilon}_{\bm{\gamma}} $ has 
the SS 
density
\begin{equation}\label{bep_sim_f}
 \bm{\epsilon}_{\bm{\gamma}}\sim f(\|\bm{\epsilon}_\bm{\gamma}\|^2)
\end{equation}
for all $\bm{\gamma}$, 
with $E[\bm{\epsilon}_{\bm{\gamma}}]=\bm{0}_n$ and $\mathrm{Var}[\bm{\epsilon}_{\bm{\gamma}}]=\bm{I}_n$.
Further $\sigma_{\bm{\gamma}}$ is an unknown scalar in the error term.
We note that, in most earlier studies, the error terms 
in linear models have been assumed to have a Gaussian distribution, 
e.g.~as in \cite{George-Foster-2000} and \cite{Liang-etal-2008}. 

In this paper, we assume that $n > p+1$ (the so called classical setup)
and $\{\bm{x}_1,\dots,\bm{x}_p\}$ are linearly independent, 
which implies that 
\begin{equation}
\mathrm{rank}\ \bm{X}_F=p, 
\quad \mathrm{rank}\ \bm{X}_{\bm{\gamma}}=q_{\bm{\gamma}}.
\end{equation}
We also assume in much of the paper 
that the null model $\mathcal{M}_N$ ($q_{\bm{\gamma}}=0$ or $ \bm{\gamma}=(0,\dots,0)'$)
is not a possible model, that is, the number of possible models is $2^p-1$, 
rather than, $2^p$, although
in Section \ref{sec:marginal}, we indicate an alternative development
which allows all $2^p$ possible models.

A Bayesian approach to this problem entails the specification of 
prior distributions on the models $\pi_{\bm{\gamma}} =\mathrm{Pr}(\mathcal{M}_{\bm{\gamma}})$, 
and on the parameters  
$\alpha,\bm{\beta}_{\bm{\gamma}}, \sigma_{\bm{\gamma}}$ of each model.  
For each such specification, 
of key interest is the posterior probability of $\mathcal{M}_{\bm{\gamma}}$
given $\bm{y}$,
\begin{equation} \label{posterior-2}
 \mathrm{Pr}(\mathcal{M}_{\bm{\gamma}}|\bm{y})
=\frac
{\pi_{\bm{\gamma}} m_{\bm{\gamma}}(\bm{y})}
{\sum_{\bm{\gamma}} \pi_{\bm{\gamma}} m_{\bm{\gamma}}(\bm{y})}
=\frac
{\pi_{\bm{\gamma}} \mathrm{BF}_{\gamma:F}}
{\sum_{\bm{\gamma}} \pi_{\bm{\gamma}} \mathrm{BF}_{\gamma:F}},
\end{equation}
where $\pi_N=\mathrm{Pr}(\mathcal{M}_N)=0$ is assumed
as mentioned in the above (although see Remark \ref{rem:null}).
In \eqref{posterior-2}, $m_{\bm{\gamma}}(\bm{y})$ is the marginal density under
$\mathcal{M}_{\bm{\gamma}}$ and 
$\mathrm{BF}_{\gamma:F}$ is the Bayes factor for comparing
each of $\mathcal{M}_{\bm{\gamma}}$ to the full model $\mathcal{M}_F$
which is defined as
\begin{equation*}
 \mathrm{BF}_{\gamma:F}
= \frac{m_{\bm{\gamma}}(\bm{y})}{m_F(\bm{y})},
\end{equation*}
where $m_F(\bm{y})$ is the marginal density under the full model.
In Bayesian model selection,
\begin{equation}\label{eq:relation_m_BF}
\mathop{\mathrm{argmax}}_{\bm{\gamma}} \mathrm{Pr}(\mathcal{M}_{\bm{\gamma}}|\bm{y})=
\mathop{\mathrm{argmax}}_{\bm{\gamma}} \pi_{\bm{\gamma}} m_{\bm{\gamma}}(\bm{y})=
\mathop{\mathrm{argmax}}_{\bm{\gamma}} \pi_{\bm{\gamma}} \mathrm{BF}_{\gamma:F},
\end{equation}
is typically selected as the best model.

In this paper, the main focus is on $ \mathrm{BF}_{\gamma:F} $, not $\pi_{\bm{\gamma}} $.
Hence our main aim is to propose and study specifications for
the prior distribution of the parameters for
each submodel $\mathcal{M}_{\bm{\gamma}}$.
In particular, the joint density we consider for $\mathcal{M}_\bm{\gamma} $ has the form
\begin{equation} \label{intro-prior}
p_{\bm{\gamma}}(\alpha,\bm{\theta}_{\bm{\gamma}}, \sigma_{\bm{\gamma}}|\nu)
\propto 
\sigma_{\bm{\gamma}}^{-\nu-1} \|\bm{\theta}_{\bm{\gamma}} \|^{-q_{\bm{\gamma}}+\nu}
\end{equation}
for 
$\bm{\theta}_{\bm{\gamma}}=(\theta_1,\dots,\theta_{q_{\bm{\gamma}}})'=(\bm{X}'_{\bm{\gamma}} \bm{X}_{\bm{\gamma}})^{1/2}\bm{\beta}_{\bm{\gamma}}$ 
and a non-random hyper-parameter $\nu$ with $0< \nu< q_{\bm{\gamma}}$. 
Since the term including $\bm{\theta}_{\bm{\gamma}}$ in the prior above, 
$ \|\bm{\theta}_{\bm{\gamma}} \|^{-q_{\bm{\gamma}}+\nu}$ for $0<\nu<\min(2,q_{\bm{\gamma}})$,
is known as a sub-harmonic function, that is,
\begin{equation*}
\textstyle{\sum_{i=1}^{q_{\bm{\gamma}}} \left\{\partial^2/\partial \theta_i^2\right\}
\|\bm{\theta}_{\bm{\gamma}} \|^{-q_{\bm{\gamma}}+\nu}>0},
\end{equation*}
we call the prior given by \eqref{intro-prior} a sub-harmonic prior.
We will show that such priors lead to robust Bayes factors, in the sense that
each Bayes factor does not depend on the form of the underlying error 
distribution. We also show that the resulting procedure has model selection consistency
and is also coherent.

The organization of this paper is as follows.
In Section \ref{sec:prior}, we give details of the prior distribution.
In Section \ref{sec:marginal},
we show that the Bayes factor with respect to the above prior 
is given by
\begin{equation}\label{BF:intro}
\mathrm{BF}_{\gamma:F}(\nu)=\mathrm{BF}^G_{\gamma:F}(\nu)
\end{equation}
where
\begin{equation} \label{BFG:intro}
\mathrm{BF}^G_{\gamma:F}(\nu)
=\frac
{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-q_{\bm{\gamma}}-1}{2}}
\{g(1-R_{\bm{\gamma}}^2)+1 \}^{-\frac{n-1}{2}}\, dg}
{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-p-1}{2}}
\{g(1-R_F^2)+1 \}^{-\frac{n-1}{2}}\, dg},
\end{equation}
for $0<\nu<\min_{\bm{\gamma}} q_{\bm{\gamma}}=1$.
In \eqref{BFG:intro}, $\mathrm{BF}^G_{\gamma:F}(\nu) $ 
is the Bayes factor for standard Gaussian errors and
$R_{\bm{\gamma}}^2$ and $R_F^2$ are the coefficient of determination 
under the submodel $\mathcal{M}_{\bm{\gamma}}$ and the full model $\mathcal{M}_F$, respectively.
From \eqref{BF:intro}, 
the Bayes factor does not depend on the same SS  sampling density.
Hence, even when there is no specific information about 
the form of the error distribution of each model 
(other than spherical symmetry), 
it is not necessary to specify the exact form of the sampling density. 
It suffices to assume it is Gaussian.
As far as we know, in the area of Bayesian variable selection with 
shrinkage priors or Zellner's $g$-priors,
the sampling density has been assumed to be Gaussian and
this kind of robustness result has not yet been studied.
Originally analogous robustness results were derived by \cite{Maruyama-2003b},
\cite{Maru-Straw-2005} and others in the problem of estimating regression
coefficients with the Stein effect
and it was this type of result that lead us to search for a similar phenomenon in 
the context of model selection. Note that we use the term ``robustness''
in this sense of distributional robustness over the class of SS  error distributions.
We specifically are not using the term to indicate a high breakdown point.
The use of the term ``robustness'' in our sense is however common
(if somewhat misleading) in the context of insensitivity to the
error distribution in the context of shrinkage literature.
In Section \ref{sec:Laplace}, by use of the Laplace approximation,
we approximate the Bayes factor given by \eqref{BF:intro} as 
$ \mathrm{BF}^G_{\gamma:F}(\nu)\approx\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$
\begin{equation}\label{eq:intro.approxBF}
\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu) =
\frac{\varphi(q_{\bm{\gamma}}-\nu,1-R^2_{\bm{\gamma}})}{\varphi(p-\nu,1-R^2_F)}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]
\end{equation}
as $n\to\infty$ where 
\begin{equation*}
\varphi(s,r)= rs^{s-1}\{(1/r-1)e\}^{-s}
\end{equation*}
and $\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]$ is the BIC based alternative 
for standard Gaussian errors
\begin{equation*}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]= \left\{\frac{(1-R_{\bm{\gamma}}^2)^{-n}n^{-q_{\bm{\gamma}}}}{(1-R_F^2)^{-n}n^{-p}}\right\}^{1/2}.
\end{equation*}
(See, e.g.~\cite{Hastie-etal-2009}, Chapter 7.) 
Since $\varphi(s,r)$ does not depend on $n$, 
\eqref{eq:intro.approxBF} is asymptotically equivalent to BIC with
a simple $O(1)$ rational correction function
depending upon $\nu$ as well as the $R$-squares and the numbers of predictors.
Actually this is a special case of Theorem \ref{thm:Laplace} in which
several Bayes factors under Gaussian errors which have been 
proposed in earlier studies, are shown to have similar asymptotic
approximations.
While the main theme in this paper is to develop the relationship 
\eqref{BF:intro} under sub-harmonic priors, we believe 
that this asymptotic equivalence is another noteworthy contribution, in particular, from 
the computational point of view.
In Section \ref{sec:consistency}, we show that our Bayes factor 
has model selection consistency 
uniformly over the class of SS  error distributions,
as $n \to \infty$ and $p$ is fixed.
It also follows from these results that several model selection 
methods recently studied in the literature for Gaussian errors
have model selection consistency for the entire class of SS  error distributions.
We provide illustrations of the method and comparisons with other methods
using both simulated and real data in Section \ref{sec:ex}.
We give concluding remarks in Section \ref{sec:cr}. 
The Appendix presents some of the more technical proofs.

\section{Prior distributions}
\label{sec:prior}
In this section, for each submodel, we give a prior joint density of a form
\begin{equation*}
p(\alpha,\bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}})= p(\alpha)p(\eta_{\bm{\gamma}})
p(\bm{\beta}_{\bm{\gamma}}|\eta_{\bm{\gamma}}) ,
\end{equation*}
where $\eta_{\bm{\gamma}}=1/\sigma_{\bm{\gamma}}^2$.
The natural choice of priors for location ($\alpha$) and scale ($\eta_{\bm{\gamma}}$) are
\begin{equation} \label{improper-prior-alpha}
p^I_\alpha(\alpha)=I_{(-\infty,\infty)}(\alpha),
\end{equation}
and 
\begin{equation} \label{improper-prior-sigma^2}
p^I_{\eta}(\eta_{\bm{\gamma}})=\eta_{\bm{\gamma}}^{-1}I_{(0,\infty)}(\eta_{\bm{\gamma}}).
\end{equation}
In \eqref{improper-prior-alpha} and \eqref{improper-prior-sigma^2},
the superscript $^I$ means that the prior density is improper.
Since \eqref{improper-prior-alpha} and \eqref{improper-prior-sigma^2} 
have invariance to location and scale transformation,
respectively, they are considered by many as non-informative objective
priors. 
The problem is that they are improper 
and hence determined only up to an arbitrary multiplicative constant.
In this paper, the use of improper priors is justified through
sequences of proper priors approaching the target improper priors
\eqref{improper-prior-alpha} and \eqref{improper-prior-sigma^2}:
\begin{equation}
p_\alpha(\alpha;h_\alpha)=\frac{1}{2h_\alpha} I_{(-h_\alpha,h_\alpha)}(\alpha)
\label{proper-prior-alpha}
\end{equation}
where $h_\alpha \to \infty$
and
\begin{equation}
p_{\eta}(\eta_{\bm{\gamma}};h_\eta) 
 =\frac{1}
{\int_{1/h_\eta}^{h_\eta} t^{-1}dt}
\eta_{\bm{\gamma}}^{-1}I_{(1/h_\eta, \ h_\eta)}(\eta_{\bm{\gamma}})
=
\frac{\eta_{\bm{\gamma}}^{-1}I_{(1/h_\eta, \ h_\eta)}(\eta_{\bm{\gamma}})}{2\log h_\eta} 
\label{proper-prior-sigma^2}
\end{equation}
where $h_\eta \to \infty$. 
See the beginning of Section \ref{sec:marginal} for details of the
justification. 

Next we give a sequence of proper conditional priors on
$\bm{\beta}_{\bm{\gamma}}$ given $\eta_{\bm{\gamma}}$, which approach
an improper conditional prior on $\bm{\beta}_{\bm{\gamma}}$ given $\eta_{\bm{\gamma}}$:
\begin{equation}\label{g-prior-1}
p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g) 
=
\int_0^\infty p_g(g|\nu;h_g) 
\phi_{q_{\bm{\gamma}}}(\bm{\beta}_{\bm{\gamma}}| \bm{0},
g\eta_{\bm{\gamma}}^{-1}(\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}})^{-1})
dg
\end{equation}
where $\nu$ and $h_g$ are non-random positive parameters. Further
\begin{equation}\label{g-prior-2}
 p_g(g|\nu;h_g) = \{\nu/2\}h_g^{-\nu/2}g^{\nu/2-1}I_{[0,h_g]}(g),
\end{equation}
and $\phi_q(\cdot| \bm{\mu}, \bm{\Sigma})$
denotes the $q$-variate Gaussian density with mean vector $\bm{\mu}$ and
covariance matrix $\bm{\Sigma}$.
The prior \eqref{g-prior-1} clearly has a hierarchical structure and 
it can be interpreted a scale mixture of Zellner's $g$-priors.
Similar priors have been considered by \cite{Liang-etal-2008} and
\cite{Maruyama-George-2010} under the Gaussian linear regression setup.
See Sub-Section \ref{sec:g-prior} below for a review of priors on $g$.
For any fixed $h_g>0$ and $\nu>0$, 
the prior $p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g)$ 
is clearly a proper probability density, that is,
\begin{equation*}
\int_{\mathcal{R}^{q_{\bm{\gamma}}}} p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g) d\bm{\beta}_{\bm{\gamma}}=1.
\end{equation*}
When $h_g \to \infty$ and $0<\nu<q_{\bm{\gamma}}$, 
the limit of a variant of 
$p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g)$ is given by
\begin{equation}\label{eq:limit.variant}
\begin{split}
& p^I_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\}) \\ 
&=
\lim_{h_g \to \infty} \left\{\frac{h_g^{\nu/2}}{\nu/2}\right\}
p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g) \\
& = \int_0^{\infty} g^{\nu/2-1} 
\frac{|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2}\eta_{\bm{\gamma}}^{q_{\bm{\gamma}}/2}}{(2\pi)^{q_{\bm{\gamma}}/2}g^{q_{\bm{\gamma}}/2}}
\exp\left(-\frac{\eta_{\bm{\gamma}}}{2g}\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}}\right)dg \\
&= \frac{\Gamma(\{q_{\bm{\gamma}}-\nu\}/2)}{2^{\nu/2}\pi^{q_{\bm{\gamma}}/2}}|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2} 
(\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}})^{-(q_{\bm{\gamma}}-\nu)/2}\eta_{\bm{\gamma}}^{\nu/2}.
\end{split}
\end{equation}
In summary, 
the proper  prior joint density under $\mathcal{M}_{\bm{\gamma}}$, 
which we will use in this paper, is given by
\begin{equation} \label{proper-joint}
p_{\bm{\gamma}}(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}}|\nu;\bm{h})
=p_\alpha(\alpha;h_\alpha) 
p_{\eta}(\eta_{\bm{\gamma}};h_\eta) 
p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\};h_g)
\end{equation}
with $\bm{h}=(h_\alpha,h_\eta,h_g)'$, which clearly satisfies
\begin{equation*}
\int_{-\infty}^\infty \int_{\mathcal{R}^{q_{\bm{\gamma}}}} \int_0^\infty 
p_{\bm{\gamma}}(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}}|\nu;\bm{h}) d\alpha d\bm{\beta}_{\bm{\gamma}} d\eta_{\bm{\gamma}}
=1,
\end{equation*}
for any fixed $\bm{h}$ and any $\nu>0$.
In Section \ref{sec:marginal}, we will use the improper joint density 
of $\alpha$, $\bm{\beta}_{\bm{\gamma}}$ and $\eta_{\bm{\gamma}}$ given by
\begin{equation} \label{improper-joint}
\begin{split}
& p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu) 
 =p^I_\alpha(\alpha)p^I_{\eta}(\eta_{\bm{\gamma}})
p^I_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\{\eta_{\bm{\gamma}},\nu\}) \\
&= \lim_{\bm{h}\to \bm{\infty}}V(\bm{h},\nu)
p_{\bm{\gamma}}(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}}|\nu;\bm{h})\\
& = 
\frac{\Gamma(\{q_{\bm{\gamma}}-\nu\}/2)}{2^{\nu/2}\pi^{q_{\bm{\gamma}}/2}}|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2} 
 (\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}})^{-(q_{\bm{\gamma}}-\nu)/2}\eta_{\bm{\gamma}}^{\nu/2-1},
\end{split}
\end{equation}
where $0<\nu<q_{\bm{\gamma}}$ and
\begin{equation}
 V(\bm{h},\nu)=\{2h_\alpha\}\{2\log h_{\eta}\}\{(2/\nu)h_g^{\nu/2}\}
= (8/\nu)h_\alpha \{\log h_{\eta}\} h_g^{\nu/2}.
\end{equation}
Hence, in a certain sense, the larger $\nu$ is (as long as $q_\bm{\gamma}-\nu>0$), 
the more objective the joint prior $p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu)$ is.

In this presentation of the improper joint density $p^I_{\bm{\gamma}}(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu)$, 
two facts,
\begin{enumerate}[{K}1.]
\item
$(\alpha,\bm{\beta}_{\bm{\gamma}})$ and $\eta_{\bm{\gamma}}$ are separable, \label{Key1}
\item 
the part depending on $\eta_{\bm{\gamma}}$ is given by the power function $\eta_{\bm{\gamma}}^{\nu/2-1}$, \label{Key2}
\end{enumerate}
will be the key for calculating the marginal density in the next section.
See also Remark \ref{rem:separable}.

If, in the above joint prior on $(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}})$, we make 
the change of variables, $\bm{\theta}_{\bm{\gamma}}=(\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}})^{1/2}\bm{\beta}_{\bm{\gamma}}$,
the joint prior of $(\alpha,\bm{\theta}_{\bm{\gamma}},\eta_{\bm{\gamma}})$ becomes
\begin{equation} \label{improper-joint-1}
 p_{\bm{\gamma}}^I(\alpha,\bm{\theta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu) 
 = 
\frac{\Gamma(\{q_{\bm{\gamma}}-\nu\}/2)}{2^{\nu/2}\pi^{q_{\bm{\gamma}}/2}}
\|\bm{\theta}_{\bm{\gamma}}\|^{-(q_{\bm{\gamma}}-\nu)}\eta_{\bm{\gamma}}^{\nu/2-1}.
\end{equation}
As noted in Section \ref{sec:intro},
the part depending on $\bm{\theta}_{\bm{\gamma}}$, $\|\bm{\theta}_{\bm{\gamma}}\|^{-(q_{\bm{\gamma}}-\nu)}$ for $0 <\nu<\min(2,q_{\bm{\gamma}})$,
is known as a sub-harmonic function, 
that is,
\begin{equation*}
 \sum_{i=1}^{q_{\bm{\gamma}}} \frac{\partial^2}{\partial \theta_i^2}
\|\bm{\theta}_{\bm{\gamma}}\|^{-(q_{\bm{\gamma}}-\nu)} 
 = (q_{\bm{\gamma}}-\nu)(2-\nu)\|\bm{\theta}_{\bm{\gamma}}\|^{-(q_{\bm{\gamma}}-\nu)-2}>0.
\end{equation*}

\subsection{Review of priors on $g$} \label{sec:g-prior}
As noted above, the prior given by \eqref{g-prior-1}
is a scale mixture of Zellner's $g$-priors.
Actually the original Zellner's $g$-priors were used
for the Gaussian linear regression setup and historically 
the hyperparameter $g$ has been a priori fixed or somehow estimated.
The first paper to effectively use a prior on $g$
was \cite{Zellner-Siow-1980}; 
they stated things in terms of multivariate Cauchy densities, 
which can always be expressed as a mixture of $g$-priors.
Here we review the prior on $g$, the second stage of $g$-priors 
for Gaussian linear regression. We hope that it makes the positioning
of our prior on $g$, \eqref{g-prior-2}, which is for SS 
errors, clearer.

As a generalization of $p_g(g|\nu;h_g)$ given by \eqref{g-prior-2},
consider the prior on $g$,
\begin{equation} \label{nu_a}
p_g(g|\{\nu,k\};h_g)=
\frac{g^{\nu/2-1}(1+g^{-1})^{-k/2}I_{[0,h_g]}(g)}{V_g(\{\nu,k\};h_g)},
\end{equation}
(so $ p_g(g|\{\nu,0\};h_g)=p_g(g|\nu;h_g)$).
In \eqref{nu_a}, $k\geq 0$, $\nu>-k$, $h_g>0$ and 
\begin{equation*}
  V_g(\{\nu,k\};h_g)=\int_0^{h_g}g^{\nu/2-1}(1+g^{-1})^{-k/2}dg.
\end{equation*}
Note that it is improper at $g=0$ when $ \nu\leq -k$.
As we will see in Section \ref{sec:marginal} and \ref{sec:Laplace},
Bayes factors are not well-defined when the prior on $g$ is improper
at $g=0$ and that is why $ \nu > -k$ is assumed.

When $k>0$ and $ -k<\nu<0$, without truncating by $h_g$
the prior of $g$ is proper under the support $g\in(0,\infty)$.
In this case, the normalizing constant becomes
 \begin{equation*}
  V_g(\{\nu,k\};\infty)=B(-\nu/2,\{\nu+k\}/2).
 \end{equation*}
Further even when $0\leq\nu< q_{\bm{\gamma}}$ (non-negative $\nu$ implies impropriety
at $g=\infty$),
the Bayes factor under $\mathcal{M}_{\bm{\gamma}}$ is well-defined as shown
in Section \ref{sec:marginal}. 

As in \eqref{g-prior-1} and \eqref{eq:limit.variant}, 
we are considering the (improper) 
case where $0<\nu<q_{\bm{\gamma}}$ and $k=0$.
\cite{Liang-etal-2008} considered the (proper) case where $-2<\nu<0$
and $k=2-\nu$.
\cite{Guo-Speckman-2009} and \cite{Celeux-Anbari-Marin-Robert-2011}
considered the (improper) case where $\nu=0$ and $k=2$.
\begin{remark}
The noteworthy difference between $k>0$ and $k=0$ is that
the (proper or improper) prior on $\bm{\beta}_{\bm{\gamma}}$ given $\eta_{\bm{\gamma}}$, 
\begin{equation}\label{bbe|eta}
\begin{split}
&p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\eta_{\bm{\gamma}},\{\nu,k\}) \\
&\stackrel{\mathrm{def.}}{=}\lim_{h_g\to\infty}
V_g(\{\nu,k\};h_g)
\int_0^{h_g} p_g(g|\{\nu,k\};h_g)\phi_{q_{\bm{\gamma}}}(\bm{\beta}_{\bm{\gamma}}| \bm{0},
g\eta_{\bm{\gamma}}^{-1} (\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}})^{-1})dg \\
&=\int_0^\infty \frac{g^{\nu/2-1}}{(1+g^{-1})^{k/2}}
\frac{\eta_{\bm{\gamma}}^{q_{\bm{\gamma}}/2}|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2}}{(2\pi g)^{q_{\bm{\gamma}}/2}}
\exp\left(-\frac{\eta_{\bm{\gamma}}}{2g}\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}}\right)dg 
\end{split}
\end{equation}
is not separable for any $k>0$. As will be seen in Section \ref{sec:marginal},
the separability of the prior is the key for calculating 
the marginal density for SS  error models other than Gaussian.
\end{remark}

\begin{remark}
Here we discuss objectivity (or at least non-subjectivity)
of the prior in terms of hyper-parameters of the prior on $g$.
Consider the (proper or improper) prior on
$\bm{\beta}_{\bm{\gamma}}$ given $\eta_{\bm{\gamma}}$, which is given by \eqref{bbe|eta}.
In order to obtain the asymptotic behavior of the density as
$ \eta_{\bm{\gamma}}^{1/2}\|\bm{\beta}_{\bm{\gamma}}\|\to\infty$,
we appeal to the Tauberian theorem for the Laplace transform 
(see \cite{geluk-dehaan-1987}).
Since $ (1+g^{-1})^{k/2}\to 1$ for any $k $ as $g\to\infty$,
\begin{equation} \label{eq:asymp.beta}
 \lim_{\eta_{\bm{\gamma}}^{1/2}\|\bm{\beta}_{\bm{\gamma}}\|\to\infty}
\frac{p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\eta_{\bm{\gamma}},\{\nu,k\})}
{\left\{\eta_{\bm{\gamma}}\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}
\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}}\right\}^{-\frac{q_{\bm{\gamma}}-\nu}{2}}
\eta_{\bm{\gamma}}^{\frac{q_{\bm{\gamma}}}{2}}} 
=
\frac{\Gamma(\{q_{\bm{\gamma}}-\nu\}/2)}{2^{\nu/2}\pi^{q_{\bm{\gamma}}/2}}|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2}
\end{equation}
when $\nu<q_{\bm{\gamma}}$.
Hence the asymptotic order of \eqref{eq:asymp.beta} is 
the same as \eqref{eq:limit.variant}.
The larger $\nu(<q_{\bm{\gamma}})$ is, the more objective the prior
$ p_{\bm{\beta}|\eta}(\bm{\beta}_{\bm{\gamma}}|\eta_{\bm{\gamma}},\{\nu,k\})$ is.
\end{remark}

\section{Marginal density and Bayes factor under sub-harmonic priors}
\label{sec:marginal}
In this section we derive the marginal density under each submodel
and the Bayes factor for comparing
each $ \mathcal{M}_{\bm{\gamma}}$ to the full model $\mathcal{M}_F$.
The marginal density of $\bm{y}$ under
$\mathcal{M}_{\bm{\gamma}} $,
is given by
\begin{equation}
\begin{split}
 m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h})= &
\int_{-\infty}^{\infty} \int_{R^{q_{\bm{\gamma}}}}
 \int_{0}^{\infty}
\eta_{\bm{\gamma}}^{n/2}f(\eta_{\bm{\gamma}}\|\bm{y}-\alpha \bm{1}_n -\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}}\|^2) \\
& \qquad \times p_{\bm{\gamma}}(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}}|\nu;\bm{h})
d \alpha \, d \bm{\beta}_{\bm{\gamma}} \, d\eta_{\bm{\gamma}} , 
\end{split}
\label{full-marginal}
\end{equation}
where the proper joint prior 
$ p_{\bm{\gamma}}(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}}|\nu;\bm{h})$ is 
given by \eqref{proper-joint}.
In \eqref{full-marginal}, $\bm{h}$
does not depend on the submodel, but is the same in all models.
In the following, instead of $m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h})$ directly,
we consider the limit of the product of $m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h})$
and $V(\bm{h},\nu)$,
\begin{equation}\label{full-marginal-M} 
\begin{split}
& M_{\bm{\gamma}}(\bm{y}|\nu) = \lim_{\bm{h}\to\bm{\infty}}
V(\bm{h},\nu)
m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h}) \\
&=\int_{-\infty}^{\infty} \int_{R^{q_{\bm{\gamma}}}}
 \int_{0}^{\infty}
\eta_{\bm{\gamma}}^{n/2}f
\left(\eta_{\bm{\gamma}}\|\bm{y}-\alpha \bm{1}_n -\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}}\|^2\right) \\
&\qquad \qquad \qquad \qquad \qquad \qquad \times p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu)
d \alpha \, d \bm{\beta}_{\bm{\gamma}} \, d\eta_{\bm{\gamma}}, 
\end{split} 
\end{equation}
which is the marginal density of $\bm{y}$ with respect to 
the improper prior $p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu)$ 
given by \eqref{improper-joint}:
\begin{equation} 
 p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu) 
 = 
\frac{\Gamma(\{q_{\bm{\gamma}}-\nu\}/2)}{2^{\nu/2}\pi^{q_{\bm{\gamma}}/2}}
\frac{|\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}|^{1/2}\eta_{\bm{\gamma}}^{\nu/2-1}} 
{ (\bm{\beta}'_{\bm{\gamma}}\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}}\bm{\beta}_{\bm{\gamma}})^{(q_{\bm{\gamma}}-\nu)/2}}.
\end{equation}
The second equality in \eqref{full-marginal-M} 
follows from the monotone convergence theorem.
Since $V(\bm{h},\nu)$ does not depend on $\mathcal{M}_\bm{\gamma}$ if
we choose the same $\nu$ in all submodels, 
the Bayes factor 
$m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h})/m_F(\bm{y}|\nu;\bm{h})$ approaches 
$M_{\bm{\gamma}}(\bm{y}|\nu)/M_F(\bm{y}|\nu)$ as $\bm{h}\to\bm{\infty}$.
Hence the use of the improper joint prior 
is justified as long as 
$M_{\bm{\gamma}}(\bm{y}|\nu)/M_F(\bm{y}|\nu)$
is well-defined.
As remarked in Section \ref{sec:intro}, the null-model is not a possible model.
Since there is no $\bm{\beta}$ and hence no $h_g$,
$M_N(\bm{y}|\nu)/M_F(\bm{y}|\nu)$ is not well-defined.
And thus $ \pi_N=\mathrm{Pr}(\mathcal{M}_N)=0$ is assumed.
However, see Remark \ref{rem:null} for an alternative development which allows
$\mathcal{M}_N$ as a possible model.

Let $ M_{\bm{\gamma}}^G(\bm{y}|\nu) $ be the marginal density under $\mathcal{M}_{\bm{\gamma}}$
with standard Gaussian errors $\bm{\epsilon}_G$. 
Before proceeding with the entire calculation of the marginal density, 
$ M_{\bm{\gamma}}(\bm{y}|\nu) $, when $\bm{\epsilon}_\bm{\gamma}$ has a general SS  distribution,
we will provide a relationship between $ M_{\bm{\gamma}}(\bm{y}|\nu) $ and $ M^G_{\bm{\gamma}}(\bm{y}|\nu) $ as follows.
\begin{lemma}\label{relationship-G}
Let $\nu$ be between $0$ and $q_{\bm{\gamma}}$. 
Assume the existence of $E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]$. Then
\begin{equation} \label{MG}
M_{\bm{\gamma}}(\bm{y}|\nu)
=\frac{E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]}{E[\|\bm{\epsilon}_G \|^\nu]}
M_{\bm{\gamma}}^G(\bm{y}|\nu).
\end{equation}
\end{lemma}
\begin{proof}
 See Appendix.
\end{proof}
Hence $ M_{\bm{\gamma}}(\bm{y}|\nu)$ depends on the error distribution $ \bm{\epsilon}_{\bm{\gamma}}$
only through the $\nu$-th moment of $\bm{\epsilon}_{\bm{\gamma}}$, 
$ E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]$. 

\begin{remark}\label{rem:separable}
This relationship of Lemma \ref{relationship-G} 
remains true under more general separable priors
$\eta_{\bm{\gamma}}^{\nu/2-1}\pi(\alpha, \bm{\beta}_{\bm{\gamma}})$, which satisfies
K\ref{Key1} and K\ref{Key2}. Unfortunately we could not find any other priors, 
which lead to analytically tractable Bayes factors under Gaussian errors,
except for our sub-harmonic priors $ p_{\bm{\gamma}}^I(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu) $
given by \eqref{improper-joint}.
It follows that for any such separable prior, the distributional
robustness of the Bayesian variable selection procedure we develop
for our sub-harmonic prior will carry over.
That is, whether MCMC or another computational method is used,
it suffices to assume that the error distribution is Gaussian 
(with the same error variance) for all submodels, but the procedure
will be simultaneously valid for all SS  error 
distributions.
Our prior has the advantage of analytic tractability, 
and as we show below, good performance.
\end{remark}

We will make use of the following result which may be founded in \cite{Liang-etal-2008}.
\begin{lemma} \label{m-G}
Let $0<\nu<q_{\bm{\gamma}}$. 
Then
\begin{equation}\label{MGG}
M^G_{\bm{\gamma}}(\bm{y}|\nu)
 = \frac{n^{1/2}\Gamma(\{n-1\}/2)}{\|\bm{y}-\bar{y}\bm{1}_n\|^{n-1}\pi^{(n-1)/2}}
\int_0^{\infty} \frac{g^{\nu/2-1}(1+g)^{(n-q_{\bm{\gamma}}-1)/2}}
{\left\{g(1-R_{\bm{\gamma}}^2)+1 \right\}^{(n-1)/2}} \, dg, 
\end{equation}
where $R_{\bm{\gamma}}^2$ is the coefficient of determination under the submodel $\mathcal{M}_{\bm{\gamma}}$.
\end{lemma}
See equation (5) of \cite{Liang-etal-2008}.
Combining Lemmas \ref{relationship-G} and \ref{m-G}, 
we have the main result of this paper.
\begin{thm}\label{thm:main-BF}
Assume the full model $\mathcal{M}_F$ and the submodel $\mathcal{M}_\bm{\gamma}$ are
given by \eqref{full-model} and \eqref{submodel-gamma}, respectively.
Also assume their error terms, $\bm{\epsilon}_F$ and $\bm{\epsilon}_\bm{\gamma}$ have the same
SS  distribution \eqref{bep_sim_f} with mean zero and the identity
covariance matrix.
Let $0<\nu <q_{\bm{\gamma}}$.
Assume that the proper joint prior densities
of $(\alpha, \bm{\beta}_F, \eta_F)$ and $(\alpha, \bm{\beta}_{\bm{\gamma}}, \eta_{\bm{\gamma}})$
are given by \eqref{proper-joint} and assume also $E[\|\bm{\epsilon}_{\bm{\gamma}}\|^\nu]<\infty$.
Then, for $\mathcal{M}_\bm{\gamma}\neq\mathcal{M}_N$, the limit of the Bayes factor for comparing
each of $ \mathcal{M}_{\bm{\gamma}}$ to the full model $ \mathcal{M}_F$ is
given by
\begin{equation}\label{eq:main-thm-BF}
\mathrm{BF}_{\gamma:F}(\nu)
 = \lim_{\bm{h}\to\bm{\infty}}
\frac{m_{\bm{\gamma}}(\bm{y}|\nu;\bm{h})}{m_F(\bm{y}|\nu;\bm{h})} 
=\mathrm{BF}^G_{\gamma:F}(\nu)
\end{equation}
where
\begin{equation} 
\label{main-BF}
\mathrm{BF}^G_{\gamma:F}(\nu)
=
\frac{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-q_{\bm{\gamma}}-1}{2}}
\{g(1-R_{\bm{\gamma}}^2)+1 \}^{-\frac{n-1}{2}}\, dg}
{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-p-1}{2}}
\{g(1-R_F^2)+1 \}^{-\frac{n-1}{2}}\, dg}.
\end{equation}
\end{thm}
\begin{remark}
\cite{Maruyama-George-2010}
considered Bayesian variable selection under Gaussian errors.
They proposed Bayes factors with a simple analytic form
under generalized ridge-type priors.
The results heavily depend on 
special features of Gaussian distributions and hence
the extension or generalization of \cite{Maruyama-George-2010}
to the general SS  case,
may not be possible, or may not lead to analytically tractable procedures
which are distributionally robust to SS  error distributions.
\end{remark}

\begin{remark}\label{rem:null}
Expression \eqref{main-BF} again shows why the null model is not allowed 
as a possibility.
For the null model $\mathcal{M}_N$, $R^2_N=0$ so the numerator of \eqref{main-BF} 
is infinite, and hence so would be $ \mathrm{BF}^G_{N:F}(\nu)$.
This situation may be avoided at a slight cost in complexity and
in interpretability of the expressions.
The required alteration in the prior distributions (proper and improper)
is to treat the intercept parameter $\alpha$ as another $\bm{\beta}$,
(and not give it a ``uniform'' prior).
This results in replacing the improper prior in \eqref{improper-joint} by
\begin{equation*} 
 p_{\bm{\gamma}}^I(\check{\bm{\beta}}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu)  = 
\frac{\Gamma(\frac{q_{\bm{\gamma}}+1-\nu}{2})}{2^{\frac{\nu}{2}}\pi^{\frac{q_{\bm{\gamma}}+1}{2}}}
|\check{\bm{X}}'_{\bm{\gamma}}\check{\bm{X}}_{\bm{\gamma}}|^{\frac{1}{2}} 
(\check{\bm{\beta}}'_{\bm{\gamma}}\check{\bm{X}}'_{\bm{\gamma}}\check{\bm{X}}_{\bm{\gamma}}\check{\bm{\beta}}_{\bm{\gamma}})
^{-\frac{q_{\bm{\gamma}}+1-\nu}{2}}\eta_{\bm{\gamma}}^{\frac{\nu}{2}-1},
\end{equation*}
where $\check{\bm{\beta}}_{\bm{\gamma}}=(\alpha,\bm{\beta}'_{\bm{\gamma}})'$ 
and $\check{\bm{X}}_{\bm{\gamma}}=(\bm{1}_n | \bm{X}_{\bm{\gamma}})$.
Similarly the marginal distribution in \eqref{MGG}
and the Bayes factor given by \eqref{main-BF}
are replaced by
\begin{equation*}
\check{M}^G_{\bm{\gamma}}(\bm{y}|\nu)
 = \frac{\Gamma(n/2)}{\|\bm{y}\|^{n}\pi^{n/2}}
\int_0^{\infty} \frac{g^{\nu/2-1}(1+g)^{(n-q_{\bm{\gamma}}-1)/2}}
{\left\{g(1-\check{R}_{\bm{\gamma}}^2)+1 \right\}^{n/2}}\, dg, 
\end{equation*}
and
\begin{equation} \label{eq:BFcheck}
\check{\mathrm{BF}}^G_{\gamma:F}(\nu)
=
\frac{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-q_{\bm{\gamma}}-1}{2}}
\{g(1-\check{R}_{\bm{\gamma}}^2)+1 \}^{-\frac{n}{2}}
\, dg}{\int_0^{\infty} g^{\frac{\nu}{2}-1}(1+g)^{\frac{n-p-1}{2}}
\{g(1-\check{R}_F^2)+1 \}^{-\frac{n}{2}}
\, dg},
\end{equation}
where $ \check{R}_{\bm{\gamma}}^2=1-\mathrm{RSS}_{\bm{\gamma}}/\|\bm{y}\|^2$,
(the ``coefficient of determination'' of the model $\mathcal{M}_{\bm{\gamma}}$ relative to
the $0$-intercept model).
Hence with the substitution $R^2_{\bm{\gamma}} \to \check{R}^2_{\bm{\gamma}}$, 
$n-1 \to n$, $q_{\bm{\gamma}} \to q_{\bm{\gamma}}+1$, $\bm{y}-\bar{y}\bm{1}_n \to \bm{y}$, 
all expressions and results in the paper remains valid and 
the result (Corollary \ref{cor:consistency}) on model selection consistency 
in Section \ref{sec:consistency} holds also for the null model.
Clearly $ \check{R}^2_{\bm{\gamma}}$ is somewhat unusual, but
if model selection consistency under the null-model is desirable,
we can use $ \check{\mathrm{BF}}^G_{\gamma:F}(\nu) $.
Actually, under the Gaussian regression setup, 
\cite{Guo-Speckman-2009} and \cite{Celeux-Anbari-Marin-Robert-2011}
recommend use of the Bayes factor as a function of $ \check{R}^2_{\bm{\gamma}}$,
which just substitutes $g^{\nu/2-1}$ with $1/(1+g)$ in 
$ \check{\mathrm{BF}}^G_{\gamma:F}(\nu) $ given by \eqref{eq:BFcheck}.
\end{remark}

\begin{remark}\label{rem:coherent}
A collection of Bayes factors is called coherent if 
\begin{equation*}
 \mathrm{BF}_{\gamma_1:\gamma_2}=\mathrm{BF}_{\gamma_1:\gamma_0}\mathrm{BF}_{\gamma_0:\gamma_2}, 
\mbox{ and }\mathrm{BF}_{\gamma_1:\gamma_2}=1/\mathrm{BF}_{\gamma_2:\gamma_1},
\end{equation*}
for all $\gamma_1$ and $\gamma_2$ (see, e.g.~\cite{Robert-2007-BC}).
By \eqref{eq:main-thm-BF}, the Bayes factors corresponding
to our sub-harmonic priors are coherent 
(with the exception of those involving the null model $\mathcal{M}_N$),
which is why we require $\pi_N=0$.
Also with the adaptation of the alternative specification given in
Remark \ref{rem:null}, coherence holds for all Bayes factors including
those involving the null model $\mathcal{M}_N$.

As in \cite{Zellner-Siow-1980} and \cite{Liang-etal-2008}, 
the posterior probability of any model is an expression of the 
form \eqref{posterior-2} given by
\begin{equation*}
 \mathrm{Pr}(\mathcal{M}_{\bm{\gamma}}|\bm{y})
=\frac
{\pi_{\bm{\gamma}} \mathrm{BF}_{\gamma:F}}
{\sum_{\bm{\gamma}} \pi_{\bm{\gamma}} \mathrm{BF}_{\gamma:F}}.
\end{equation*}
In our development, we choose the full model $\mathcal{M}_F$ as the base model
rather than the null model $\mathcal{M}_N$ employing the encompassing approach of
\cite{Zellner-Siow-1980}. Without employing the adaptation of our prior
described in Remark \ref{rem:null},
the choice of the full model $\mathcal{M}_F$ as the base model, 
as opposed to the null model $\mathcal{M}_N$, is forced on us since $\pi_N=0$. 
\cite{Liang-etal-2008} argue that the null model is the 
superior choice as a base model in their setup (which also involves $g$-priors
or mixture thereof) due to incoherence of the Bayes factors 
if $\mathcal{M}_F$ is chosen as the base model (in their setup).
This incoherence arises in the setup of \cite{Liang-etal-2008} because
prior distribution on the full model $\mathcal{M}_F$ depends on each nested
alternative $\mathcal{M}_\gamma$. This incoherence is not a problem in our setup
since the choice of prior for each submodel depends only on the submodel,
and we have taken care that all relevant (conditional) posteriors are
well defined. 
As noted above, by \eqref{eq:main-thm-BF} our Bayes factors are coherent.
In fact our development (aside from eliminating the null model from the 
consideration) is very close in spirit to the null-based Bayes factors 
approach in \cite{Liang-etal-2008}. 
\end{remark}


\begin{remark}\label{rem:stein-effect}
By Theorem \ref{thm:main-BF},
even when there is no specific information about 
the error distribution of each model (other than spherical symmetry), 
but we assume they are all the same, 
it is not necessary to specify the exact form of the sampling density. 
It suffices to assume they are all Gaussian.
As far as we know, in the area of Bayesian variable selection with shrinkage priors,
the sampling density has been assumed to be Gaussian and
this kind of robustness result has not yet been studied.
Analogous robustness results have been derived by \cite{Cellier-et-al-1989},
\cite{Maruyama-2003b}, \cite{Maru-Straw-2005} and others
in the problem of estimating regression coefficients with the Stein effect.
\end{remark}



\subsection{BIC under spherically symmetric error distributions}
BIC (\cite{Schwarz-1978}) is a popular criterion for model selection. 
See e.g.~\cite{Hastie-etal-2009} Chapter 7. 
We will show in this subsection 
that BIC has a similar distributional robustness property to the
above Bayes model selection procedure.
In Section \ref{sec:Laplace}, we will develop Laplace approximations
to our Bayes factors which relate them to BIC.
In Section \ref{sec:consistency} we will show that both the BIC and our Bayes model 
selection procedures are consistent for the entire class of
SS  models.

BIC for the model $\mathcal{M}_{\bm{\gamma}}$ is defined as
\begin{equation} \label{BIC-1}
[\mathrm{BIC}]_\bm{\gamma}=-2\ln
\left\{ 
\max_{\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}}\eta_{\bm{\gamma}}^{n/2}
f\left(
\eta_{\bm{\gamma}}\|\bm{y}-\alpha \bm{1}_n -\bm{X}_{\bm{\gamma}} \bm{\beta}_{\bm{\gamma}} \|^2
\right)
n^{-q_{\bm{\gamma}}/2}\right\},
\end{equation}
and is derived by eliminating $O(1)$ terms
from the approximate marginal densities.
Here we denote 
\begin{equation}
 M_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})=\exp(-[\mathrm{BIC}]_\bm{\gamma}/2).
\end{equation}
In general, 
the maximization with respect to unknown parameters in \eqref{BIC-1}
is not always tractable. However
when $\bm{\epsilon}_{\bm{\gamma}} $ has a unimodal 
SS  distribution,
the maximization is achieved by $\hat{\alpha}=\bar{y}$, 
$\hat{\bm{\beta}}_{\bm{\gamma}}=(\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}})^{-1}\bm{X}'_{\bm{\gamma}}\bm{y}$, and 
\begin{equation}
 1/\hat{\eta}_{\bm{\gamma}}=c\| \bm{y} -\hat{\alpha}\bm{1}_n -\bm{X}_{\bm{\gamma}}\hat{\bm{\beta}}_{\bm{\gamma}}\|^2
=c \|\bm{y} -\bar{y}\bm{1}_n\|^2(1-R_{\bm{\gamma}}^2)
\end{equation}
where $c$ is the sole solution of 
\begin{equation} \label{c-gamma}
 n/2+cf'(c)/f(c)=0.
\end{equation}
Hence $ M_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})$ may be expressed as
\begin{equation}\label{MG-BIC}
M_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})
=\frac{c^{-n/2}f(c)}{c_G^{-n/2}f_G(c_G)}
M^G_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})
\end{equation}
where $ M^G_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})$ is $ M_{\bm{\gamma}}(\bm{y}|\mathrm{BIC}) $ 
with Gaussian errors, specifically
\begin{equation} \label{MGG-BIC}
\begin{split}
M^G_{\bm{\gamma}}(\bm{y}|\mathrm{BIC}) &=
c_G^{-n/2}f_G(c_G)\{ \|\bm{y} -\bar{y}\bm{1}_n\|^2(1-R_{\bm{\gamma}}^2)\}^{-n/2}
n^{-q_{\bm{\gamma}}/2} \\
&=  n^{-n/2}f_G(n)\{ \|\bm{y} -\bar{y}\bm{1}_n\|^2(1-R_{\bm{\gamma}}^2)\}^{-n/2}
n^{-q_{\bm{\gamma}}/2}
\end{split}
\end{equation}
(since $c_G$ is given by $n$).
Clearly \eqref{MG-BIC} and \eqref{MGG-BIC} correspond to \eqref{MG} and
\eqref{MGG}, respectively.
Hence we have the following result.
\begin{thm}\label{thm:BICrobust}
Assume the full model $\mathcal{M}_F$ and the submodel $\mathcal{M}_\bm{\gamma}$ are
given by \eqref{full-model} and \eqref{submodel-gamma}, respectively.
Also assume their error terms, $\bm{\epsilon}_F$ and $\bm{\epsilon}_\bm{\gamma}$ have a
unimodal SS  distribution \eqref{bep_sim_f} with the mean zero and the identity
covariance matrix.
Then the Bayes factor based on BIC 
for comparing each of $ \mathcal{M}_{\bm{\gamma}}$ to the full model $ \mathcal{M}_F$ is given by
\begin{equation} \label{main-BIC}
\mathrm{BF}_{\gamma:F}[\mathrm{BIC}]
=\frac{M_{\bm{\gamma}}(\bm{y}|\mathrm{BIC})}{M_F(\bm{y}|\mathrm{BIC})}
=\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]
\end{equation}
where $\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]$ is the BIC based Bayes factor under 
Gaussian errors,
\begin{equation}\label{BF-G-BIC}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]
= \left\{\frac{(1-R_{\bm{\gamma}}^2)^{-n}n^{-q_{\bm{\gamma}}}}{(1-R_F^2)^{-n}n^{-p}}\right\}^{1/2}.
\end{equation}
\end{thm}
Obviously \eqref{BF-G-BIC} corresponds to \eqref{main-BF}.
By Theorem \ref{thm:BICrobust}, 
the Bayes factor based on BIC is also independent of 
the error distribution provided each distribution is unimodal and is
the same for all models (c.f.~Theorem \ref{thm:main-BF}).
Note that $\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]$ is well defined if $\mathcal{M}_\bm{\gamma}=\mathcal{M}_N$.
\begin{remark}\label{rem:generalization}
 In an earlier version of this paper we developed the results in the more
general context wherein the SS  distribution of 
$\bm{\epsilon}_\bm{\gamma} \sim f_\bm{\gamma}(\|\bm{\epsilon}_\bm{\gamma}\|^2)$ could depend on $\bm{\gamma}$, i.e.~it
could be different for each submodel $\mathcal{M}_\bm{\gamma}$. 
So $\bm{\epsilon}_F\sim f_F(\|\bm{\epsilon}_F\|^2)$ as well.
All of the above results can be developed for the more general case. 
The only essential changes
are that \eqref{eq:main-thm-BF} in Theorem \ref{thm:main-BF} becomes
\begin{equation}\label{eq:main-thm-BF-1}
\mathrm{BF}_{\gamma:F}(\nu)
=\frac{E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]}{E[\|\bm{\epsilon}_F \|^\nu]}
\mathrm{BF}^G_{\gamma:F}(\nu)
\end{equation}
and that \eqref{main-BIC} in Theorem \ref{thm:BICrobust} becomes
\begin{equation}\label{main-BIC-1} 
\mathrm{BF}_{\gamma:F}[\mathrm{BIC}]
=\frac{c^{-n/2}_{\bm{\gamma}} f_{\bm{\gamma}}(c_{\bm{\gamma}})}{c^{-n/2}_F f_F(c_F)}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]
\end{equation}
where $c_\bm{\gamma}$ and $c_F$ are the corresponding solution of \eqref{c-gamma}
respectively. We  investigated the ranges of these ``correction terms'' in 
\eqref{eq:main-thm-BF-1} and \eqref{main-BIC-1} respectively
when $\bm{\epsilon}_\bm{\gamma}$ and $\bm{\epsilon}_F$ have possibly different SS  $t$-distributions
with at least $3$-degrees of freedom (so that the variances exist),
both analytically and numerically. We found that 
$\mathrm{BF}_{\gamma:F}(\nu)$ was independent of $n$ and reasonably stable for all $\nu$
in range $(0,1)$ but that stability was greater for $\nu$ close to $0$.
This trade-off between stability (favoring $\nu\approx 0$) and 
objectivity (favoring larger $\nu$) led us initially 
to prefer the midpoint of the allowable values in $(0,1)$, 
namely $ \nu=1/2$ as the default choice.
However the examples presented in Section \ref{sec:ex} indicate that
the performance of the method seems insensitive to the choice of $\nu$
in the range of $(0,1)$.
The correction factor for $\mathrm{BF}_{\gamma:F}[\mathrm{BIC}]$ on the other hand depends on
$n$ and is considerably less stable to changes in the distributions of
$\bm{\epsilon}_\bm{\gamma}$ than $\mathrm{BF}_{\gamma:F}(\nu)$.

It is interesting to note in connection with the above that
the correction term $E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]/E[\|\bm{\epsilon}_F \|^\nu]$
for $\mathrm{BF}_{\gamma:F}(\nu)$ approaches $1$ as $\nu\to 0$.
Hence choices of $\nu$ close to $0$ are essentially completely robust
to choice of SS  error distribution for the submodels.

Note also that if we force $\pi_\bm{\gamma}=0$ for all submodels such that
$q_\bm{\gamma}\leq 2$, then the allowable range of $\nu$ is $(0,3)$ and
hence $\nu=2$ becomes a possible choice. In this case, again,
the correction term $E[\|\bm{\epsilon}_{\bm{\gamma}} \|^\nu]/E[\|\bm{\epsilon}_F \|^\nu]=1$
regardless of the choice of error distributions, since
we have assumed the variance of each component is $1$. Hence, again, 
in this case, the Bayes factor is completely robust to choice of SS 
error distribution. Additionally, the case $\nu=2$
corresponds to the harmonic prior
\begin{equation*}
p^{I}_{\bm{\gamma}}(\alpha,\bm{\beta}_{\bm{\gamma}},\eta_{\bm{\gamma}}|\nu) \propto 
\|\bm{\theta}_{\bm{\gamma}}\|^{2-q_{\bm{\gamma}}}
\end{equation*}
where $\bm{\theta}_{\bm{\gamma}}=(\bm{X}'_{\bm{\gamma}}\bm{X}_{\bm{\gamma}})^{1/2}\bm{\beta}_{\bm{\gamma}}$
and $q_{\bm{\gamma}}\geq 3$.
It is well-known that the harmonic prior plays important roles 
in estimation problems with the Stein effect. 
See \cite{Maruyama-2003b} for the detail.
It is interesting to observe the additional advantage of 
the harmonic prior in the model choice problem.
See Section \ref{sec:cr} for some additional discussion of such priors.

\end{remark}




\section{The Laplace approximation of BF under Gaussian errors}
\label{sec:Laplace}
In Section \ref{sec:marginal}, we saw that the Bayes factor $\mathrm{BF}_{\gamma:F}(\nu) $
under SS  errors is equal to
$\mathrm{BF}^G_{\gamma:F}(\nu) $, which is the Bayes factor under Gaussian errors.
In this section, we consider the so-called Laplace approximation
of some Bayes factors under Gaussian errors.
We will approximate 
not only the function $\mathrm{BF}^G_{\gamma:F}(\nu) $ but also
Bayes factors with respect to more general priors where the prior on $g$ is
\eqref{nu_a};
\begin{equation}\label{prior-Laplace.sec}
p_g(g|\{\nu,k\};\infty) \propto g^{\nu/2-1}(1+g^{-1})^{-k/2}.
\end{equation}
When the same prior on $g$ is used for $ \mathcal{M}_{\bm{\gamma}}$ and $ \mathcal{M}_F$, 
improper choices of $\nu$ ($0\leq \nu<q_{\bm{\gamma}}$) as well as proper
choices of $\nu$ 
($-k<\nu<0$) are valid for use.
Under Gaussian errors, the Bayes factor for comparing
each of $ \mathcal{M}_{\bm{\gamma}}$ to $ \mathcal{M}_F$ 
is well-defined as
\begin{equation} 
\label{BF-1}
\begin{split}
&\mathrm{BF}^G_{\gamma:F}[\nu,k] \\
&=\frac{\int_0^{\infty} g^{\nu/2-1}(1+g^{-1})^{-k/2}(1+g)^{\frac{n-q_{\bm{\gamma}}-1}{2}}
\{g(1-R_{\bm{\gamma}}^2)+1 \}^{-\frac{n-1}{2}}\, dg}
{\int_0^{\infty} g^{\nu/2-1}(1+g^{-1})^{-a/2}(1+g)^{\frac{n-p-1}{2}}
\{g(1-R_F^2)+1 \}^{-\frac{n-1}{2}}\, dg}
\end{split}
\end{equation}
where $k\geq 0$, $-k<\nu<q_{\bm{\gamma}}$.

First we provide a summary of Laplace approximations to the integral 
based on \cite{Tierney-Kadane-1986}.
For integrals of the form 
\begin{equation*}
 \int_{-\infty}^\infty \exp(h(\tau,n))d\tau,
\end{equation*}
we make the use of the fully exponential Laplace approximation, based on
expanding a smooth unimodal function $h(\tau,n)$ in a Taylor series expansion
about $\hat{\tau}$, the mode of $h(\tau,n)$.
The Laplace approximation is given by 
\begin{equation}\label{Laplace1}
\lim_{n \to \infty} 
\frac
{\int_{-\infty}^\infty \exp(h(\tau,n))d\tau}
{(2\pi)^{1/2}\hat{\sigma}_h\exp(h(\hat{\tau},n))}
=1
\end{equation}
where
\begin{equation*}
 \hat{\sigma}_h=
\left\{-\frac{\partial^2 h(\tau,n)}{\partial\tau^2} \Big|_{\tau=\hat{\tau}} \right\}^{-1/2}.
\end{equation*}
In the following, we will use the symbol $f(n) \approx g(n) $ ($n \to \infty$)
if 
\begin{equation}
 \lim_{n \to \infty} \frac{f(n)}{g(n)}=1.
\end{equation}
Hence the approximation given by \eqref{Laplace1} is written as
\begin{equation}
\int_{-\infty}^\infty \exp(h(\tau,n))d\tau \approx 
(2\pi)^{1/2}\hat{\sigma}_h \exp(h(\hat{\tau},n)) , \ (n \to \infty).
\end{equation}

The next result gives approximations of the Bayes factor
\eqref{main-BF} in terms of the Bayes factor based on
BIC given in \eqref{BF-G-BIC}.
\begin{thm} \label{thm:Laplace}
Let the prior be given by \eqref{prior-Laplace.sec}. 
Assume that $\{\nu,k\}$ does not depend on $n$.
\begin{enumerate}
 \item Assume $-k<\nu<q$, and $0<r<1$. Then 
\begin{equation}\label{eq:thmLaplace1}
 \int_0^{\infty} \frac{g^{\nu/2-1}}{(1+g^{-1})^{k/2}}
\frac{(1+g)^{\frac{n-q-1}{2}}}{{(1+rg)^{\frac{n-1}{2}}}} dg \\
 \approx  
\left\{\frac{4\pi \varphi(q-\nu,r)}{n^{q-\nu}r^n}\right\}^{1/2},
\end{equation}
where $ \varphi(s,r)= rs^{s-1}\{(1/r-1)e\}^{-s}$.
\label{thm:Laplace:part1}
\item Assume that $-k<\nu<q_{\bm{\gamma}}$ and 
that $\{\nu,k\}$ does not depend on $\mathcal{M}_{\bm{\gamma}}$.
Assume also $0<R_\bm{\gamma}^2 \leq R^2_F<1$.
Then $ \mathrm{BF}^G_{\gamma:F}[\nu,k] \approx \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ where
\begin{equation}
 \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)  =
\left\{
\frac{\varphi(q_{\bm{\gamma}}-\nu,1-R_{\bm{\gamma}}^2)}{\varphi(p-\nu,1-R_F^2)}\right\}^{1/2}\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}] 
\end{equation} 
and $ \mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]$ is the $\mathrm{BIC}$ based alternative under Gaussian
errors
\label{thm:Laplace:part2}
\end{enumerate}
\end{thm}
\begin{proof}
See Appendix.
\end{proof}
Clearly the function $\varphi$ does not depend on $n$ and hence 
Theorem \ref{thm:Laplace} shows that $ \mathrm{BF}^G_{\gamma:F}[\nu,k]$ 
is asymptotically equivalent to BIC with a simple $O(1)$ correction function
depending $\nu$ 
as well as $\{p,q_{\bm{\gamma}}\}$ and the $R$-squares.
Although several fully Bayes factors for the variable selection problem
have been proposed in the literature, the relationship between the approximate Bayes
factors and naive BIC has not been shown to the authors' knowledge.
In this sense, while the main contributions in this paper 
are given in Section \ref{sec:marginal}, 
Theorem \ref{thm:Laplace} 
may be a practically useful contribution because of the simplicity of the 
approximate Bays factor.

In this section, we have considered general Bayes factors under 
Gaussian errors. Remember that $\mathrm{BF}_{\gamma:F}(\nu)$ the Bayes factor w.r.t.~sub-harmonic
priors, under SS  errors, is equal to $\mathrm{BF}^G_{\gamma:F}[\nu,0]$
for $0<\nu<1$. 
Under Gaussian errors, 
\cite{Liang-etal-2008} recommended
the use of $ \mathrm{BF}^G_{\gamma:F}[\nu,2-\nu]$ with $ -2<\nu<0$.
\cite{Guo-Speckman-2009} and \cite{Celeux-Anbari-Marin-Robert-2011}
recommended the use of $ \mathrm{BF}^G_{\gamma:F}[0,2]$.
By Theorem \ref{thm:Laplace},
these Bayes factors may be approximated as follows.
\begin{corollary}\label{cor:approxBF}
\begin{equation*}
\begin{split}
&   \mathrm{BF}_{\gamma:F}(\nu)=\mathrm{BF}^G_{\gamma:F}[\nu,0] \approx \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu) 
\mbox{ for  } 0<\nu<1. \quad \mbox{sub-harmonic prior}\\
&   \mathrm{BF}^G_{\gamma:F}[0,2] \approx \widetilde{\mathrm{BF}}^G_{\gamma:F}(0). \quad
\mbox{a version of \cite{Guo-Speckman-2009}}\\
&   \mathrm{BF}^G_{\gamma:F}[\nu,2-\nu] \approx \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu) 
\mbox{ for  } -2<\nu<0. \quad \mbox{\cite{Liang-etal-2008}}
\end{split}
\end{equation*}
\end{corollary}
In Section \ref{sec:ex}, we will see how 
approximate Bayes factors $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$
work numerically and 
how sensitive they are to to the choice of $\nu$.

\section{Model selection consistency} \label{sec:consistency}
In this section, we consider model selection consistency in the case
where $p$ is fixed and as $n$ approaches infinity.
Let $\mathcal{M}_T$ be the true model,
\begin{equation*}
 \bm{y}=\alpha_T \bm{1}_n+\bm{X}_T\bm{\beta}_T+\sigma_T \bm{\epsilon}.
\end{equation*}
Consistency for model choice is defined as
\begin{equation*}
\mathop{\mathrm{plim}}_{n \to \infty}
\mathrm{Pr}(\mathcal{M}_T |\bm{y})=1, 
\end{equation*}
where plim denotes convergence in probability and the probability
distribution is the sampling distribution under the true model
$\mathcal{M}_T$.
We will show that Bayes factors considered in the previous sections
has a model selection consistency under generally SS  errors.
The consistency property is clearly equivalent to
\begin{equation} \label{equiv-consistency}
\mathop{\mathrm{plim}}_{n \to \infty}
\mathrm{BF}_{\gamma:T}=
\mathop{\mathrm{plim}}_{n \to \infty}
\frac{\mathrm{BF}_{\gamma:F}}
{\mathrm{BF}_{T:F}}
=0  \quad \forall \bm{\gamma} \neq T.
\end{equation}
For model selection consistency, we make the following assumptions;
\begin{enumerate}[{A}1.]
\item  
$U_n=\|\bm{\epsilon}\|^2/n$ is bounded 
in probability from below and from above, that is,
for any $c>0$ and any positive integer $n$, 
there exists an $M$ such that
\begin{equation*}
\textstyle{\mathrm{Pr}\left(M^{-1}< U_n < M\right) >1-c}.
\end{equation*} \label{AS1} 
\item 
The limit of the correlation matrix of $x_1,\dots,x_p$, 
$\lim_{n \to \infty}\bm{X}'_F\bm{X}_F/n$, exists and is positive definite. 
\label{AS2}
\end{enumerate}
A\ref{AS1} seems more general than necessary.
It appears that, by the law of large numbers,
$U_n$ ought to converge to $1$ in probability, but
this is not necessarily true if the error distribution is not Gaussian.
In the case of a scale mixture of Gaussians, $U_n$ approaches, in law,
a random variable $\xi$ which has the distribution of the mixing variable
of the variance. Even when the error distribution is not a scale
mixture of Gaussians, A\ref{AS1} appears to be a reasonable
and minimal assumption.
A\ref{AS2} is the standard assumption which also appears in \cite{Knight-Fu-2000}
and \cite{Zou-2006}. 
Under these mild assumptions, 
we have following preliminary results for proving the consistency.

\begin{lemma} \label{lem:prel}
Assume A\ref{AS1} and A\ref{AS2}. 
\begin{enumerate}
 \item \label{1:lem:prel}
Assume $\mathcal{M}_\bm{\gamma}\neq\mathcal{M}_N$.
For any $0<k<1$ and any positive integer $n$, 
there exists a $c_1(\bm{\gamma}, k)>2$ such that
\begin{equation}
 \mathrm{Pr}\left(\frac{1}{c_1(\bm{\gamma}, k)} < R^2_{\bm{\gamma}}
<1-\frac{1}{c_1(\bm{\gamma}, k)}\right)>1-k.
\end{equation}
\item \label{2:lem:prel}
Let $\bm{\gamma} \supsetneq T$. 
Then $ (1-R_T^2)/(1-R_{\bm{\gamma}}^2)\geq 1$. Further 
for any $0<k<1$ and any positive integer $n$, 
there exists a $c_2(\bm{\gamma},T,k)>0 $ such that
\begin{equation}
 \mathrm{Pr}\left(1 \leq
\left(\frac{1-R_T^2}{1-R_{\bm{\gamma}}^2}\right)^n<
1+c_2(\bm{\gamma},T,k)\right)>1-k.
\end{equation}
\item \label{3:lem:prel}
Let $\bm{\gamma} \nsupseteq T$. Then 
for any $0<k<1$ and any positive integer $n$, 
there exists a $c_3(\bm{\gamma},T,k)>1$ such that
\begin{equation}
\mathrm{Pr}\left(\frac{1-R_T^2}{1-R_{\bm{\gamma}}^2}<1-\frac{1}{c_3(\bm{\gamma},T,k)} \right)
>1-k.
\end{equation}
\end{enumerate}
\end{lemma}
\begin{proof}
 See Appendix.
\end{proof}
First we give a consistency result on BIC.
\begin{thm}\label{thm:BICconsistency}
Assume A\ref{AS1} and A\ref{AS2}. 
The Bayes factor based on BIC under Gaussian errors
\begin{equation*}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]
= \left\{\frac{(1-R_{\bm{\gamma}}^2)^{-n}n^{-q_{\bm{\gamma}}}}{(1-R_F^2)^{-n}n^{-p}}\right\}^{1/2}
\end{equation*}
is consistent for model selection under SS  errors (including $\mathcal{M}_\bm{\gamma}=\mathcal{M}_N$).
\end{thm}
\begin{proof}
we show that
\begin{equation} \label{main-plim}
\mathop{\mathrm{plim}}_{n \to \infty}
\frac{\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]}{\mathrm{BF}^G_{T:F}[\mathrm{BIC}]}=
\mathop{\mathrm{plim}}_{n \to \infty}
\left\{ n^{q_T-q_{\bm{\gamma}}} 
\left(\frac{1-R^2_T}{1-R^2_{\bm{\gamma}}}\right)^{n}
\right\}^{1/2}=0.
\end{equation}
Consider the following two situations:
\begin{enumerate}
\item 
$\bm{\gamma} \supsetneq T$: 
By part \ref{2:lem:prel} of Lemma \ref{lem:prel}, 
$\{(1-R^2_T)/(1-R^2_{\bm{\gamma}})\}^{n}$ is bounded in probability.
Since $ q_{\bm{\gamma}}> q_T$, \eqref{main-plim} is satisfied.
\item
$\bm{\gamma} \nsupseteq T$: 
\quad
By part \ref{3:lem:prel} of Lemma \ref{lem:prel}, 
$ (1-R^2_T)/(1-R^2_{\bm{\gamma}})$ 
is strictly less than $1$ in probability. 
Hence $ \{(1-R^2_T)/(1-R^2_{\bm{\gamma}})\}^n$
converges to zero in probability
exponentially fast with respect to $n$.
Therefore, no matter what value $ q_T-q_{\bm{\gamma}}$ takes,
\eqref{main-plim} is satisfied.
\end{enumerate}
These complete the proof.
\end{proof}
Note that in Theorem \ref{thm:BICconsistency}
 we do not exclude the null model $\mathcal{M}_N$ and hence
BIC has model selection consistency even when the null model is true.
When we consider consistency of the Bayes factors
treated in the previous sections, 
$\mathrm{BF}_{\gamma:F}(\nu) $, $\mathrm{BF}^G_{\gamma:F}[\nu,k] $, $ \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$,
we have to exclude the null model $\mathcal{M}_N$, but they all still have
model selection consistency among non-null models.
\begin{corollary} \label{cor:consistency}
Assume A\ref{AS1} and A\ref{AS2}. $\{\nu,k\}$ is assumed independent of $n$ and
$\mathcal{M}_\bm{\gamma}$. 
Assume also $\mathcal{M}_N$ is excluded from possible models.
Then 
\begin{enumerate}
\item $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ for $ \nu<1$
is consistent for model selection under SS  errors.
\item $\mathrm{BF}^G_{\gamma:F}[\nu,k]$ for $-k<\nu<1$ 
is consistent for model selection under SS  errors.
\item 
$\mathrm{BF}_{\gamma:F}(\nu)$ for $0<\nu<1$ 
is consistent for model selection under SS  errors.
\end{enumerate}
\end{corollary}
\begin{proof}
By part \ref{1:lem:prel} of Lemma \ref{lem:prel}, when 
$ \mathcal{M}_\bm{\gamma}\neq\mathcal{M}_N$, both $R^2_\bm{\gamma}$ and $R^2_F$ 
are positive and strictly less than 1 
with probability 1. Hence both $\varphi(q_\bm{\gamma}-\nu,1-R^2_\bm{\gamma})$ and
$\varphi(p-\nu,1-R^2_F)$ 
where
\begin{equation*}
  \varphi(s,r)= rs^{s-1}\{(1/r-1)e\}^{-s}
\end{equation*}
are positive and bounded from above
with probability 1 provided $\nu<1$ and $\nu$ is independent of $n$ and $\mathcal{M}_\bm{\gamma}$.
(On the other hand, since $R_N^2 \equiv 0$, $\varphi(q_\bm{\gamma}-\nu,1-R^2_N)$
is not defined.)
As in Theorem \ref{thm:Laplace}, 
\begin{equation*}
 \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)  =
\left\{
\frac{\varphi(q_{\bm{\gamma}}-\nu,1-R_{\bm{\gamma}}^2)}{\varphi(p-\nu,1-R_F^2)}\right\}^{1/2}
\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}] .
\end{equation*} 
Hence consistency of $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ follows from
consistency of BIC.

Further as $n\to\infty$, we have
\begin{equation*}
 \mathrm{BF}^G_{\gamma:F}[\nu,k] \approx \widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)
\end{equation*}
by Theorem \ref{thm:Laplace}
provided $-k<\nu<1$ and $\{\nu,k\}$ are independent of $n$ and $\mathcal{M}_\bm{\gamma}$.
Hence consistency of $\mathrm{BF}^G_{\gamma:F}[\nu,k]$ follows from
consistency of $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$.

Remember that $\mathrm{BF}_{\gamma:F}(\nu)$, the Bayes factor w.r.t.~sub-harmonic
priors, under SS  errors, is equal to $\mathrm{BF}^G_{\gamma:F}[\nu,0]$
for $0<\nu<1$. 
Hence consistency of $\mathrm{BF}_{\gamma:F}(\nu)$ follows from consistency of $\mathrm{BF}^G_{\gamma:F}[\nu,k]$.
\end{proof}
\begin{remark}
\cite{Liang-etal-2008} established model selection consistency 
for $\nu<0$ and $k=2-\nu$ for Gaussian errors. Corollary \ref{cor:consistency}
in conjunction with Theorem \ref{thm:Laplace}
extends their result to the entire class of SS  distributions for 
a broader class of $\nu$ and $k$.
Also as noted in Remark \ref{rem:null}, $\check{\mathrm{BF}}^G_{\gamma:F}(\nu)$ results in 
model selection consistency for all models including the null model
for all SS  distributions. Additionally a development along the
lines of Theorem \ref{thm:Laplace} allows an analogous BIC based approximation
to $\check{\mathrm{BF}}^G_{\gamma:F}(\nu)$.

This extension also allows an extension of model selection consistency
to all SS  error distributions for the method of 
\cite{Guo-Speckman-2009,Celeux-Anbari-Marin-Robert-2011} based on
$\check{R}^2_\bm{\gamma}$. Further as noted above, an alternative robust sub-harmonic
prior based method, the model selection consistency for these
$\check{R}^2_\bm{\gamma}$ based methods also apply for the null model.

It should be emphasized in each of the above cases that is is the Bayes
factor method developed for the Gaussian case that is shown to have
model selection consistency for the entire class of SS  error distributions.
These Gaussian based Bayes factors, however, are not Bayes factors
for error distributions which are not Gaussian, the sole exception 
being our robust Bayes factors which are based on separable priors 
in the sense described earlier, and which are simultaneously 
Bayes factors relative to the same prior.
\end{remark}

\begin{remark}
The issue of model selection consistency in our setup, is somewhat
complicated by the wide choice of possible error distributions.
If all errors are normally distributed, then under our assumptions A\ref{AS2} 
on the design matrix $\bm{X}_F$, 
imply that each $R_{\bm{\gamma}}^2$ approaches a constant,
and that $\|\bm{\epsilon}\|^2/n \to 1$.
If on the other hand, all models are variance mixtures of Gaussians
with mixture variance distributed as a positive random variable $\xi$,
then $\|\bm{\epsilon}\|^2/n \to \xi$ a random variable, and $R^2_{\bm{\gamma}}$
also approaches a random variable which is bounded above and below
in probability provided that $\xi$ is similarly bounded.

In general philosophical terms, it might be better to assume that
the sequence of error terms $\bm{\epsilon}=(\epsilon_1,\dots,\epsilon_n)'$ 
are exchangeable for all $n$.
By De finetti's Theorem, this would imply that the error terms
all have a variance mixture of normal distributions.
We have chosen a slightly weaker requirement on the sequence of error
distributions, namely, that $\|\bm{\epsilon}\|^2/n$ remains bounded above and
below in probability, which extracts the necessary limiting behavior
of the error terms to ensure consistency of model selection.
Interestingly, although we attain model selection consistency with 
these assumptions, it is not necessarily true that
$1=\mathrm{var}\epsilon_i=\mathrm{var}\xi$ is consistently estimated
by $\|\bm{\epsilon}\|^2/n$.
\end{remark}

\section{Examples}\label{sec:ex}
In this section, 
we provide illustrations of the method using both simulated and real data.
In each example, we compare several different versions of the Laplace
approximated Bayes factors $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ and $\mathrm{BF}^G_{\gamma:F}[\mathrm{BIC}]$.
The values of $\nu$ are $-2,-1,0,0.5,0.95$. These choices correspond to our default
choice, $\nu=1/2$ and $\nu=0.95$ which also satisfies our robustness condition $0<\nu<1$.
The choice $\nu=0$ approximates $\mathrm{BF}^G_{\gamma:F}[0,2]$ of \cite{Guo-Speckman-2009} and
$\nu=-2$ and $\nu=-1$ approximates two choices of \cite{Liang-etal-2008}
as presented in Corollary \ref{cor:approxBF}.
\subsection{Simulation Studies}\label{sec:sim}
We compare numerical performance of our 
with BIC in a small simulation study.
We generated 16 possible correlated predictors ($p=16$) as follows:
\begin{equation*}
\begin{split}
& \overbrace{x_1, x_2}^{\mathrm{cor}=0.5},
\underbrace{x_3, x_4}_{\mathrm{cor}=-0.4},
 \overbrace{x_5, x_6}^{\mathrm{cor}=0.3},
 \underbrace{x_7, x_8}_{\mathrm{cor}=-0.2},
\overbrace{x_9, x_{10}}^{\mathrm{cor}=0.1}
\sim N(0,1) \\
& x_{11}, x_{12}, x_{13}, x_{14},
x_{15}, x_{16}  \sim N(0,1). 
\end{split}
\end{equation*}
Here ``cor'' denotes the correlation of two Gaussian random variables.
Also
$(x_1,x_2)$, $(x_3,x_4)$, $(x_5,x_6)$, $(x_7,x_8)$, $(x_9,x_{10})$,
$x_{11},x_{12}, x_{13}, x_{14}, x_{15},x_{16}$ are assumed to be independent.
After generating pseudo random $x_1,\dots,x_{16}$, we centered and scaled them
as noted in Section \ref{sec:intro}.
We set $n=30$ and consider 4 cases where the true predictors are
\begin{description}
 \item[$q_{T}=16$] \
$x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8, x_9, x_{10},
x_{11}, x_{12}, x_{13}, x_{14}, x_{15}, x_{16}$ 
 \item[$q_{T}=12$] \
$x_1, x_2, x_3, x_4, x_5, x_6, x_7, x_8,  x_{9}, x_{10}, x_{11}, x_{12}$ 
 \item[$q_{T}=8$] \ \, 
$x_1, x_2, \qquad \quad x_5, x_6, \qquad \quad x_9, x_{10}, x_{11} , x_{12}$ 
 \item[$q_{T}=4$] \  \,
 $x_1, x_2, \qquad \quad x_5, x_6$ 
\end{description}
(where $q_T$ denotes the number of true predictors) and the true model
is given by
\begin{equation} \label{simu-true-model}
\bm{y}=\bm{1}_{30} + 2 \sum_{i \in \{\mbox{true}\}}\bm{x}_i+ 
\sigma \times
\begin{cases}
 N_{30}(\bm{0},\bm{I}_{30}), \\ \mbox{Multi-}t(\bm{0},\bm{I}_{30};3,30),
\end{cases}
\end{equation}
with $\sigma=0.5,\ 1, \ 2$.
Tables \ref{normal} and \ref{multi-t} show 
how often the true model ranks first and how often 
it is in the top 3 among $2^{16}-1$ candidates
when the number of replicates is $N=200$. 
The error distributions are Gaussian (Table \ref{normal}) and multivariate-$t$
with 3 degrees of freedom (Table \ref{multi-t}).
For the case of normally distributed errors (Table \ref{normal}),
the Bayes factor methods performed well and stably for $\sigma=0.5$ and
$\sigma=1$ and did reasonably well for $\sigma=2$ for the smaller
true models ($q_T=4,8$). BIC seemed, generally, to have a preference
for larger models, and performed much less well than the Bayes factor method
for $\sigma=0.5$ and $\sigma=1$ for models of smaller size ($q_T=4,8,12$)
For $\sigma=2$, BIC did substantially better than BF for the largest
model ($q_T=16$) and somewhat better for $q_T=12$.
Performance of $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ seemed relatively insensitive
to the choice of $\nu$. 
When $q_T \neq 16$, the choice of $\nu$ makes little difference. 
But when $q_T = 16$, positive $\nu=(0.5,0.95)$ seems
to perform better especially for larger $\sigma$.

Interestingly, for the case of a multivariate-$t$ error distribution with $3$
degrees of freedom (the minimum so that a variance exists),
the numerical results were quite similar to those in the normal case
for both $\widetilde{\mathrm{BF}}^G_{\gamma:F}(\nu)$ and BIC, both quantitatively and qualitatively.
One possible aspect of the relative insensitivity of the results to 
choice of $\nu$ in heavy tailed case is the extension of model selection
consistency for the entire class of SS  errors to a broad class 
mixture of $g$-prior based methods given by Corollary \ref{cor:consistency}.

\begin{table} 
\setlength{\tabcolsep}{4pt}
\caption{Frequency of the true model (Gaussian error)}
\begin{center}
 \begin{tabular}{c|cc|cc|cc|cc} \toprule
$q_T$ & \multicolumn{2}{c}{16} & \multicolumn{2}{c}{12}
& \multicolumn{2}{c}{8} & \multicolumn{2}{c}{4} \\ \midrule 
rank & 1 & 1-3 & 1 & 1-3 &  1 & 1-3 &  1 & 1-3  \\
   \midrule \midrule
\multicolumn{3}{c}{$\sigma=0.5$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$ & 1.00 & 1.00 & 0.94 &  1.00 & 0.94 & 0.99 & 0.89 & 0.99  \\
$\widetilde{\mathrm{BF}}^G(0.5)$ & 1.00 & 1.00 & 0.95 & 1.00 & 0.94 & 0.99 & 0.88 & 0.99 \\
$\widetilde{\mathrm{BF}}^G(0)$  & 1.00 & 1.00 & 0.95 & 1.00 & 0.94 & 0.99 & 0.88 & 0.99   \\
$\widetilde{\mathrm{BF}}^G(-1)$ & 1.00 & 1.00 & 0.96 & 1.00 & 0.94 & 1.00 & 0.89 & 0.99  \\ 
$\widetilde{\mathrm{BF}}^G(-2)$ & 1.00 & 1.00 & 0.96 & 1.00 & 0.94 & 1.00 & 0.89 & 0.99  \\ 
\cmidrule(r){1-1} 
BIC & 1.00 & 1.00 & 0.44 & 0.62 & 0.28 & 0.46 & 0.20 & 0.35  \\
\midrule \midrule
\multicolumn{3}{c}{$\sigma=1$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$  & 0.85 & 0.92 & 0.87 & 0.99 & 0.86 & 0.97 & 0.76 & 0.95  \\
$\widetilde{\mathrm{BF}}^G(0.5)$  & 0.83 & 0.90 & 0.88 & 0.99 & 0.87 & 0.97 & 0.76 & 0.94 \\
$\widetilde{\mathrm{BF}}^G(0)$  & 0.80 & 0.87 & 0.89 & 1.00 & 0.87 & 0.97 & 0.75 & 0.94  \\
$\widetilde{\mathrm{BF}}^G(-1)$  & 0.73 & 0.81 & 0.89 & 1.00 & 0.88 & 0.97 & 0.74 & 0.94  \\
$\widetilde{\mathrm{BF}}^G(-2)$  & 0.55 & 0.72 & 0.90 & 1.00 & 0.89 & 0.98 & 0.76 & 0.95  \\
\cmidrule(r){1-1} 
BIC & 1.00 & 1.00 & 0.44 & 0.62 & 0.28 & 0.46 & 0.20 & 0.35  \\
\midrule \midrule
\multicolumn{3}{c}{$\sigma=2$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$ & 0.06 & 0.11 & 0.26 & 0.42 & 0.50 & 0.74 & 0.51 & 0.73 \\
$\widetilde{\mathrm{BF}}^G(0.5)$ & 0.05 & 0.10 & 0.25 & 0.41 & 0.51 & 0.74 & 0.50 & 0.72 \\
$\widetilde{\mathrm{BF}}^G(0)$ & 0.05 & 0.10 & 0.24 & 0.41 & 0.52 & 0.73 & 0.49 & 0.72 \\
$\widetilde{\mathrm{BF}}^G(-1)$ & 0.04 & 0.06 & 0.22 & 0.39 & 0.43 & 0.74 & 0.48 & 0.72 \\
 $\widetilde{\mathrm{BF}}^G(-2)$ & 0.02 & 0.03 & 0.17 & 0.32 & 0.47 & 0.72 & 0.43 & 0.71 \\
\cmidrule(r){1-1} 
BIC & 0.62 & 0.77 & 0.31 & 0.48 & 0.24 & 0.40 & 0.19 & 0.34 \\
\bottomrule
 \end{tabular}
\end{center}
\label{normal}
\end{table}%

\begin{table} 
\setlength{\tabcolsep}{4pt}
\caption{Frequency of the true model (multi-$t$ error)}
\begin{center}
 \begin{tabular}{c|cc|cc|cc|cc} \toprule
$q_T$ & \multicolumn{2}{c}{16} & \multicolumn{2}{c}{12}
& \multicolumn{2}{c}{8} & \multicolumn{2}{c}{4} \\ \midrule 
rank & 1 & 1-3 & 1 & 1-3 &  1 & 1-3 &  1 & 1-3  \\
   \midrule \midrule
\multicolumn{3}{c}{$\sigma=0.5$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$ & 0.94 & 0.95 & 0.92 & 0.98 & 0.85 & 0.97 & 0.84 & 0.96 \\
$\widetilde{\mathrm{BF}}^G(0.5)$ & 0.94 & 0.95 & 0.92 & 0.98 & 0.86 & 0.97 & 0.84 & 0.96 \\
$\widetilde{\mathrm{BF}}^G(0)$  & 0.93 & 0.95 & 0.93 & 0.98 & 0.86 & 0.97 & 0.84 & 0.96 \\
$\widetilde{\mathrm{BF}}^G(-1)$ & 0.93 & 0.94 & 0.93 & 0.98 & 0.87 & 0.97 & 0.83 & 0.96 \\
$\widetilde{\mathrm{BF}}^G(-2)$ & 0.91 & 0.92 & 0.95 & 0.98 & 0.89 & 0.98 & 0.84 & 0.96 \\
\cmidrule(r){1-1} 
BIC & 0.98 & 0.99 & 0.47 & 0.63 & 0.28 & 0.45 & 0.26 & 0.38 \\
\midrule \midrule
\multicolumn{3}{c}{$\sigma=1$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$  & 0.67 & 0.70 & 0.75 & 0.84 & 0.72 & 0.88 & 0.70 & 0.85 \\
$\widetilde{\mathrm{BF}}^G(0.5)$  & 0.64 & 0.68 & 0.76 & 0.84 & 0.71 & 0.88 & 0.70 & 0.85 \\
$\widetilde{\mathrm{BF}}^G(0)$  & 0.62 & 0.67 & 0.77 & 0.84 & 0.72 & 0.89 & 0.69 & 0.84 \\
$\widetilde{\mathrm{BF}}^G(-1)$  & 0.58 & 0.63 & 0.76 & 0.83 & 0.75 & 0.89 & 0.69 & 0.84 \\
$\widetilde{\mathrm{BF}}^G(-2)$  & 0.49 & 0.58 & 0.76 & 0.83 & 0.75 & 0.89 & 0.69 & 0.84 \\
\cmidrule(r){1-1} 
BIC & 0.89 & 0.93 & 0.45 & 0.60 & 0.27 & 0.44 & 0.25 & 0.37 \\
\midrule \midrule
\multicolumn{3}{c}{$\sigma=2$} & \multicolumn{6}{c}{} \\ 
\cmidrule(r){1-4} 
$\widetilde{\mathrm{BF}}^G(0.95)$ & 0.14 & 0.20 & 0.28 & 0.37 & 0.37 & 0.51 & 0.45 & 0.58 \\
$\widetilde{\mathrm{BF}}^G(0.5)$ & 0.14 & 0.18 & 0.29 & 0.36 & 0.37 & 0.51 & 0.45 & 0.57 \\
$\widetilde{\mathrm{BF}}^G(0)$ & 0.13 & 0.17 & 0.29 & 0.36 & 0.38 & 0.51 & 0.45 & 0.57 \\
$\widetilde{\mathrm{BF}}^G(-1)$ & 0.09 & 0.13 & 0.29 & 0.33 & 0.38 & 0.51 & 0.44 & 0.56 \\
$\widetilde{\mathrm{BF}}^G(-2)$ & 0.07 & 0.11 & 0.24 & 0.34 & 0.35 & 0.49 & 0.40 & 0.55 \\
\cmidrule(r){1-1} 
BIC & 0.47 & 0.59 &0.28 & 0.39 & 0.21 & 0.32 & 0.20 & 0.30 \\
\bottomrule
\end{tabular}
\end{center}
\label{multi-t}
\end{table}%

\subsection{Analysis of real data} \label{sec:real}
In this section, we apply our methods (approximate Bayes factor and BIC)
to Hald data set presented and analyzed in \cite{Casella-Moreno-2006}
and to the US Crime data set in \cite{Raftery-Madigan-Hoeting-1997}.
See those papers for detailed descriptions of the data sets.
Table \ref{table:hald} and \ref{table:uscrime} present posterior probabilities 
based on $\widetilde{\mathrm{BF}}^G(\nu)$ of the top three selected models 
(assuming equal prior probabilities on all models) for several different choices
of $\nu$ ($0.95,0.5,0,-1,-2$). 
BIC was also included in the study. In each case,
the first, second and third ranked choices based on $\widetilde{\mathrm{BF}}^G(\nu)$
were identical regardless of the choices of $\nu$.
Also in each case the top ranked submodel based on $\widetilde{\mathrm{BF}}^G(\nu)$
was regarded as reasonable in the earlier papers. As in the simulation study,
and as noted in several previous studies, BIC seems to choose bigger models.
In particular for the Hald data, the top choice $\{x_1,x_2\}$ agrees
with that of \cite{Casella-Moreno-2006} and also of \cite{Berger-Pericchi-1996}
and \cite{Draper-Smith-1998}.

For the US Crime data, our top ranked model agrees with that of the Occam's
window posterior in Table 2 of \cite{Raftery-Madigan-Hoeting-1997}.
Interestingly our second ranked model includes $x_{15}$ which does not
occur in any of \citeapos{Raftery-Madigan-Hoeting-1997} Occam's window model choices,
but which does occur in several models chosen by such classical methods
as Mallow's $C_p$, adjusted $R^2$, etc.~in their Table 1.

\begin{table}
\setlength{\tabcolsep}{4pt}
\caption{Hald data: posterior probabilities of top 3 selected models}
\begin{center}
\begin{tabular}{cccccccc} \toprule
& & & \multicolumn{5}{c}{$\widetilde{\mathrm{BF}}^G$} \\
&  & $\nu$ & 0.95 & 0.5 & 0 & -1 & -2 \\
1 & $\{1,2, \, \ \quad\}$ & & 0.66 & 0.63 & 0.61 & 0.57 & 0.54 \\
2 & $\{1,\qquad 4\}$ & & 0.16 & 0.17 & 0.17 & 0.18 & 0.20 \\
3 & $\{1,2,\quad 4\}$ & & 0.06 & 0.07 & 0.07 & 0.08 & 0.08 \\
\midrule
& & & BIC & \multicolumn{4}{c}{} \\
1 & $\{1,2, \, \ \quad\}$ & & 0.25 & \multicolumn{4}{c}{} \\
2 & $\{1,2,\quad 4\}$ & & 0.23 & \multicolumn{4}{c}{} \\
3 & $\{1,2,3, \ \ \}$ & & 0.23 & \multicolumn{4}{c}{} \\
\bottomrule
\end{tabular}
 \end{center}
\label{table:hald}
\end{table}

\begin{table}
\setlength{\tabcolsep}{4pt}
\caption{US crime data: posterior probabilities of top 3 selected models}
\begin{center}
\begin{tabular}{cccccccc} \toprule
& & & \multicolumn{5}{c}{$\widetilde{\mathrm{BF}}^G$} \\
& & $\nu$ & 0.95 & 0.5 & 0 & -1 & -2  \\
1 & $\{1,3,4, \quad 9,11, \quad \ 13,14 \quad \ \}$   & & 0.020 & 0.019 & 0.018 & 0.016  & 0.015 \\
2 & $\{1,3,4,\quad 9,11, \quad \ 13,14,15\}$ & & 0.018 & 0.018 & 0.017 & 0.015 & 0.014 \\
3 & $\{1,3,\quad 5,9,11, \quad \ 13,14 \quad \ \}$    & & 0.013 & 0.013 & 0.012 & 0.011 & 0.010 \\
\midrule
& & & BIC & \multicolumn{4}{c}{} \\
1 & $\{1,3,4,\quad 9,11, \quad \ 13,14,15\}$ & & 0.035 & \multicolumn{4}{c}{} \\
2 & $\{1,3,4, \quad 9,11, \quad \ 13,14 \quad \ \}$ & & 0.026 & \multicolumn{4}{c}{} \\
3 & $\{1,3,4,\quad 9,11, 12, 13,14,15\}$ & & 0.019 & \multicolumn{4}{c}{} \\
\bottomrule
\end{tabular}
\end{center}
\label{table:uscrime}
\end{table}

\section{Concluding remarks}\label{sec:cr}
Bayesian model selection for linear regression models with Gaussian errors
has been popular area of study for some time. There is also a substantial
literature devoted to studying the extension of Stein-type shrinkage estimators
from models with Gaussian errors to those with general SS 
errors. In particular, it has long been observed that certain shrinkage
estimators which improve over the least squares (LS) estimator for
Gaussian models also improve over the LS estimator simultaneously for all
SS  error models 
(See for example, \cite{Cellier-et-al-1989}).
\cite{Maruyama-2003b} and \cite{Maru-Straw-2005} found, in addition, that
certain separable priors (in the sense described in Section \ref{sec:prior})
leads to generalized Bayes shrinkage estimators that do not depend on
the form of the underlying SS  distribution and that
also simultaneously improve on the LS estimator, sometimes dramatically so.
The original aim of this research was to see if similar separable priors
could be found that have this distributional robustness property in the
variable selection problem (and, of course, also to perform well).

The generalized Bayes priors developed in sections \ref{sec:prior} and
\ref{sec:marginal} turned out to satisfy our requirements and also to be closely 
related to other so called $g$-priors (or mixtures of $g$-priors) in the
literature (See e.g.~\cite{Liang-etal-2008,Guo-Speckman-2009}).
Our original development required that the null model be excluded from 
the class of possible models, however we observed (see Remark \ref{rem:null})
that a slightly altered version which treats the intercept term the same as
all other regression coefficients (as opposed to giving it a uniform prior),
allows the null model to be included as well.
This alternative development gives Bayes factor which depend on the
$\check{R}^2$ relative to the null model where all coefficients including $\alpha$
is $0$. This dependence on $\check{R}^2$ is related to the model selection procedures of 
\cite{Guo-Speckman-2009,Celeux-Anbari-Marin-Robert-2011}.

Our prior distribution on the regression parameters is sub-harmonic for
each (non-null) model, and as such, leads to admissible estimators under
quadratic loss in each of the non-null models when the  variance is
known (see \cite{Maruyama-Takemura-2008}) regardless of the SS 
error distribution. Also if $q \geq 3$ and $\nu=2$ the generalized
Bayes estimator for each submodel is minimax and dominates the James-Stein
estimator (See \cite{Maruyama-2003b}). Hence although we require
$ 0<\nu<1$, we nevertheless expect good performance of the corresponding
Bayes estimators.

The expression of our Bayes factors, e.g.~\eqref{eq:main-thm-BF},
are relatively simple involving the ratio of two $1$-dimensional integrals.
To further simplify calculations we investigated Laplace approximations
to our Bayes factors, and more generally, to a collection of Bayes factors arising from
mixtures of $g$-priors that have recently appealed 
(See \cite{Liang-etal-2008, Guo-Speckman-2009}). 
We show in Section 4 that in each case the Bayes factor can be approximated
as the Bayes factor for the Gaussian model based on BIC times a simple
rational function depending $q$, $\nu$ and the $R^2$ of the models.

Using these Laplace approximations we are able to establish model selection
consistency of our robust  procedure for the entire class of SS  distributions
and to extend the model consistency results of several earlier papers
for the Gaussian case to the entire class of SS  distributions.

A small simulation study and an analysis of the Hald data 
(See \cite{Casella-Moreno-2006}) 
and the US Crime data (See \cite{Raftery-Madigan-Hoeting-1997}) indicates that our method 
performs well. It gives results consistent with the results of the cited 
papers for the real data sets and performs comparably
and sometimes better than several of the mixture of $g$-prior methods.

