\section{The LArTPC technique}
Liquid Argon Time Projection Chamber (LArTPC) technology~\cite{rubbia} provides mm-scale position resolution, calorimetry, three dimensional imaging, and efficient particle identification for neutrino detection and reconstruction. In a LArTPC, the interaction products exiting the neutrino vertex ionize argon atoms as they traverse the detector. As argon is a noble element, this ionization is free to drift across the medium. An electric field is imposed in the detector and the ionization tracks are drifted towards readout wire planes that are oriented with respect to one another at an angle. The ionization induces a current on the induction plane and is collected by the collection plane. Combining the information from both planes and detecting the events in time allows a three dimensional image of the event with calorimetric information to be obtained. Argon, with a relatively large density, short radiation length, and high scintillation yield, is ideal for neutrino detection and event containment. A schematic depicting the LArTPC concept for neutrino detection can be seen in Figure~\ref{lartpc}.

\begin{figure}[h!]
\centering
\begin{tabular}{c}
\epsfig{file=./lartpc.pdf,width=4.9in,angle=0}
\end{tabular}
\vspace{-2.5cm}
\caption{The LArTPC concept for neutrino detection.}
\label{lartpc}
\end{figure}

\section{ArgoNeuT}
ArgoNeuT (NSF/DOE) sat just upstream of the on-axis MINOS~\cite{MINOS} near detector in the NuMI beamline, about 1~km from the target station and 100~m underground at Fermilab (see Figure~\ref{detectorpic}). The experiment features a 170~L active volume TPC with a fully contained recirculation and purification system. The main detector characteristics are summarized in Table~\ref{vitals}. ArgoNeuT collected about 1.35E20 protons on target (POT) with the ``low energy" NuMI configuration (0.1E20 POT in neutrino mode and 1.25E20 POT in anti-neutrino mode) from September 2009 to February 2010. Run operations were largely stable and shift-free over the $>$5 month span. Along with a Charged Current Quasi-Elastic (CCQE) cross section measurement (one of many potential cross section measurements) ArgoNeuT will be demonstrating the dE/dx and topological-based particle identification capabilities of LArTPCs with a special emphasis on electron/photon differentiation, a vital step toward the realization of a future background-free electron-neutrino appearance search. ArgoNeuT is also developing reconstruction techniques specific to LArTPCs in collaboration with MicroBooNE~\cite{microboone} and LBNE~\cite{lbne} through the ``LArSoft" collaboration and performing research and design for future LArTPCs. Indeed, ArgoNeuT is one of the key steps in the phased US LArTPC program moving toward construction of a $>$10~kiloton detector to be used for long-baseline and atmospheric neutrino oscillation physics, a proton decay search, supernova burst and diffuse neutrino detection, and more~\cite{modular,landd,glacier,flare,icanoe2,Conrad:2010mh,Cocco:2004ac,Rubbia:2004yq,Arneodo:2001tx}.

\begin{table}[h]
\centering
\begin{tabular}{|c|c|}
	\hline
