Prior-based Bayesian information criterion

ABSTRACT

We present a new approach to model selection and Bayes factor determination, based on Laplace expansions (as in BIC), which we call Prior-based Bayes Information Criterion (PBIC). In this approach, the Laplace expansion is only done with the likelihood function, and then a suitable prior distribution is chosen to allow exact computation of the (approximate) marginal likelihood arising from the Laplace approximation and the prior. The result is a closed-form expression similar to BIC, but now involves a term arising from the prior distribution (which BIC ignores) and also incorporates the idea that different parameters can have different effective sample sizes (whereas BIC only allows one overall sample size n). We also consider a modification of PBIC which is more favourable to complex models.

KEYWORDS:

• Bayes factors
• model selection
• Cauchy priors
• consistency
• effective sample size
• Fisher information
• Laplace expansions
• robust priors

1. Background

1.1. The original BIC (Schwarz, 1978)

Suppose that we observe ${\mathit{X}}_i=(X_{i1},\dots ,X_{ip})\sim g({\mathbf{x}}_i\mid \mathit{\theta })$ for $i=1,\dots ,n$. Here $\mathit{\theta }=({\theta }_1,\dots ,{\theta }_p)$ is a unknown vector and, in Schwartz's derivation of BIC, $g(\mathbf{x}\mid \mathit{\theta })$ is an exponential family. Then the log-likelihood function is $l\left(\mathit{\theta }\right)=\log f(\mathbf{x}\mid \mathit{\theta })=\log \left(\prod _{i=1}^{n}g({\mathbf{x}}_i\mid \mathit{\theta })\right),$ where $\mathbf{x}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n)$. The goal of Schwarz (1978) is to find a simple approximation to the marginal density $m\left(\mathbf{x}\right)=\int f(\mathbf{x}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta },$ where $\pi \left(\mathit{\theta }\right)$ is a prior density for the unknown $\mathit{\theta }$, and use the approximation for model comparison.

Result 1.1

Stone, 1979

Let $\hat{\mathit{\theta }}$ be the MLE of $\mathit{\theta }$. Then, under reasonable conditions and as $n\to \mathrm{\infty }$, $\text{BIC}\equiv -2l\left(\hat{\mathit{\theta }}\right)+p\log n=-2\log m\left(\mathbf{x}\right)+c+o\left(1\right),$ where c is a constant.

Schwartz then suggested comparing two models $M_1$ and $M_2$, using $\mathrm{\Delta }\text{BIC}={\text{BIC}}_2-{\text{BIC}}_1,$

preferring $M_2$ ($M_1$) as this is negative (positive). Clearly this is equivalent to basing the model comparison on the Bayes factor (odds) of $M_2$ to $M_1$, with the approximation (1) $\begin{array}{rl}{B}_{21}& \equiv \frac{m_2\left(\mathbf{x}\right)}{m_1\left(\mathbf{x}\right)}=\frac{\exp \left(-\frac{1}{2}{\text{BIC}}_2\right)}{\exp \left(-\frac{1}{2}{\text{BIC}}_1\right)}\exp \left(\frac{1}{2}(c_2-c_1)\right)\\ & \times (1+o(1\left)\right)\approx \frac{\exp \left(-\frac{1}{2}{\text{BIC}}_2\right)}{\exp \left(-\frac{1}{2}{\text{BIC}}_1\right)}.\end{array}$(1)

1.2. Problems with general use of BIC

BIC is an excellent tool for the class of problems for which it was developed. Unfortunately, it is today used ubiquitously, for completely different classes of problems. We here outline some of the issues with using BIC inappropriately.

Problem 1. The term $\exp \mathbf{\left(}\frac{\mathbf{1}}{\mathbf{2}}\mathbf{\right(}{\mathbf{c}}_{\mathbf{2}}\mathbf{-}{\mathbf{c}}_{\mathbf{1}}\mathbf{\left)}\mathbf{\right)}$ in (1(1) $\begin{array}{rl}{B}_{21}& \equiv \frac{m_2\left(\mathbf{x}\right)}{m_1\left(\mathbf{x}\right)}=\frac{\exp \left(-\frac{1}{2}{\text{BIC}}_2\right)}{\exp \left(-\frac{1}{2}{\text{BIC}}_1\right)}\exp \left(\frac{1}{2}(c_2-c_1)\right)\\ & \times (1+o(1\left)\right)\approx \frac{\exp \left(-\frac{1}{2}{\text{BIC}}_2\right)}{\exp \left(-\frac{1}{2}{\text{BIC}}_1\right)}.\end{array}$(1) ) is ignored by BIC.

This could have been a serious problem even with proper use of BIC, except that there happens to be pseudo-prior distributions that yield BIC itself (Raftery, 1999), i.e. for which the term $\exp \left(\frac{1}{2}\right(c_2-c_1\left)\right)=1$. These pseudo-priors are not real priors, in that they are centred at the mle's of each model, which is a problematical double use of the data. Nevertheless it is comforting that there is at least some type of prior distribution that yields BIC exactly.

Problem 2. What is $\mathbf{n}$?

1. A common mistake in specifying n: Note that, in Schwartz's setup, there are n vector observations of dimension p, so that there are a total of np real observations. It is common to mistakenly use $n^{\ast }=np$ as the sample size in BIC, rather than the correct n.

2. Different parameters can have different n.

Example 1.2

Group means

For $i=1,\dots ,p$ and $l=1,\dots ,r$, suppose we observe $X_{il}={\mu }_i+{\epsilon }_{il},$ where ${\epsilon }_{il}\sim N(0,{\sigma }^2)$. If ${\sigma }^2$ were known, this would be exactly the setup of Schwartz, and the sample size for $\mathit{\mu }=({\mu }_1,\dots ,{\mu }_p)$ would be r. In effect, each ${\mu }_i$ has a sample size of r associated with it. But, if ${\sigma }^2$ is unknown, the parameter is $\mathit{\theta }=({\mu }_1,\dots ,{\mu }_p,{\sigma }^2)$ and it is not reasonable to also associate the sample size of r to ${\sigma }^2$, in that we know there are $p(r-1)$ degrees of freedom associated with the mle of ${\sigma }^2$.

An alternative argument is to note that the observed information matrix $\hat{\mathit{I}}=\left({\hat{I}}_{jk}\right)$, with $(j,k)$ entry ${\hat{I}}_{jk}=-\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\theta }_j\mathrm{\partial }{\theta }_k}\log f(\mathbf{x}\mid \mathit{\theta }){|}_{\mathit{\theta }=\hat{\mathit{\theta }}}$ is given by $\hat{\mathit{I}}=\left(\begin{array}{cc}\frac{r}{{\hat{\sigma }}^2}{I}_{p\times p}& 0\\ 0& \frac{pr}{2{\hat{\sigma }}^4}\end{array}\right),$ where ${\hat{\sigma }}^2=\left(1/pr\right)\sum _{i=1}^p\sum _{l=1}^r(X_{il}-{\overline{X}}_i{)}^2$. The information matrix suggests that the effective sample size for each ${\mu }_i$ is r, while the effective sample size for ${\sigma }^2$ is pr. Whether we use $p(r-1)$ or pr for the sample size associated with ${\sigma }^2$ will not typically make much difference, whereas the difference with using r, instead, will be quite large.

1. Different observations can have different observed information content.

Example 1.3

Suppose each independent observation, $X_i,i=1,\dots ,n$, has probability 1/2 of arising from the $N(\theta ,1)$ distribution and probability 1/2 of arising from the $N(\theta ,1000)$ distribution. Clearly half the observations are essentially worthless, and the ‘effective sample size’ is n/2.

Example 1.4

Findley's BIC counterexample

One of the famous counter examples against inappropriate use of BIC is in Findley (1991). Suppose the observations are (2) $\begin{array}{rl}{X}_i& =\frac{1}{\sqrt{i}}\cdot \theta +{\epsilon }_i,\text{where }{\epsilon }_i\sim N(0,1),\\ i& =1,\dots ,n,\end{array}$(2) and we are comparing the models $H_0:\theta =0$ and $H_1:\theta \ne 0$. It turns out that the mle for θ is consistent under $H_1$ (a necessary condition to apply BIC), but that BIC is inconsistent if $0<|\theta |<1$, in that BIC will then declare $H_0$ to be the true model as $n\to \mathrm{\infty }$. The problem here is that, even though the information about θ goes to ∞ as n grows, it grows much more slowly than n (actually, the information grows at roughly $\log n$ rate), and BIC erroneously assigns the rate to be n.

Problem 3. What is p?

Just as n is often not clearly defined for use in BIC, the parameter dimension p is often not clearly defined (see also Pauler, 1998).

Example 1.5

Random effect group means