Cryostat volume  &500 liters  \\ \hline
TPC volume  &175 liters  \\ \hline
\# Electronic channels  &480 (240/plane)  \\ \hline
Wire spacing  &4 mm  \\ \hline
Electronics style (temperature)  &JFET (293 K)  \\ \hline
Max. drift length (time)  &50 cm (330 $\mu$s)  \\ \hline
Light collection?  &No  \\
\hline
\end{tabular}
\caption{The ArgoNeuT detector's main specifications.}
\label{vitals}
\end{table}
\vspace{-.5cm}
\begin{figure}[h!]
\centering
\begin{tabular}{c c}
\hspace{-1.8cm}
\epsfig{file=./detectorpic.pdf,width=4.4in,angle=0}
&
\hspace{-4.cm}
\epsfig{file=./fnal.pdf,width=4.4in,angle=0}
\end{tabular}
\vspace{0.0in}
\caption{(Left) ArgoNeuT sitting just in front of the MINOS near detector in the NuMI beamline. (Right) The location of ArgoNeuT in the MINOS hall, relative to the Fermilab accelerator complex.}
\label{detectorpic}
\end{figure}

\section{The CCQE cross section and vertex activity}
The Charged Current Quasi-Elastic (CCQE) interaction is the most consequential neutrino interaction channel for accelerator-based neutrino oscillation experiments. The $\nu_{\mu}$ ($\overline{\nu}_{\mu}$) CCQE interaction $\nu_{\mu}n\longrightarrow\mu^{-}p$ ($\overline{\nu}_{\mu}p\longrightarrow\mu^{+}n$) is used to measure the muon-neutrino flux, constraining the expected electron-neutrino flux at the near and far detectors. Long baseline electron neutrino appearance searches (MINOS~\cite{MINOS}, T2K~\cite{t2k}, NOvA~\cite{nova}, LBNE~\cite{lbne}, etc.) all employ a near and far detector in order to minimize neutrino flux uncertainties. However, due to the non-zero near-far detector systematics (i.e. differences between the near and far detector) inherent in all long baseline experiments, cross section assumptions need to be made\footnote{Such systematics come from the differences in the neutrino energy spectra seen at the near and far detector. Also, systematics can come from inherent differences in detector parameters (detector medium, size, acceptance, etc.) between near and far.}. Knowledge of the neutrino cross section is required even in the case of a near-far configuration.  The simple two particle event topology and dominant cross section in the most relevant energy region makes the CCQE interaction the golden channel for long baseline neutrino oscillation experiments. Despite being the golden channel, however, the CCQE cross section is known with only 20-30\% precision over most energies. Furthermore, only one (preliminary) measurement of the CCQE cross section on argon has ever been taken (at $<E_{\nu}>$=28~GeV)~\cite{icaruswanf,wanfprd}.

The neutrino energy can be reconstructed with the lepton energy and angle with respect to beam axis alone in the case of a CCQE interaction. LArTPCs have the added advantage of the ability to reconstruct the proton(s) that come out of the interaction and, with knowledge of the lepton kinematics, fully reconstruct the neutrino energy. The additional kinematic constraints from the reconstructed proton(s) aid in neutrino energy resolution and background suppression. The muon and proton(s) are identified unambiguously in a $\nu_{\mu}$ CCQE event with a combination of dE/dx and range measurements in ArgoNeuT. The downstream magnetized MINOS near detector is used to sign-select and, with ArgoNeuT, fully reconstruct the muon while the proton(s) (contained $>$50\% of the time) is reconstructed with ArgoNeuT alone.

CCQE cross section measurements are subject to large systematics coming from final state interaction (FSI) modeling, flux uncertainties, background neutrino charged current pion interaction cross section uncertainties, and more. Charged current pion events ($\nu_{\mu}p\longrightarrow\mu^{-}p\pi^{+}$, $\nu_{\mu}n\longrightarrow\mu^{-}p\pi^{0}$, or $\nu_{\mu}n\longrightarrow\mu^{-}n\pi^{+}$) can appear ``CCQE-like" if the pion is absorbed (e.g. $\pi N \rightarrow NN$) during an interaction within the nucleus. More CCQE backgrounds can come from interactions involving final state neutral particles (neutrons, gammas) that escape detection and/or neutral current events with a final state pion incorrectly identified as a muon that renders the event ``CCQE-like". In true CCQE events the nucleon is also subject to various FSI processes including pion production, charge exchange, and elastic and inelastic scattering within the parent nucleus. Flux-integrated differential cross sections with muon kinematics and cross sections based on post-FSI observable particles only (rather than pre-FSI, which introduces a model dependence in an attempt to track back to the neutrino interaction itself) are now starting to be reported alongside absolute total pre-FSI cross sections~\cite{miniccqe}.

Recent MiniBooNE~\cite{miniccqe} and NOMAD~\cite{nomad} CCQE cross section measurements disagree under the relativistic Fermi gas model (Ref.~\cite{smithmoniz} is commonly used) by up to 30\% or more, despite using the same nuclear target ($^{12}$C). Surprisingly, the MiniBooNE CCQE cross section measurement (the larger of the two reported cross sections) surpasses the free neutron cross section. The discrepancy may be due to a CCQE multinucleon channel in which two correlated same-flavor nucleons are ejected~\cite{martini} and the experiments' varying levels of blindness to low energy protons coming out of the vertex or ``vertex activity". The channel may be the result of both short range correlations and 2 body currents between nucleons. Protons are not considered in the MiniBooNE CCQE measurement and protons with momentum less than $\sim$300 MeV/c are not reconstructed in the NOMAD analysis. MiniBooNE's complete insensitivity to the proton(s) is advantageous in the sense that the extracted cross section measurements do not rely as much on proton FSI modeling uncertainties~\cite{miniccqe}. However, the experiment is unable to measure a possible multinucleon excitation channel.

Most neutrino detectors have a difficult time resolving vertex activity involving short tracks below Cherenkov threshold. ArgoNeuT, with its mm-scale resolution, three dimensional imaging, and few-10s of MeV energy threshold, will analyze the kinematics of vertex activity in CCQE interactions in an attempt to understand the multinucleon CCQE channel and final state interactions in general. Figure~\ref{ccqe} shows a CCQE interaction in ArgoNeuT with clear evidence of vertex activity through a set of $<$100~MeV protons. With the help of the downstream magnetized MINOS detector, ArgoNeuT will be able to compare neutrino and anti-neutrino CCQE events in terms of vertex activity. A multinucleon channel involving a 2 proton (neutron) final state before final-state-interactions in CCQE neutrino (anti-neutrino) events will be searched for. Comparing neutrino and anti-neutrino CCQE-like events may provide some sensitivity to this channel. Event characteristics such as energy deposited in a sphere of cm-scale radius around the vertex, number of protons, and proton kinematics in general and the differences between neutrino and anti-neutrino events in terms of these variables will be studied. Furthermore, the possibility of separating neutrino and anti-neutrino CCQE events with LArTPCs using vertex activity alone will be investigated.

\begin{figure}[h!]
\centering
\begin{tabular}{c}
\hspace{-.3cm}
\epsfig{file=./ccqe.pdf,width=5.0in,angle=0}
\end{tabular}
\vspace{-2.5cm}
\caption{A zoomed-in neutrino CCQE candidate with clear evidence of multiple, highly ionizing protons and a muon (the topmost track) emanating from the vertex. The colors are indicative of the amount of energy deposited (blue$<$yellow$<$red). For a sense of scale, note that the wire spacing is 4~mm.}
\label{ccqe}
\end{figure}

\section{Reconstructing the muon and finding the neutrino vertex}
The initial goals of ArgoNeuT's automated neutrino reconstruction software are to reconstruct the kinematics of the outgoing muon and find the neutrino interaction vertex. These steps are the key to beginning a CCQE analysis with ArgoNeuT. The muon reconstruction chain proceeds as follows:
\begin{itemize}
\item Event filter, to remove empty beam spills from the sample.
\item Deconvolution of the TPC and electronics response from the raw TPC-wire data, to produce standardized unipolar signal shapes.
\item Identification of hits above a threshold having a position, width, and projected ionization associated with each particle track.
\item DBSCAN (Density-Based Spatial Clustering of Applications with Noise)-based~\cite{dbscan} cluster finding, to group proximal hits together to form clusters and to identify noise hits.
\item Hough-transform-based~\cite{hough} cluster finding and track fitting, to find and fit line-like clusters/tracks.
\item Three dimensional track matching with the wire signal and time information from both wire planes.
\item Hit analysis and calorimetry, to characterize the hits that have been associated with long line-like tracks.
\end{itemize}

The vertex finding process takes hits that have been associated with DBSCAN clusters as an input and employs a Harris/Stephens-based~\cite{harris} algorithm to find the vertex and endpoints in the event. The algorithm considers the hits from a geometric point of view only with an implementation based on~Ref.~\cite{morgan}. These vertices/endpoints, after being given a ``strength", are matched with the endpoints of lines found with the Hough-transform-based line-finding/fitting algorithm to create ``strong vertices". An example of the reconstruction and endpoint/vertex finding algorithms at work can be seen in Figure~\ref{vertex}.

\begin{figure}[h!]
\centering
\begin{tabular}{c}
\hspace{-.3cm}
\epsfig{file=./vertex.pdf,width=5.5in,angle=0}
\end{tabular}
\vspace{-2.5cm}
\caption{A few examples of the ArgoNeuT automated reconstruction software. (Top Left) A probable neutral current $\pi^{0}$ event ($\nu N\longrightarrow\nu N\pi^{0}$, $\pi^{0}\longrightarrow\gamma\gamma$) involving two energetic gamma-induced electromagnetic showers. The hit finding and density-based clustering algorithms show the hits and clusters associated with the event. (Top Right) A probable charged current resonant pion event ($\nu_{\mu}N\longrightarrow\mu^{-}N\pi^{+}$) reconstructed in three dimensions. (Bottom) The same event as seen on the induction plane with overlaid reconstructed line-like tracks and endpoints/vertices.}
\label{vertex}
\end{figure}

\section{Conclusion}
ArgoNeuT recently completed its $>$5 month physics run in the NuMI beamline at Fermilab. The first few thousand low energy neutrino and anti-neutrino events ever detected with a LArTPC are currently being analyzed with a fully automated software capable of reconstructing the muon and finding the neutrino interaction vertex. Among other channels, ArgoNeuT is investigating the CCQE cross section with a focus on the vertex activity associated with such events. With the unique ability to sign-select muons via the downstream MINOS near detector, ArgoNeuT will perform a direct neutrino/anti-neutrino comparison of low energy protons emanating from the interaction site. This measurement may be sensitive to a multinucleon excitation channel involving two same-flavor nucleons (protons, in the case of neutrino CCQE). This study will also seek to illuminate the influence of nuclear effects on the absolute and differential CCQE cross sections and on the experimental definition of CCQE itself.

\section*{References}
\bibliographystyle{iopart-num}

\section{Introduction}

The \texttt{iopart-num} Bib\TeX{} style is intended for use in
preparing manuscripts for Institute of Physics Publishing journals,
including Journal of Physics.  It provides numeric citation with
Harvard-like formatting, based upon the specification in ``How to
prepare and submit an article for publication in an IOP journal using
\LaTeXe'' by Graham Douglas (2005).

The \texttt{iopart-num} package is available on the Comprehensive
\TeX{} Archive Network (CTAN) as \texttt{/biblio/bibtex/contrib/iopart-num}.

\section{General instructions}

To use the \texttt{iopart-num} style, include the command
\verb+\bibliographystyle{iopart-num}+ in the
document preamble.  The reference section is then inserted into the
document with the command \verb+


\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.