Consider hierarchical or random effect versions of the group means problem, where it is assumed that ${\mu }_i\sim N(\xi ,{\tau }^2),$ with ξ and ${\tau }^2$ being unknown. The number of parameters might appear to be p+3 (the means, along with ${\sigma }^2$, ξ and ${\tau }^2$), but one could, alternatively, integrate out $\mathit{\mu }=({\mu }_1,\dots ,{\mu }_p)$ (since it has a known distribution) obtaining $\begin{array}{rl}f(\mathbf{x}\mid {\sigma }^2,\xi ,{\tau }^2)& =\int f(\mathbf{x}\mid \mathit{\mu },\xi ,{\sigma }^2)\pi (\mathit{\mu }\mid \xi ,{\tau }^2)\mathrm{d}\mathit{\mu }\\ & \propto \frac{1}{{\sigma }^{-p(r-1)}}\exp \left\{\frac{{\hat{\sigma }}^2}{2{\sigma }^2}\right\}\prod _{i=1}^p\\ & \times \exp \left\{-\frac{({\overline{x}}_i-\xi {)}^2}{2(\frac{{\sigma }^2}{r}+{\tau }^2)}\right\}.\end{array}$ The marginal likelihood will be the integral of this, with respect to a prior $\pi ({\sigma }^2,\xi ,{\tau }^2)$, so that, if one is really viewing BIC as an approximation to the marginal likelihood, it would be correct to set p=3.

Problem 4. What if p grows with n?

BIC is based on an asymptotic argument with p fixed and n growing, but often p is growing with n; BIC then does not apply. If one were to erroneously apply BIC in such a situation, one could end up with inconsistency, as shone by Stone (1979) for the group means example, with known variance ${\sigma }^2=1$ for simplicity. Indeed, in comparing models $H_0:\mathit{\mu }=\mathbf{0}$ and $H_1:\mathit{\mu }\ne \mathbf{0}$ for the group means problem with r=2, $\mathrm{\Delta }\text{BIC}={\text{BIC}}_1-{\text{BIC}}_0=-2\sum _{i=1}^p{\overline{x}}_i^2+p\log 2,$ which, under $H_0$, behaves like $p(\log 2-1)\to -\mathrm{\infty }$ as p grows, thus incorrectly selecting model $H_1$.

1.3. Variants of BIC

Noting the limitations of BIC, researchers have proposed a host of generalisations, many of which have performed better than BIC under specific scenarios. Many of these methods arise from the variations in retaining the number of terms in the Laplace approximation of the Bayes Factor (Kass & Raftery, 1995). One variant – called the HBIC – (Haughton, 1988) retains the third term in the Laplace approximation of the Bayes Factor. A simulation study by Haughton, Oud, and Jansen (1997) shows that HBIC performs better in model selection for structural equation models than does the usual BIC. Following HBIC, Bollen, Ray, Zavisca, and Harden (2012) developed a similar criterion, called the Information matrix-based Bayesian Information Criterion (IBIC), which retains more terms in the Bayes Factor approximation and outperforms BIC and HBIC in many scenarios. Bollen et al. (2012) also proposed another criterion, named the scaled unit information prior (SPBIC), which generalises the interpretation of the unit information prior in the context of BIC. For approximation of Bayes factors as the model dimension grows, Berger, Ghosh, and Mukhopadhyay (2003) proposed another approximation, named GBIC. Following Berger et al. (2003), a generalisation of BIC for the general exponential family was proposed by Chakrabarti and Ghosh (2006), and a new BIC for change point analysis was proposed by Shen and Ghosh (2011). Some other extensions of BIC include techniques for comparing graphical models (Foygel & Drton, 2010), singular models (Drton & Plummer, 2017), and sparse models (Zak-Szatkowska & Bogdan, 2011).

1.4. Overview of the paper

Section 2 presents a proposal to generalise BIC, in order to overcome the problems mentioned above. It is based on use of a specific (robust) prior distribution in the computation of the approximate marginal likelihood of a model. Section 3 discusses a critical aspect of the definition of PBIC, namely the need to determine the ‘effective sample size’ corresponding to each parameter in a model. Section 4 presents an alternative called PBIC*. It employs an empirical Bayes prior in computation of the marginal likelihood approximation, resulting in answers more favourable to complex models. Section 5 illustrates the use of PBIC and PBIC* in the normal linear model; it is of interest that PBIC and PBIC* correspond to exact marginal likelihoods here. Illustrations in the section are simple linear regression, testing the equality of normal means with known unequal variances, Findley's counterexample, and the group means problem, where consistency results for PBIC and PBIC* are established as $p\to \mathrm{\infty }$.

2. The PBIC solution

We propose a solution to these problems that depends only on software that can compute mle's and observed information matrices. The basis of the solution is a modified Laplace approximation to $m\left(\mathbf{x}\right)$ for reasonable default priors.

2.1. Two important preliminaries

One should analytically integrate out any parameter that has a distribution given other parameters, if it is possible to do so. For example, in the hierarchical group means example, base the analysis on the marginal likelihood $f(\mathbf{x}\mid {\sigma }^2,\xi ,{\tau }^2)$, rather than the full likelihood.

We will be utilising the Laplace approximation, which is most accurate (Kass & Vaidyanathan, 1992; Tierney, Kass, & Kadane, 1989) if the parameter space is transformed to be all of ${\mathrm{\Re }}^p$. Transformation to ${\mathrm{\Re }}^p$ will also be necessary for the subsequent step of the analysis. As an illustration, in the (non-multilevel) group means example, transform to $\nu =\log {\sigma }^2$. Then $\mathit{\theta }=({\mu }_1,\dots ,{\mu }_p,\nu )\in {\mathrm{\Re }}^{(p+1)}$. Note that one then works with the transformed mle $\log {\hat{\sigma }}^2$ and the transformed observed information matrix ${\hat{\mathit{I}}}^{\ast }(\mathit{\mu },\nu )=\left(\begin{array}{cc}\frac{r}{{\hat{\sigma }}^2}{I}_{p\times p}& 0\\ 0& \frac{pr}{2}\end{array}\right).$ In the multilevel group means example, both ${\sigma }^2$ and ${\tau }^2$ would need to be transformed in this fashion.

2.2. PBIC and PBIC* definitions

Suppose $\mathit{\theta }=({\mathit{\theta }}_{\left(1\right)},{\mathit{\theta }}_{\left(2\right)})$, where ${\mathit{\theta }}_{\left(2\right)}$ denotes the parameters that are common to all models under consideration (e.g. an intercept in linear regression). Changing notation, let p denote the dimension of ${\mathit{\theta }}_{\left(1\right)}$ and q denote the dimension of ${\mathit{\theta }}_{\left(2\right)}$; note that p will typically vary from model to model, while q is fixed. Partition the observed information matrix for a model accordingly, as (3) $\begin{array}{rl}\hat{\mathit{I}}& =\left(\begin{array}{cc}{\hat{\mathit{I}}}_{11}& {\hat{\mathit{I}}}_{12}\\ {\hat{\mathit{I}}}_{21}& {\hat{\mathit{I}}}_{22}\end{array}\right),\text{and define }\\ {\mathbf{\Sigma }}^{-1}& ={\hat{\mathit{I}}}_{11}-{\hat{\mathit{I}}}_{12}{\hat{\mathit{I}}}_{22}^{-1}{\hat{\mathit{I}}}_{12}^t.\end{array}$(3) (If there are no common parameters to all models, then $\mathbf{\Sigma }={\hat{\mathit{I}}}^{-1}$.) Change variables to $\mathit{\xi }=\mathit{O}{\mathit{\theta }}_{\left(1\right)}$, where $\mathit{O}$ is an orthogonal matrix such that $\mathbf{\Sigma }={\mathit{O}}^t\mathit{D}\mathit{O}$, with $\mathit{D}=\text{diag}\left\{d_i\right\}$ for $i=1,\dots ,p$, and define $\hat{\mathit{\xi }}=\mathit{O}{\hat{\mathit{\theta }}}_{\left(1\right)}$ (the transformed mle). The choice of $\mathit{O}$ does not affect the definition below. For each transformed parameter ${\xi }_j$, let $n_j^e$ be the effective sample size corresponding to that parameter. This is the most difficult aspect of the construction, but equals the intuitive choices of parameter sample size discussed in the earlier examples; formal definitions will be presented in Section 3. Then PBIC is defined as (4) $\begin{array}{rl}\text{PBIC}& \equiv -2l\left(\hat{\mathit{\theta }}\right)+\log |{\hat{\mathit{I}}}_{22}|+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i},\end{array}$(4) where $v_i={\hat{\xi }}_i^2/\left[d_i\right(1+n_i^e\left)\right]$. For a certain natural prior distribution, PBIC will be shown to be accurate, as an approximation to $-2\log m\left(\mathit{x}\right)$, up to a $o\left(1\right)$ term as $n\to \mathrm{\infty }$ (for fixed dimension p). Note that, if there are no common parameters to all models, then (5) $\begin{array}{rl}\text{PBIC}& =-2l\left(\hat{\mathit{\theta }}\right)+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i}.\end{array}$(5) In the classic case considered by Schwartz, all $n_i^e$ would equal a common n, and the first two terms in this expression are then BIC (up to a o(1) term); the ‘constant’ ignored by BIC is the final term in (5(5) $\begin{array}{rl}\text{PBIC}& =-2l\left(\hat{\mathit{\theta }}\right)+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i}.\end{array}$(5) ).

To summarise results in one place, here is an alternative version of the approximation, one which is more favourable to complex models; its development is given in Section 4: (6) $\begin{array}{rl}{\text{PBIC}}^{\ast }& \equiv -2l\left(\hat{\mathit{\theta }}\right)+\log |{\hat{\mathit{I}}}_{22}|+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}}.\end{array}$(6) Note that, if dealing with only normal mean parameters, PBIC and PBIC* are exact as an approximation to $-2\log m\left(\mathit{x}\right)$, as discussed below. This would mean, for instance, that when dealing with $p\to \mathrm{\infty }$, there would be no need to worry about the accuracy of the approximations.

Here are the steps in the derivation of PBIC.

2.2.1. Laplace approximation

By a Taylor's series expansion of $l\left(\mathit{\theta }\right)$ about the mle $\hat{\mathit{\theta }}$, with ∇ denoting the gradient and $\hat{\mathit{I}}$ being the observed information matrix as defined earlier, (7) $\begin{array}{rl}m\left(\mathbf{x}\right)& =\int f(\mathbf{x}\mid \mathit{\theta })\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & =\int e^{l\left(\mathit{\theta }\right)}\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & =\int \exp \left[l\left(\hat{\mathit{\theta }}\right)+(\mathit{\theta }-\hat{\mathit{\theta }}{)}^t\mathrm{\nabla }l(\hat{\mathit{\theta }})\phantom{\frac{1}{2}}\right.\\ & -\left.\frac{1}{2}(\mathit{\theta }-\hat{\mathit{\theta }}{)}^t\hat{\mathit{I}}(\mathit{\theta }-\hat{\mathit{\theta }})\right]\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }(1+o(1\left)\right)\\ & =e^{l\left(\hat{\mathit{\theta }}\right)}\int e^{-\left(1/2\right)(\mathit{\theta }-\hat{\mathit{\theta }}{)}^t\hat{\mathit{I}}(\mathit{\theta }-\widehat{\mathit{\theta })}}\pi \left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }(1+o(1\left)\right),\end{array}$(7) where $o\left(1\right)$ denotes a term that goes to zero as the sample size n grows. Technical conditions for the validity of this Laplace approximation can be found in, e.g. Tierney et al. (1989), Kass and Vaidyanathan (1992); the key assumption needed is that $\hat{\mathit{\theta }}$ occurs on the interior of the parameter space, so that $\mathrm{\nabla }l\left(\hat{\mathit{\theta }}\right)=0$. (If not true, the analysis must proceed as in Dudley & Haughton, 1997; Haughton, 1991, 1993). Also, the presence of $o\left(1\right)$ assumes that p is fixed as n grows. We will nevertheless use this approximation, even as p grows with n, relying on the considerable evidence that the Laplace approximation is quite generally accurate.

Note that we do not use the more common version of the Laplace expansion which involves $\pi \left(\mathit{\theta }\right)$ in the Taylor's expansion because we will be choosing $\pi \left(\mathit{\theta }\right)$ so that the integral in this expression can be evaluated in closed form. In particular, this means that, if we are dealing with the situation where $\mathit{\theta }$ is the mean parameter of a normal model, then the computations herein will be entirely closed form, with no approximation being involved (and no need to then worry about p growing with n).

2.2.2. Choosing a good prior $\pi \left(\mathit{\theta }\right)$

Assume that the transformations in Section 2.1 have been made.

Step 1. Recall that $\mathit{\theta }=({\mathit{\theta }}_{\left(1\right)},{\mathit{\theta }}_{\left(2\right)})$, where ${\mathit{\theta }}_{\left(2\right)}$ denotes the common parameters to all models. We will utilise a prior distribution $\pi \left(\mathit{\theta }\right)=\left(2\pi {)}^{-q}\pi \right({\mathit{\theta }}_{\left(1\right)}),$ where $\pi \left({\mathit{\theta }}_{\left(1\right)}\right)$ is defined later. The key point is that, since ${\mathit{\theta }}_{\left(2\right)}$ is common to all models, it can be assigned a constant prior density (see, e.g. Bayarri, Berger, Forte, & García-Donato, 2012; Berger, Pericchi, & Varshavsky, 1998); choosing the constant to be $(2\pi {)}^{-q}$ is to simplify the resulting expression. With the definitions given in (3(3) $\begin{array}{rl}\hat{\mathit{I}}& =\left(\begin{array}{cc}{\hat{\mathit{I}}}_{11}& {\hat{\mathit{I}}}_{12}\\ {\hat{\mathit{I}}}_{21}& {\hat{\mathit{I}}}_{22}\end{array}\right),\text{and define }\\ {\mathbf{\Sigma }}^{-1}& ={\hat{\mathit{I}}}_{11}-{\hat{\mathit{I}}}_{12}{\hat{\mathit{I}}}_{22}^{-1}{\hat{\mathit{I}}}_{12}^t.\end{array}$(3) ), integrating out ${\mathit{\theta }}_{\left(2\right)}$ results in the expression $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}\int \exp \left(-\frac{1}{2}(\mathit{\theta }-\hat{\mathit{\theta }}{)}^t\hat{\mathit{I}}(\mathit{\theta }-\hat{\mathit{\theta }})\right)\\ & \times \left(2\pi {)}^{-q}d{\mathit{\theta }}_{\left(2\right)}\pi \right({\mathit{\theta }}_{\left(1\right)}\left)d{\mathit{\theta }}_{\left(1\right)}\right(1+o\left(1\right))\\ & =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \int \frac{\exp \left(-\frac{1}{2}({\mathit{\theta }}_{\left(1\right)}-{\hat{\mathit{\theta }}}_{\left(1\right)}{)}^t{\mathbf{\Sigma }}^{-1}({\mathit{\theta }}_{\left(1\right)}-{\hat{\mathit{\theta }}}_{\left(1\right)})\right)}{|\mathbf{\Sigma }{|}^{1/2}}\\ & \times \pi \left({\mathit{\theta }}_{\left(1\right)}\right)\mathrm{d}{\mathit{\theta }}_{\left(1\right)}(1+o(1\left)\right).\end{array}$

Step 2. Change variables to $\mathit{\xi }=\mathit{O}{\mathit{\theta }}_{\left(1\right)}$, where $\mathit{O}$ is an orthogonal matrix such that $\mathbf{\Sigma }={\mathit{O}}^t\mathit{D}\mathit{O}$, with $\mathit{D}=\text{diag}\left(d_i\right)$ for $i=1,\dots ,p$. (The choice of $\mathit{O}$ does not matter in the following.) Note that $\hat{\mathit{\xi }}=\mathit{O}{\hat{\mathit{\theta }}}_{\left(1\right)}$.

For this model, we will utilise a prior distribution that is independent in the ${\xi }_i$, i.e. $\pi \left(\mathit{\xi }\right)=\prod _{i=1}^p{\pi }_i\left({\xi }_i\right)$. Then we can write (8) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\int \frac{1}{\sqrt{d_i}}{e}^{-\left(({\xi }_i-{\hat{\xi }}_i{)}^2/2d_i\right)}{\pi }_i\left({\xi }_i\right)\mathrm{d}{\xi }_i\right]\\ & \times (1+o(1\left)\right).\end{array}$(8) For ${\pi }_i\left({\xi }_i\right)$, in a similar situation, Jeffreys (1961) recommended the Cauchy$(0,b_i)$ density $\left(1/\pi \sqrt{b_i}\right)\left(1/\right(1+{\xi }_i^2/b_i\left)\right)$, where $b_i$ is chosen to represent unit information for ${\xi }_i$ (see Kass & Wasserman, 1995; also to be discussed later). A prior that yields almost the same results is (9) ${\pi }_i^R\left({\xi }_i\right)={\int }_0^{1}N\left({\xi }_i|0,\frac{1}{2{\lambda }_i}(d_i+b_i)-d_i\right)\frac{1}{2\sqrt{{\lambda }_i}}d{\lambda }_i,$(9) which is well-defined if $b_i\ge d_i$. Interestingly, this prior is very similar to the Cauchy prior no matter what $d_i$ happens to be (as shown in the Appendix), so we will interpret this prior (and $b_i$) exactly as we would with the Cauchy prior. The attraction of ${\pi }^R$ is that the ensuing computations can be done in closed form. That one can have all the advantages that Jeffreys pointed out are possessed by the Cauchy prior for model selection, while maintaining closed form expressions, is a significant advantage when dealing with large model spaces. This prior was extensively discussed in Berger (1985), as a robust prior (hence the R label) for estimation problems, but its even greater value for model selection was not recognised. (This type of prior was first utilised in Strawderman (1971) in shrinkage estimation.) See also Bayarri et al. (2012), where a multivariate version of this prior is utilised for model selection in normal linear models.

With the prior in (9(9) ${\pi }_i^R\left({\xi }_i\right)={\int }_0^{1}N\left({\xi }_i|0,\frac{1}{2{\lambda }_i}(d_i+b_i)-d_i\right)\frac{1}{2\sqrt{{\lambda }_i}}d{\lambda }_i,$(9) ), the integral in (8(8) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\int \frac{1}{\sqrt{d_i}}{e}^{-\left(({\xi }_i-{\hat{\xi }}_i{)}^2/2d_i\right)}{\pi }_i\left({\xi }_i\right)\mathrm{d}{\xi }_i\right]\\ & \times (1+o(1\left)\right).\end{array}$(8) ) is straightforward to evaluate in closed form (first integrate over ${\xi }_i$, then over ${\lambda }_i$) yielding (10) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\frac{1}{\sqrt{(d_i+b_i)}}\frac{\left(1-e^{-{\hat{\xi }}_i^2/(d_i+b_i)}\right)}{\sqrt{2}{\hat{\xi }}_i^2/(d_i+b_i)}\right]\\ & \times (1+o(1\left)\right).\end{array}$(10)

Step 3. Define the unit information, $b_i$, by (11) $\begin{array}{rl}{b}_i& =n_i^e{d}_i,\text{where}\\ n_i^e& =\text{effective sample size for }{\xi }_i\text{ and recall }\\ v_i& =\frac{{\hat{\xi }}_i^2}{d_i(1+n_i^e)}.\end{array}$(11) Definitions of the effective sample size will be given in Section 3. It will be the case that $n_i^e\ge 1$ so that $b_i\ge d_i$ (the condition mentioned earlier for ${\pi }^R$ to be well defined). Then (10(10) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\frac{1}{\sqrt{(d_i+b_i)}}\frac{\left(1-e^{-{\hat{\xi }}_i^2/(d_i+b_i)}\right)}{\sqrt{2}{\hat{\xi }}_i^2/(d_i+b_i)}\right]\\ & \times (1+o(1\left)\right).\end{array}$(10) ) becomes $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}\frac{|\hat{\mathit{I}}{|}^{-1/2}}{|\mathit{D}{|}^{1/2}}\prod _{i=1}^p\\ & \times \frac{1}{\sqrt{1+n_i^e}}\frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i}(1+o(1\left)\right).\end{array}$ Since $|\hat{\mathit{I}}|=|{\mathbf{\Sigma }}^{-1}||{\hat{\mathit{I}}}_{22}|=|{\mathit{D}}^{-1}||{\hat{\mathit{I}}}_{22}|$, we thus have that $-2\log m\left(\mathbf{x}\right)=\text{PBIC}+o\left(1\right),$ with PBIC defined in (4(4) $\begin{array}{rl}\text{PBIC}& \equiv -2l\left(\hat{\mathit{\theta }}\right)+\log |{\hat{\mathit{I}}}_{22}|+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i},\end{array}$(4) ).

3. Defining ‘effective sample size’ $n_j$, for parameter ${\xi }_j$

The most difficult aspect of dealing with PBIC turns out to be defining the effective sample size corresponding to a parameter. We first present a solution for linear models, and then suggest a possible solution for the general case.

3.1. Effective sample sizes in linear models

Suppose that all models under consideration are linear models of the form (12) $\mathit{Y}={\mathit{X}}^{\ast }\mathit{\alpha }+\mathit{X}\mathit{\beta }+\mathit{\epsilon },\text{where }\mathit{\epsilon }\sim N(\mathbf{0},\mathbf{\Gamma }),\mathbf{\Gamma }\text{known},$(12) with dimensions ${\mathit{Y}}_{[n\times 1]},{\mathit{X}}_{[n\times q]}^{\ast },{\mathit{\alpha }}_{[q\times 1]},{\mathit{X}}_{[n\times p]},{\mathit{\beta }}_{[p\times 1]},{\mathit{\epsilon }}_{[n\times 1]}$ and ${\mathbf{\Gamma }}_{[n\times n]}$. Here ${\mathit{X}}^{\ast }\mathit{\alpha }$ is a common term present in all models (e.g. an intercept in linear regression), but $\mathit{X}\mathit{\beta }$ will differ from model to model. This fits into the framework for PBIC by defining ${\mathit{\theta }}_{\left(1\right)}=\mathit{\beta }$ and ${\mathit{\theta }}_{\left(2\right)}=\mathit{\alpha }$.

Since $\mathit{\alpha }$ will be integrated out in PBIC, only the effective sample size for linear functions of $\mathit{\beta }$ will be needed. The first step of the process is to orthogonalise the parameters by transforming $\mathit{\alpha }$ to ${\mathit{\alpha }}^{\ast }=\mathit{\alpha }+({{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1}{\mathit{X}}^{\ast }{)}^{-1}{{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1}\mathit{X}\mathit{\beta }$ and defining (13) $\tilde{\mathit{X}}=(\mathit{I}-{\mathit{X}}^{\ast }({{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1}{\mathit{X}}^{\ast }{)}^{-1}{{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1})\mathit{X}.$(13) Since ${\mathit{X}}^{\ast }\mathit{\alpha }+\mathit{X}\mathit{\beta }={\mathit{X}}^{\ast }{\mathit{\alpha }}^{\ast }+\tilde{\mathit{X}}\mathit{\beta }$, the linear part of the model has not changed in this reparameterisation, but now ${\tilde{\mathit{X}}}^t{\mathbf{\Gamma }}^{-1}{\mathit{X}}^{\ast }=\mathbf{0}$, so that ${\mathit{\alpha }}^{\ast }$ and $\mathit{\beta }$ are orthogonal. There are two important aspects of this. First, since ${\mathit{X}}^{\ast }$ has not been altered, the new ${\mathit{\alpha }}^{\ast }$ can still be considered common parameters in each model, and will be integrated out in PBIC, so that their changed definition is irrelevant. Second, $\mathit{\beta }$ has not been transformed, crucial because we wish effective sample sizes for linear functions of $\mathit{\beta }$

Write $\mathbf{\Gamma }=\mathit{\sigma }\mathit{R}\mathit{\sigma }$, with $\mathit{\sigma }=\text{diag}\{{\sigma }_1,\dots ,{\sigma }_p\}$, where $\mathit{R}$ is the correlation matrix, and define ${\mathit{C}}_{[p\times p]}$ to be the diagonal matrix with entries $c_{ii}=\underset{j}{max}\left\{|{\tilde{X}}_{ji}|/{\sigma }_j\right\}$. Berger, Bayarri, and Pericchi (2014) gave, as the general definition of the effective sample size (called TESS), for any scalar linear transformation $\xi =\mathit{v}\mathit{\beta }$ ($\mathit{v}$ is $[1\times p]$) of $\mathit{\beta }$, (14) $n^e=\frac{|\mathit{v}{|}^2}{\mathit{v}\mathit{C}({\tilde{\mathit{X}}}^t{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}{)}^{-1}\mathit{C}{\mathit{v}}^t}.$(14)

Example 3.1

Group means example

Assume $Y_{ij}={\mu }_i+{\epsilon }_{ij}$ for $i=1,\dots ,p$ groups, and $j=1,\dots ,r_i$ replicates in the ith group, and that the ${\epsilon }_{ij}$ are i.i.d. $N(0,{\sigma }^2)$. Computation yields that TESS for ${\mu }_i$ is $n_i^e=r_i$, as is to be expected. Note that $r_i$ could be 1, which can be seen to be the lower bound on TESS for linear models when $\mathbf{\Gamma }={\sigma }^2\mathit{I}$.

Example 3.2

Orthogonal and related designs

Assume that $\mathit{X}$ has orthogonal columns with entries $\pm a_i\ne 0$, and that $\mathbf{\Gamma }={\sigma }^2\mathit{I}$. Simple computation here shows that $n_i^e=n$ for each ${\beta }_i$.

Note that the effective sample size here is n, in contrast to the group means problem where the effective sample size can be as low as $r_i=1$. Indeed, it can be shown that, when $\mathbf{\Gamma }={\sigma }^2\mathit{I}$, TESS will always be between 1 and n, with both limits attainable.

Example 3.3

Heteroscedastic independent observations

Assume $Y_i=\mu +{\epsilon }_i$, ${\epsilon }_i$ independent, ${\epsilon }_i\sim N(0,{\sigma }_i^2),i=1,\dots ,n$. Here the effective sample size is $n^e=\frac{\sum _{i=1}^{n}1/{\sigma }_i^2}{\underset{i}{max}\left\{1/{\sigma }_i^2\right\}}.$ Consider the particular case where, for $i=1,\dots ,n_1$, we have $Y_i\sim N(\mu ,{\sigma }_1^2)$, whereas for the remaining $n_2=n-n_1$ observations, $Y_i\sim N(\mu ,{\sigma }_2^2)$, where ${\sigma }_2^2$ is much larger than ${\sigma }_1^2$; thus, intuitively, only the first $n_1$ observations count. Then, unless $n_2$ is large, $n^e=\frac{n_1/{\sigma }_1^2+n_2/{\sigma }_2^2}{1/{\sigma }_1^2}=n_1+n_2\frac{{\sigma }_1^2}{{\sigma }_2^2}\approx n_1.$

3.2. A general definition of effective sample size

Suppose one has independent observations $({\mathit{x}}_1,\dots ,{\mathit{x}}_n)$. A possible general definition for the ‘effective sample size’ follows from considering the information associated with observation ${\mathit{x}}_i$ arising from the single-observation expected information matrix ${\mathit{I}}_i^{\ast }={\mathit{O}}^{\mathrm{\prime }}\left(I_{i,jk}^{\ast }\right)\mathit{O}$, where $I_{i,jk}^{\ast }=-\text{E}\left[\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\theta }_j\mathrm{\partial }{\theta }_k}\log f_i({\mathit{x}}_i\mid \mathit{\theta })\right]{|}_{\mathit{\theta }=\hat{\mathit{\theta }}}.$ Since $I_{jj}^{\ast }=\sum _{i=1}^n{I}_{i,jj}^{\ast }$ is the expected information about ${\xi }_j$, a reasonable way to define the effective sample size, $n_j^e$, is

• define information weights $w_{ij}=I_{i,jj}^{\ast }/\sum _{k=1}^n{I}_{k,jj}^{\ast };$

• define the effective sample size for ${\xi }_j$ as $n_j^e=\frac{I_{jj}^{\ast }}{\sum _{i=1}^n{w}_{ij}{I}_{i,jj}^{\ast }}=\frac{{\left(I_{jj}^{\ast }\right)}^2}{\sum _{i=1}^n{\left(I_{i,jj}^{\ast }\right)}^2}.$

Intuitively, $\sum w_{ij}{I}_{i,jj}^{\ast }$ is a weighted measure of the information ‘per observation’, and dividing the total information about ${\xi }_j$ by this information per case seems plausible as an effective sample size.

Unfortunately, this does not seem to always be an effective definition; for instance, it does not reduce to TESS for all linear models. This should thus be viewed as primarily a starting point for future investigations of effective sample size in nonlinear models.

4. PBIC*: a version more favourable to complex models

Recall, from Raftery (1999), that BIC can be thought of as arising from unit information priors for each model that are centred at the model likelihood. This choice of prior seems highly favourable to more complex models, since the prior gives virtually all of its mass to a modest neighbourhood of the likelihood for each model.

In contrast, PBIC utilises unit information priors that are centred at 0 and, hence, can give little mass to the region of high model likelihood. The fat tails of the prior do result in reasonable answers (cf. Bayarri et al., 2012; Jeffreys, 1961), but it is of interest to investigate an intermediate solution.

The intermediate solution is to keep the prior centred at 0, but choose the scales of the prior, $b_i$, so that the prior will extend out to the likelihood. In our setup, this can be implemented by choosing the $b_i$ so as to maximise $m\left(\mathit{x}\right)$ in (10(10) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\frac{1}{\sqrt{(d_i+b_i)}}\frac{\left(1-e^{-{\hat{\xi }}_i^2/(d_i+b_i)}\right)}{\sqrt{2}{\hat{\xi }}_i^2/(d_i+b_i)}\right]\\ & \times (1+o(1\left)\right).\end{array}$(10) ); thus we are effectively choosing the prior in our class that is most favourable to each model. Clearly this does allow the prior to give more mass to the region of high model likelihood, but does not allow complete concentration of mass in this region.

Since this prior is maximising the marginal likelihood among the given class, it can be viewed as the empirical Bayes prior in the class. It was also a choice popularised in the ‘robust Bayes’ literature (cf. Berger & Sellke, 1987), and was used in Bollen et al. (2012) to develop a related generalisation of BIC.

The $b_i$ that maximises (10(10) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\frac{1}{\sqrt{(d_i+b_i)}}\frac{\left(1-e^{-{\hat{\xi }}_i^2/(d_i+b_i)}\right)}{\sqrt{2}{\hat{\xi }}_i^2/(d_i+b_i)}\right]\\ & \times (1+o(1\left)\right).\end{array}$(10) ) can easily be seen to be $\begin{array}{rl}{\hat{b}}_i& =max\{d_i,\frac{{\hat{\xi }}_i^2}{w}-d_i\},\text{with }w\text{s.t.}{e}^w-1\\ & =2w,\text{or}w\approx \mathrm{1.3.}\end{array}$ Unfortunately, the resulting version of BIC has serious problems; in particular it will typically not be consistent as $n\to \mathrm{\infty }$ in that, if ${\xi }_i$ is zero, the prior will concentrate about zero at such a fast rate that the models with and without ${\xi }_i$ are essentially equivalent (and one will fail to select the model without ${\xi }_i$ with probability approaching 1). This same lack of consistency afflicts the developments in Bollen et al. (2012) and the robust Bayesian choices.

The obvious solution is simply to prevent ${\hat{b}}_i$ from becoming too small, and the obvious constraint is to restrict it to the region $[n_i^e{d}_i,\mathrm{\infty })$. This yields the recommended choice (15) $b_i^{\ast }=max\left\{n_i^e{d}_i,\frac{{\hat{\xi }}_i^2}{1.3}-d_i\right\}.$(15) This will avoid inconsistency as $n\to \mathrm{\infty }$ in that, as long as $b_i^{\ast }\to c$ with c a non-zero constant, the resulting prior behaves asymptotically when ${\xi }_i=0$ as a fixed prior, and fixed priors will yield consistency as $n\to \mathrm{\infty }$. (Consistency when the effective sample sizes do not grow is a more delicate issue, discussed in Section 5.5.)

Replacing $b_i$ with $b_i^{\ast }$, (10(10) $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\\ & \times \left[\prod _{i=1}^p\frac{1}{\sqrt{(d_i+b_i)}}\frac{\left(1-e^{-{\hat{\xi }}_i^2/(d_i+b_i)}\right)}{\sqrt{2}{\hat{\xi }}_i^2/(d_i+b_i)}\right]\\ & \times (1+o(1\left)\right).\end{array}$(10) ) becomes $\begin{array}{rl}m\left(\mathbf{x}\right)& =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\left[\prod _{i=1}^p\frac{1}{\sqrt{d_i(1+n_i^e)max\{1,v_i/1.3\}}}\right.\\ & \times \left.\frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2}min\{v_i,1.3\}}\right](1+o(1\left)\right)\\ & =e^{l\left(\hat{\mathit{\theta }}\right)}|\hat{\mathit{I}}{|}^{-1/2}\left[\prod _{i=1}^p\frac{1}{\sqrt{d_i(1+n_i^e)}}\right.\\ & \times \left.\frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}}\right](1+o(1\left)\right).\end{array}$ The resulting approximation to $-2\log m\left(\mathit{x}\right)$ is given in (6(6) $\begin{array}{rl}{\text{PBIC}}^{\ast }& \equiv -2l\left(\hat{\mathit{\theta }}\right)+\log |{\hat{\mathit{I}}}_{22}|+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}}.\end{array}$(6) ).

5. PBIC and PBIC* for the linear model

5.1. The expressions

Consider the normal linear model framework in (12(12) $\mathit{Y}={\mathit{X}}^{\ast }\mathit{\alpha }+\mathit{X}\mathit{\beta }+\mathit{\epsilon },\text{where }\mathit{\epsilon }\sim N(\mathbf{0},\mathbf{\Gamma }),\mathbf{\Gamma }\text{known},$(12) ) and assume the orthogonalisation discussed there has been carried out. This does not change PBIC, but is more convenient because we can ignore the common orthogonal parameter ${\mathit{\alpha }}^{\ast }$ (see Bayarri et al., 2012; Berger et al., 1998; Jeffreys, 1961 for justification), and focus only on the other parameters $\mathit{\beta }$, with the associated model (16) $\mathit{Y}=\tilde{\mathit{X}}\mathit{\beta }+\mathit{\epsilon },\text{where}\mathit{\epsilon }\sim N(\mathbf{0},\mathbf{\Gamma }),\mathbf{\Gamma }\text{known},$(16) with $\tilde{\mathit{X}}$ given by (13(13) $\tilde{\mathit{X}}=(\mathit{I}-{\mathit{X}}^{\ast }({{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1}{\mathit{X}}^{\ast }{)}^{-1}{{\mathit{X}}^{\ast }}^t{\mathbf{\Gamma }}^{-1})\mathit{X}.$(13) ).

Following the PBIC algorithm, note that ${\mathbf{\Sigma }}^{-1}={\tilde{\mathit{X}}}^{\prime }{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}$. Change variables to $\mathit{\xi }=\mathit{O}\mathit{\beta }$, where $\mathit{O}$ is an orthogonal matrix such that $\mathbf{\Sigma }={\mathit{O}}^t\mathit{D}\mathit{O}$, with $\mathit{D}=\text{diag}\left(d_i\right)$ for $i=1,\dots ,p$. Then, for each ${\xi }_j={\mathit{O}}_j\mathit{\beta }$, define $n_j^e$ using (14(14) $n^e=\frac{|\mathit{v}{|}^2}{\mathit{v}\mathit{C}({\tilde{\mathit{X}}}^t{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}{)}^{-1}\mathit{C}{\mathit{v}}^t}.$(14) ) with $\mathit{v}={\mathit{O}}_j$, and let ${\hat{\xi }}_j={\mathit{O}}_j\hat{\mathit{\beta }}$, where $\hat{\mathit{\beta }}=({\tilde{\mathit{X}}}^{\prime }{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}{)}^{-1}{\tilde{\mathit{X}}}^{\prime }{\mathbf{\Gamma }}^{-1}\mathit{Y}$. Finally, recalling that $v_i={\hat{\xi }}_i^2/\left[d_i\right(1+n_i^e\left)\right]$, PBIC and PBIC* are given by (17) $\begin{array}{rl}\text{PBIC}& =S^2+C+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i}\end{array}$(17) (18) $\begin{array}{rl}\text{PBIC*}& =S^2+C+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}},\end{array}$(18) where $S^2$ is the usual residual sum of squares corresponding to (12(12) $\mathit{Y}={\mathit{X}}^{\ast }\mathit{\alpha }+\mathit{X}\mathit{\beta }+\mathit{\epsilon },\text{where }\mathit{\epsilon }\sim N(\mathbf{0},\mathbf{\Gamma }),\mathbf{\Gamma }\text{known},$(12) ) and $C=\log \left(|\mathbf{\Gamma }|\right)+n\log \left(2\pi \right)=\log \left(|\mathbf{\Gamma }|\right)+n\log \left(2\pi \right).$ Note that C is the same constant in any model under consideration, and hence it can be ignored in comparing models or determining Bayes factors.

In what follows we describe some important Linear Model examples. There are more, including correlated observations and autoregressive models, in Berger et al. (2014).

5.2. Simple linear regression

Let $Y_i=\alpha +X_i\beta +{\epsilon }_i,{\epsilon }_i\stackrel{i.i.d.}{\sim }N(0,{\sigma }^2)$, so that $\begin{array}{rl}\mathit{Y}& =\left(\begin{array}{cc}1& X_1\\ ⋮& ⋮\\ 1& X_n\end{array}\right)\left(\begin{array}{c}\alpha \\ \beta \end{array}\right)\\ & +\left(\begin{array}{c}{\epsilon }_1\\ ⋮\\ {\epsilon }_n\end{array}\right),\text{where}\mathit{\epsilon }\sim N(\mathbf{0},{\sigma }^2\mathit{I}).\end{array}$ Suppose we are considering two models $M_0:\beta =0$ and $M_1:\beta \ne 0$. Computation under $M_1$ yields $\tilde{\mathit{X}}=(X_1-\overline{X},\dots ,X_n-\overline{X}{)}^{\prime }$, so that $\mathrm{\Sigma }={\sigma }^2/s_x^2={\sigma }^2/\sum _{i=1}^n(X_i-\overline{X}{)}^2$. Also, from (14(14) $n^e=\frac{|\mathit{v}{|}^2}{\mathit{v}\mathit{C}({\tilde{\mathit{X}}}^t{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}{)}^{-1}\mathit{C}{\mathit{v}}^t}.$(14) ), (19) $n^e=\frac{\sum _{i=1}^n(X_i-\overline{X}{)}^2}{\underset{i}{max}\left\{\right(X_i-\overline{X}{)}^2\}}=\frac{s_x^2}{\underset{i}{max}\left\{\right(X_i-\overline{X}{)}^2\}}.$(19) Finally, $d=\mathrm{\Sigma }={\sigma }^2/s_x^2$, $v={\hat{\beta }}^2/\left[d\right(1+n^e\left)\right]$, and $\begin{array}{rl}{S}^2& =\frac{1}{{\sigma }^2}\left(|\mathit{Y}{|}^2-\frac{\left(\sum _{i=1}^n\right(x_i-\overline{x})y_i{)}^2}{\sum _{i=1}^n(x_i-\overline{x}{)}^2}\right)\\ & =\frac{1}{{\sigma }^2}(|\mathit{Y}{|}^2-s_x^2{\hat{\beta }}^2)\end{array}$ complete the terms needed to define PBIC and PBIC* under $M_1$. Under $M_0$, we only need $S^2=\left(1/{\sigma }^2\right)|\mathit{Y}{|}^2$; thus, with $v={\hat{\beta }}^2/\left[{\sigma }^2\right(s_x^{-2}+\left(\underset{i}{max}\right\{(X_i-\overline{X}{)}^2\}{)}^{-1}\left)\right]$, $\begin{array}{rl}\mathrm{\Delta }\text{PBIC}& =-\frac{s_x^2{\hat{\beta }}^2}{{\sigma }^2}+\log \left(1+\frac{s_x^2}{\underset{i}{max}\left\{\right(X_i-\overline{X}{)}^2\}}\right)\\ & -2\log \frac{\left(1-e^{-v}\right)}{\sqrt{2}v}.\end{array}$$\mathrm{\Delta }\text{PBIC*}$ is the obvious modification of this.

5.3. Testing equality of two means with unequal variances

Consider comparing two normal means via the test $H_0:{\mu }_1={\mu }_2$ versus $H_1:{\mu }_1\ne {\mu }_2$, where the associated known variances, ${\sigma }_1^2$ and ${\sigma }_2^2$ are not equal. The linear model is thus $\begin{array}{rl}\mathit{Y}& =\mathit{X}\mathit{\mu }+\mathit{\epsilon }=\left(\begin{array}{cc}1& 0\\ ⋮& ⋮\\ 1& 0\\ 0& 1\\ ⋮& ⋮\\ 0& 1\end{array}\right)\left(\begin{array}{c}{\mu }_1\\ {\mu }_2\end{array}\right)+\left(\begin{array}{c}{\epsilon }_{11}\\ ⋮\\ {\epsilon }_{2n_2}\end{array}\right),\\ & \times \mathit{\epsilon }\sim N(\mathbf{0},\text{diag}\{\underset{n_1}{\underbrace{{\sigma }_1^2,\dots ,{\sigma }_1^2}},\underset{n_2}{\underbrace{{\sigma }_2^2,\dots ,{\sigma }_2^2}}\left\}\right).\end{array}$ Defining $\alpha =({\mu }_1+{\mu }_2)/2$ and $\beta =({\mu }_1-{\mu }_2)/2$ places this in the linear model comparison framework, where we are comparing $M_0:\beta =0$ versus $M_1:\beta \ne 0$ with the covariate matrix $\mathit{B}=\mathit{X}{\left(\begin{array}{cc}1& 1\\ 1& -1\end{array}\right)}^{-1}=\frac{1}{2}\left(\begin{array}{cc}1& 1\\ ⋮& ⋮\\ 1& 1\\ 1& -1\\ ⋮& ⋮\\ 1& -1\end{array}\right).$ Under $M_1$, computation yields $\begin{array}{rl}\tilde{\mathit{X}}& ={\left(\frac{n_2}{n^{\ast }{\sigma }_2^2},\dots ,\frac{n_2}{n^{\ast }{\sigma }_2^2},-\frac{n_1}{n^{\ast }{\sigma }_1^2},\dots ,-\frac{n_1}{n^{\ast }{\sigma }_1^2}\right)}^{\prime }\text{with}\\ n^{\ast }& ={\left(\frac{n_1}{{\sigma }_1^2}+\frac{n_2}{{\sigma }_2^2}\right)}^{-1},\end{array}$ so that $d=\mathrm{\Sigma }=\left(\frac{{\sigma }_1^2}{n_1}+\frac{{\sigma }_2^2}{n_2}\right).$ Also, from (14(14) $n^e=\frac{|\mathit{v}{|}^2}{\mathit{v}\mathit{C}({\tilde{\mathit{X}}}^t{\mathbf{\Gamma }}^{-1}\tilde{\mathit{X}}{)}^{-1}\mathit{C}{\mathit{v}}^t}.$(14) ), $\begin{array}{rl}{n}^e& =\frac{\left(\frac{{\sigma }_1^2}{n_1}+\frac{{\sigma }_2^2}{n_2}\right)}{max\left\{{\sigma }_1^2/n_1^2,{\sigma }_2^2/n_2^2\right\}}\\ & =min\left\{\frac{n_1^2}{{\sigma }_1^2},\frac{n_2^2}{{\sigma }_2^2}\right\}\left(\frac{{\sigma }_1^2}{n_1}+\frac{{\sigma }_2^2}{n_2}\right),\end{array}$ and $v={\hat{\beta }}^2/\left[d\right(1+n^e\left)\right]$.

A special case is the standard test of equality of means when ${\sigma }_1^2={\sigma }_2^2={\sigma }^2$. Then $n^e=min\left\{n_1\left(1+\frac{n_1}{n_2}\right),n_2\left(1+\frac{n_2}{n_1}\right)\right\}.$ While this may look unusual, looking at the extremes indicates why this is reasonable. Indeed, as say $n_1\to \mathrm{\infty }$, note that $n^e\to n_2$. In this scenario, we perfectly learn ${\mu }_1$, so the test of mean equality is really just a test that ${\mu }_2$ equals this known mean, based on $n_2$ observations. Attempting to utilise BIC with an adhoc choice of n, such as $(n_1+n_2)/2$, would clearly be a disaster here.

5.4. Findley's counterexample to BIC

For the simple linear model in (2(2) $\begin{array}{rl}{X}_i& =\frac{1}{\sqrt{i}}\cdot \theta +{\epsilon }_i,\text{where }{\epsilon }_i\sim N(0,1),\\ i& =1,\dots ,n,\end{array}$(2) ), computation yields that, under $H_1:\theta \ne 0$, $\begin{array}{rl}d& =\mathrm{\Sigma }={\left(\sum _{i=1}^n\frac{1}{i}\right)}^{-1},n^e=\sum _{i=1}^n\frac{1}{i},\\ S^2& =|\mathit{Y}{|}^2-{\hat{\theta }}^2\sum _{i=1}^n\frac{1}{i}.\end{array}$ It follows that $\begin{array}{rl}\mathrm{\Delta }\text{PBIC}& =-{\hat{\theta }}^2\sum _{i=1}^n\frac{1}{i}+\log (1+\sum _{i=1}^n\frac{1}{i})\\ & -2\log \frac{\left(1-e^{-v}\right)}{\sqrt{2}v},v=\frac{{\hat{\theta }}^2}{d(1+n^e)}.\end{array}$ Since $\sum _{i=1}^n\left(1/i\right)=\log n+O\left(1\right)$ and ${\hat{\theta }}^2\to {\theta }^2$ (because the mle is consistent here), $\begin{array}{rl}\mathrm{\Delta }\text{PBIC}& =-{\theta }^2(\log n+O(1\left)\right)+\log (\log n)\\ & -2\log \frac{\left(1-e^{-{\theta }^2}\right)}{\sqrt{2}{\theta }^2}+o\left(1\right).\end{array}$ Under $H_0:\theta =0$, $\mathrm{\Delta }\text{PBIC}=\log (\log n)+\log 2+o\left(1\right)\to \mathrm{\infty }$ and, under $H_1:\theta \ne 0$, $\mathrm{\Delta }\text{PBIC}\to -\mathrm{\infty }$. Thus PBIC is consistent as $n\to \mathrm{\infty }$. Essentially the same argument shows that PBIC* is consistent.

5.5. Consistency of PBIC and PBIC* as $p\to \mathrm{\infty }$ in the group means problem

Bayes model selection rules for fixed priors and fixed p are virtually always consistent as the sample size $n\to \mathrm{\infty }$. This type of consistency transfers over to rules such as PBIC and PBIC* because the priors from which they arise converge to fixed priors as $n\to \mathrm{\infty }$ with p fixed.

There is nothing within Bayesian theory, however, that guarantees consistency of Bayes rules when the dimension p also grows. Indeed, it turns out that consistency is then a very delicate property, that can easily be violated by even standard Bayes rules. The group means problem provides a simple illustration.

Example 5.1

Consider the group means problem with known ${\sigma }^2=1$ and effective sample size $n_i=r$ fixed, and reduce to the sufficient statistics ${\overline{X}}_i\sim N({\mu }_i,1/r)$ for $i=1,\dots ,p$. Consider comparison of the null model $M_0:{\mu }_1=\cdots ={\mu }_p=0$ with the full model $M_1:\text{ all }{\mu }_i\text{ non-zero}$. Suppose the ${\mu }_i$ are independently assigned $N(0,{\tau }_i^2)$ priors. Then it is easy to show that consistency obtains under $M_1$ as $p\to \mathrm{\infty }$ if and only if $V\equiv \underset{p\to \mathrm{\infty }}{lim}\left(1/p\right)\sum _i^p{\overline{X}}_i^2$ satisfies $V\ge \underset{p\to \mathrm{\infty }}{lim}\left(1/p\right)\sum _i^p{\tau }_i^2$, assuming the limits exist. (This example was brought to our attention by J. K. Ghosh.)

After reflecting upon this, it might seem surprising that any prior could achieve consistency as $p\to \mathrm{\infty }$. However, Berger et al. (2003) computed Laplace approximations to the marginal density, for this problem, that produced consistent Bayes factors when p grows with n. They used a multivariate Cauchy prior, which does not result in a closed form Bayes factor, as arises with PBIC and PBIC*. The next theorem indicates the situation involving consistency for PBIC and PBIC*.

Theorem 5.2

For the group means problem with fixed r, PBIC and PBIC* are consistent under $M_0$ as $p\to \mathrm{\infty }$. Under $M_1$ and assuming that ${\tau }^2=\underset{p\to \mathrm{\infty }}{lim}\left(1/p\right)\sum _i^p{\mu }_i^2$ exists, PBIC and PBIC* are (20) $\begin{array}{rl}& \text{consistent if }{\tau }^2\\ & >\frac{\log 2+\log (1+r)+1}{r};\text{inconsistent if }{\tau }^2\\ & <\frac{\log 2+\log (1+r)-1}{r}.\end{array}$(20)

Proof.

We utilise (17(17) $\begin{array}{rl}\text{PBIC}& =S^2+C+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-v_i}\right)}{\sqrt{2}{v}_i}\end{array}$(17) ) and (18(18) $\begin{array}{rl}\text{PBIC*}& =S^2+C+\sum _{i=1}^p\log (1+n_i^e)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}},\end{array}$(18) ) as the definitions of PBIC and PBIC*, but will ignore C since it is common to all models. Note that the $n_i^e=r$, $S_1^2=\sum _{i=1}^p\sum _{j=1}^r(x_{ij}-{\overline{x}}_i{)}^2$, $S_0^2=S_1^2+r\sum _{i=1}^p{\overline{x}}_i^2$, $v_i=r{\overline{x}}_i^2/(r+1)$ under $M_1$ and $v_i=0$ under $M_0$. Thus PBIC and PBIC* become, with subscripts denoting the model, $\begin{array}{rl}{\text{PBIC}}_0& ={\text{PBIC*}}_0=S_0^2=S_1^2+r\sum _{i=1}^p{\overline{x}}_i^2,\\ {\text{PBIC}}_1& =S_1^2+p\log (r+1)-2\sum _{i=1}^p\log \frac{1-e^{-v_i}}{\sqrt{2}{v}_i}\\ & =S_1^2+p\log \left[2\right(r+1\left)\right]-2\sum _{i=1}^p\log \frac{1-e^{-v_i}}{v_i},\\ {\text{PBIC*}}_1& =S_1^2+p\log (r+1)\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{2v_{i}min\{v_i,1.3\}}}\\ & =S_1^2+p\log \left[2\right(r+1\left)\right]\\ & -2\sum _{i=1}^p\log \frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{v_{i}min\{v_i,1.3\}}}.\end{array}$ It is straightforward to show that $\frac{1-e^{-v_i}}{v_i}<1\text{and}\frac{\left(1-e^{-min\{v_i,1.3\}}\right)}{\sqrt{v_{i}min\{v_i,1.3\}}}<1,$ so that $\mathrm{\Delta }\text{PBIC}={\text{PBIC}}_1-{\text{PBIC}}_0$ and $\mathrm{\Delta }\text{PBIC*}={\text{PBIC*}}_1-{\text{PBIC*}}_0$ satisfy $\begin{array}{rl}& \mathrm{\Delta }\text{PBIC}\left(\mathrm{\Delta }\text{PBIC*}\right)\\ & >p\log \left[2\right(r+1\left)\right]-r\sum _{i=1}^p{\overline{x}}_i^2\equiv A\left(p\right).\end{array}$ Under $M_0$, $r\sum _{i=1}^p{\overline{x}}_i^2\sim {\chi }_p^2$, so that $\begin{array}{rl}A\left(p\right)& =p\log \left[2\right(r+1\left)\right]\\ & -p\left(1+O\left(\frac{1}{\sqrt{p}}\right)\right)\to \mathrm{\infty }\text{as}p\to \mathrm{\infty },\end{array}$ establishing consistency under $M_0$.

To show inconsistency under $M_1$, note that $r\sum _{i=1}^p{\overline{x}}_i^2\sim {\chi }_p^2\left({\lambda }_p\right)$, with non-centrality parameter ${\lambda }_p=r\sum _{i=1}^p{\mu }_i^2$. Thus $\begin{array}{rl}A\left(p\right)& =p\log \left[2\right(r+1\left)\right]\\ & -(p+{\lambda }_p)\left(1+O\left(\frac{1}{\sqrt{p+{\lambda }_p}}\right)\right)\to \mathrm{\infty }\end{array}$ if ${\tau }^2=\underset{p\to \mathrm{\infty }}{lim}{\lambda }_p/\left[rp\right]<(\log [2(1+r)-1\left]\right)/r$, establishing the inconsistency result.

To investigate consistency of PBIC and PBIC* under $M_1$, note that (21) $\frac{\left(1-e^{-min\{v,1.3\}}\right)}{\sqrt{vmin\{v,1.3\}}}\ge \frac{1-e^{-v}}{v}\ge \frac{1}{1+v}.$(21) Also, because of concavity, $E[\log (1+v\left)\right]\le \log (1+E[v\left]\right)$. Thus, $\begin{array}{rl}E\left[\log \frac{1-e^{-v_i}}{v_i}\right]& \ge -E[\log (1+v_i\left)\right]\\ & \ge -\log (1+E[v_i\left]\right)\\ & =-\log \left(1+\frac{1+r{\mu }_i^2}{1+r}\right).\end{array}$ Using this inequality and the fact that $\prod _i{\omega }_i^{1/p}\le (\sum {\omega }_i)/p)$, it follows that $\begin{array}{rl}& \frac{2}{p}E\left[\sum _{i=1}^p\log \frac{1-e^{-v_i}}{v_i}\right]\\ & \ge -2\log \prod _{i=1}^p{\left(1+\frac{1+r{\mu }_i^2}{1+r}\right)}^{1/p}\\ & \ge -2\log \left[\frac{1}{p}\sum _{i=1}^p\left(1+\frac{1+r{\mu }_i^2}{1+r}\right)\right]\\ & =-2\log \left(\frac{2+r+{\lambda }_p/p}{1+r}\right).\end{array}$ Hence, by the law of large numbers, $\begin{array}{rl}& \underset{p\to \mathrm{\infty }}{lim}\frac{1}{p}\mathrm{\Delta }\text{PBIC}\\ & \le \log \left[2\right(r+1\left)\right]+2\log \left(\frac{2+r+r{\tau }^2}{1+r}\right)\\ & -\underset{p\to \mathrm{\infty }}{lim}\left(1+\frac{{\lambda }_p}{p}\right)\left(1+O\left(\frac{1}{\sqrt{p+\lambda }}\right)\right)\\ & =\log \left[2\right(r+1\left)\right]+2\log \left(\frac{2+r+r{\tau }^2}{1+r}\right)\\ & -(1+r{\tau }^2)(1+o(1\left)\right).\end{array}$ Let $B(r,{\tau }^2)$ denote the right hand side above (without the $o\left(1\right)$ term). If $B(r,{\tau }^2)<0$, then $\mathrm{\Delta }\text{PBIC}$ goes to $-\mathrm{\infty }$ as $p\to \mathrm{\infty }$, and we have consistency.

Differentiating with respect to ${\tau }^2$ shows that $B(r,{\tau }^2)$ is decreasing in ${\tau }^2$, so that, if we can find a value of ${\tau }^2$ for which $B(r,{\tau }^2)<0$, then any larger value of ${\tau }^2$ will also work. As a candidate, consider ${\tau }_c^2=[c+\log (1+r\left)\right]/r$. Then $\begin{array}{rl}B(r,{\tau }_c^2)& =\log \left[2\right(r+1\left)\right]\\ & +2\log \left(\frac{2+r+c+\log (1+r)}{1+r}\right)\\ & -(1+c+\log (1+r\left)\right).\end{array}$ Differentiating this with respect to r shows that it is decreasing in r so that all we need to show is that ${\tau }_c^2$ works for r=1. Indeed, $\begin{array}{rl}B(1,{\tau }_c^2)& =\log \left[4\right]+2\log \left(\frac{3+c+\log 2}{2}\right)\\ & -(1+c+\log 2)<0,\end{array}$ if c>1.67. Since $1+\log 2=1.693>1.67$, the condition for consistency of PBIC in the theorem is established. And because of (21(21) $\frac{\left(1-e^{-min\{v,1.3\}}\right)}{\sqrt{vmin\{v,1.3\}}}\ge \frac{1-e^{-v}}{v}\ge \frac{1}{1+v}.$(21) ), the same condition ensures that PBIC* is consistent.

Note that, if r is moderately large, PBIC and PBIC* are consistent under $M_1$, unless ${\tau }^2$ is extremely close to 0, i.e. unless the non-zero means are extremely close to 0; it is not surprising that it is difficult to distinguish between $M_1$ and $M_0$ in this situation.

There is a gap in the theorem between the consistency and inconsistency conditions under $M_1$. The gap is quite large for small r, but shrinks as r grows. A more refined analysis would reduce the gap, but the theorem does convey the basic messages about consistency.

More generally, $M_0$ could be a group means model containing some zero and some non-zero means. If $M_0$ is nested in $M_1$ and the number of additional non-zero means in $M_1$ goes to ∞, then the theorem still applies, since the common non-zero means will be integrated out at the beginning and will not affect the analysis.

Acknowledgements

This research was begun under the auspices of the 2004–2005 SAMSI program on Latent Variables in the Social Sciences.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

M. J. Bayarri's research was supported by the Spanish Ministry of Education and Science [grant number MTM2010-19528]; James Berger's research was supported by USA National Science Foundation [grant numbers DMS-1007773 and DMS-1407775]; Woncheol Jang's research was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea government (MSIP), No. 2014R1A4A1007895 and No. 2017R1A2B2012816; Luis Pericchi's research was supported by grant CA096297/CA096300 from the USA National Cancer Institute of the National Institutes of Health.

Appendix

To see that the prior in (9(9) ${\pi }_i^R\left({\xi }_i\right)={\int }_0^{1}N\left({\xi }_i|0,\frac{1}{2{\lambda }_i}(d_i+b_i)-d_i\right)\frac{1}{2\sqrt{{\lambda }_i}}d{\lambda }_i,$(9) ) is almost the same as ${\pi }^C$, the Cauchy(0,b) prior (we drop the i subscripts in this appendix), consider the extremes.

Theorem A.1

For $b\ge d,$ $\begin{array}{rl}\underset{|\xi |\to \mathrm{\infty }}{lim}\frac{{\pi }^C\left(\xi \right)}{{\pi }^R\left(\xi \right)}& =\frac{2\sqrt{b}}{\sqrt{\pi (b+d)}}\in (0.80,1.13),\\ \frac{{\pi }^C\left(0\right)}{{\pi }^R\left(0\right)}& =\frac{2d}{\sqrt{b\pi }(\sqrt{b+d}-\sqrt{b-d})}\\ & \in (0.80,1.13).\end{array}$

Proof.

Note that $\begin{array}{rl}{\pi }^R\left(0\right)& =\frac{1}{2\sqrt{\pi }}{\int }_0^1\frac{1}{\sqrt{d+b-2\lambda d}}\\ d\lambda & =\frac{\sqrt{b+d}-\sqrt{b-d}}{2d\sqrt{\pi }}.\end{array}$ Hence $\frac{{\pi }^C\left(0\right)}{{\pi }^R\left(0\right)}=\frac{2d}{\sqrt{b\pi }(\sqrt{b+d}-\sqrt{b-d})}.$ It is straightforward to show that $\sqrt{b}(\sqrt{b+d}-\sqrt{b-d})$ is decreasing in $b\ge d,$ with a maximum of $\sqrt{2}d$ and minimum of d. Thus $\sqrt{2/\pi }\le {\pi }^C\left(0\right)/{\pi }^R\left(0\right)\le \sqrt{4/\pi }$, which (to 2 decimal places) is the result above.

To prove the result as $|\xi |\to \mathrm{\infty }$, separately integrate over ${\mathrm{\Gamma }}_1=(0,|\xi {|}^{-3/2})$ and ${\mathrm{\Gamma }}_2=(|\xi {|}^{-3/2},1)$ in (9(9) ${\pi }_i^R\left({\xi }_i\right)={\int }_0^{1}N\left({\xi }_i|0,\frac{1}{2{\lambda }_i}(d_i+b_i)-d_i\right)\frac{1}{2\sqrt{{\lambda }_i}}d{\lambda }_i,$(9) ). For $\lambda \in {\mathrm{\Gamma }}_1$, note that $(d+b-2\lambda d{)}^{-1}=(d+b{)}^{-1}+O\left(|\xi {|}^{-3/2}\right)$, so that $\begin{array}{rl}(d+b-2\lambda d{)}^{-1/2}& =(d+b{)}^{-1/2}+O(|\xi {|}^{-3/2}),\\ \exp \left(-\frac{{\xi }^2\lambda }{d+b-2\lambda d}\right)& =\exp \left(-\frac{{\xi }^2}{d+b}\left[\lambda +O\left(|\xi {|}^{-3}\right)\right]\right)\\ & =\exp \left(-\frac{{\xi }^2}{d+b}\right)\left(1+O\left(|\xi {|}^{-1}\right)\right).\end{array}$ Hence the integral over ${\mathrm{\Gamma }}_1$ is $\begin{array}{rl}& \frac{1}{2\sqrt{\pi }}{\int }_0^{|\xi {|}^{-3/2}}\left(\frac{1}{\sqrt{d+b}}+O\left(|\xi {|}^{-3/2}\right)\right)\\ & \times \exp \left(-\frac{{\xi }^2\lambda }{d+b}\right)\left(1+O(|\xi {|}^{-1}\right)\mathrm{d}\lambda \\ & =\frac{\sqrt{d+b}}{2\sqrt{\pi }{\xi }^2}\left(1-\exp \left(-\frac{\sqrt{|\xi |}}{d+b}\right)\right)\left(1+O(|\xi {|}^{-1}\right).\end{array}$

Noting that $\exp (-{\xi }^2\lambda /[d+b-2\lambda d\left]\right)$ is decreasing in λ, it is immediate that the integral over ${\mathrm{\Gamma }}_2$ is bounded above by $\begin{array}{rl}& \frac{\exp (-\sqrt{|\xi |}/[d+b\left]\right)}{2\sqrt{\pi }d}{\int }_{|\xi {|}^{-3/2}}^1\frac{1}{\sqrt{d+b-2\lambda d}}\mathrm{d}\lambda \\ & =\frac{\exp (-\sqrt{|\xi |}/[d+b\left]\right)}{2\sqrt{\pi }d}\\ & \times \left(\sqrt{d+b-2|\xi {|}^{-3/2}d}-\sqrt{b-d}\right)=o\left(|\xi {|}^{-2}\right).\end{array}$

It follows that $\begin{array}{rl}& \underset{|\xi |\to \mathrm{\infty }}{lim}\frac{{\pi }^C\left(\xi \right)}{{\pi }^R\left(\xi \right)}=\underset{|\xi |\to \mathrm{\infty }}{lim}\frac{2\sqrt{\pi }\left[{\pi }^{-1}{\xi }^{-2}\sqrt{b}\right(1+O\left(|{\xi }^{-2}\right)\left)\right]}{\begin{array}{c}\sqrt{d+b}{\xi }^{-2}(1-\exp (-\sqrt{|\xi |}/[d+b\left]\right))\\ (1+O(|{\xi }^{-1}\left)\right)+o\left(|\xi {|}^{-2}\right)\end{array}}\\ & =\frac{2\sqrt{b}}{\sqrt{\pi (b+d)}}.\end{array}$

It is straightforward to show that $\sqrt{2/\pi }\le 2\sqrt{b}/\sqrt{\pi (b+d)}\le \sqrt{4/\pi }$, yielding (to two decimal places) the conclusion.