Laplace Power-expected-posterior priors for generalized linear models with
applications to logistic regression

Anupreet Porwal and Abel Rodriguez

Department of Statistics, University of Washington Seattle, WA, 98195, USA

###### Abstract

Power-expected-posterior (PEP) methodology, which borrows ideas from the
literature on power priors, expected-posterior priors and unit information
priors, provides a systematic way to construct objective priors. The basic
idea is to use imaginary training samples to update a noninformative prior
into a minimally-informative prior. In this work, we develop a novel
definition of PEP priors for generalized linear models that relies on a
Laplace expansion of the likelihood of the imaginary training sample. This
approach has various computational, practical and theoretical advantages over
previous proposals for non-informative priors for generalized linear models.
We place a special emphasis on logistic regression models, where sample
separation presents particular challenges to alternative methodologies. We
investigate both asymptotic and finite-sample properties of the procedures,
showing that is both asymptotic and intrinsic consistent, and that its
performance is at least competitive and, in some settings, superior to that of
alternative approaches in the literature.

##  1 Introduction

Generalized linear models (GLMs, e.g., see McCullagh & Nelder, 2019) are one
of the main workhorses of statistical analysis. Indeed, GLMs are widely used
both to model data directly and as building blocks for more complex
hierarchical models. However, in spite of their broad adoption, prior
elicitation for general GLMs in the absence of subjective information remains
an open problem, particularly in settings where the main goal is variable
selection. Because standard non-informative priors for GLMs that work well for
parameter estimation are often improper, they cannot be used in model
selection problems because they typically lead to ill-defined Bayes factors
(e.g., see Berger et al., 2001).

Within the subclass of Gaussian linear models, the literature on so-called
“objective” or “default” priors is extensive. Examples include point-mass
spike-&-slab priors (Mitchell & Beauchamp, 1988; Geweke, 1996), $g$-priors
(Zellner, 1986), mixtures of $g$-priors (Zellner & Siow, 1980; Liang et al.,
2008), unit information priors (Kass & Wasserman, 1995), intrinsic Bayes
factors (Berger & Pericchi, 1996a), fractional Bayes factors (O’Hagan, 1995;
De Santis & Spezzaferri, 2001), non-local priors (Johnson & Rossell, 2010,
2012) and power-expected-posterior priors (Fouskakis et al., 2015), among
other approaches. See Bayarri et al. (2012) for a review and a discussion of
desirable properties. The literature on default priors for GLMs is more
limited, with three main approaches dominating. These include those introduced
by Bové et al. (2011) and Li & Clyde (2018), both of which consider
modifications of mixtures of g-priors that are suitable for GLMs, and
Fouskakis et al. (2018), who considers extensions of power-expected-posterior
priors that relies on unnormalized power likelihoods. One feature shared by
all three approaches is that they can be thought of as being based on the idea
of calibrating (possibly improper) priors using either real or imaginary
training samples (e.g., see Berger & Pericchi, 1996b and Pérez & Berger,
2002).

In this paper we introduce a variant of the power-expected-posterior (PEP)
prior for GLMs that we call the Laplace PEP, or LPEP. While the formulation is
general, this manuscript emphasizes the development of the LPEP for logistic
regression models. This is because this subclass of models provides the best
illustration of the theoretical and practical advantages of our approach over
existing ones. For example, we note that the prior described in Li & Clyde
(2018) is improper when the maximum likelihood estimator of the regression
coefficients under the observed data does not exist. In the case of logistic
regression, this happens when there is separation among the groups (e.g., see
Albert & Anderson, 1984, Lesaffre & Albert, 1989 and Heinze & Schemper, 2002).
Separation is reasonably common in practical applications, especially in
problems with relatively small samples and several unbalanced and highly
predictive risk factors (e.g., see Section 6.1). A similar issue arises with
the PEP priors introduced in Fouskakis et al. (2018) since the imaginary
training samples are not restricted to yield finite maximum likelihood
estimators. Furthermore, both versions of the PEP prior proposed by Fouskakis
et al. (2018) are computationally intractable, requiring the use of reversible
Jump Markov Chain Monte Carlo algorithms (Green, 1995; Dellaportas et al.,
2002). Our LPEP is carefully designed to deal with separation issues, and it
is well defined as long as at least one training sample exists that yields
finite maximum likelihood estimates under the full model. Furthermore, the
LPEP can be written as a location-and-scale mixture of Gaussian priors and is
therefore computationally tractable and easy to incorporate into standard
Markov chain Monte Carlo algorithms that rely on data augmentation (e.g.,
Polson et al., 2013). Finally, the mixture structure also simplifies the
theoretical study of the prior, allowing us to show that the procedure is both
asymptotically consistent and intrinsically consistent.

It is important to stress that the focus of this manuscript is on priors for
variable selection that place positive probability on specific coefficients
being exactly zero. An alternative approach is to use continuous shrinkage
priors. Examples include the Bayesian Lasso (Park & Casella, 2008), the
Horseshoe prior (Carvalho et al., 2010), the Normal-Gamma prior (Brown &
Griffin, 2010), the Dirichlet-Laplace prior (Bhattacharya et al., 2015),
global-local shrinkage priors (Polson & Scott, 2012), the Beta-prime prior
(Bai & Ghosh, 2018), the tail-adaptive shrinkage prior (Lee et al., 2020) and
the Horseshoe-pit prior (Denti et al., 2021). Continuous shrinkage priors tend
to have computational advantages and are very effective in predictive
settings. However, because they place probability zero on any one value of the
parameter space, variable selection can be performed only by thresholding the
posterior distributions of the model coefficients. While ad-hoc techniques
have been devised for this purpose (e.g., see Li & Pati, 2017), thresholding
tends to work well mostly in settings where enough prior information is
available to establish practical significance.

The remainder of the paper is organized as follows: Section 2 reviews the
general definition and properties of PEP priors. Section 3 defines the general
LPEP and conditions under which it is well defined. Section 4 discusses in
detail the special case of the LPEP for logistic regression, including
theoretical properties such as intrinsic and model selection consistency and
computational strategies. Sections 5 and 6 presents empirical results from
simulation studies as well as three different real datasets. Finally, Section
7 discusses our results and future directions for research.

##  2 Power-expected-posterior priors: A brief review

Power-expected-posterior (PEP) priors (Fouskakis et al., 2015) extend the
expected-posterior (EP) priors introduced by Pérez & Berger (2002) by
controlling the amount of information contained in the prior using the prior
power approach originally developed by Ibrahim et al. (2000) and Chen et al.
(2000) in the context of subjective priors based on historical information.

Briefly, let ${\boldsymbol{y}}$ denote the $n$-dimensional vector containing
the observed data, $\boldsymbol{\gamma}$ index the model space, and
$\boldsymbol{\beta}_{\boldsymbol{\gamma}}$ represent vector of parameters
under model $\boldsymbol{\gamma}$. We start with a (potentially improper)
prior
$\pi_{\boldsymbol{\gamma}}^{N}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
under model $\boldsymbol{\gamma}$ and introduce an $n^{*}$-dimensional vector
of imaginary training samples arising from a distribution
$m^{*}(\mathbf{y}^{*})$. The EP prior is then constructed as

$\displaystyle\pi_{\boldsymbol{\gamma}}^{EP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\int\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)}{\int
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}}m^{*}({\boldsymbol{y}}^{*})d{\boldsymbol{y}}^{*}$

In words, the EP priors use the imaginary training sample $\mathbf{y}^{*}$ to
update the original prior
$\pi_{\boldsymbol{\gamma}}^{N}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$, and
addresses the possible effect of using any particular training sample by
averaging over the distribution $m^{*}(\mathbf{y}^{*})$. The use of a common
$m^{*}({\boldsymbol{y}}^{*})$ properly calibrates the priors across the
different values of $\boldsymbol{\gamma}$, even in situations where
$m^{*}({\boldsymbol{y}}^{*})$ is improper. Pérez & Berger (2002) discuss
various possible choices of $m^{*}({\boldsymbol{y}}^{*})$ in both informative
and non-informative settings.

Note that an implicit assumption in the formulation of the PEP is that the
training sample $\mathbf{y}^{*}$ must be large enough so that the posterior
based on it is proper, i.e.,

$\displaystyle\int
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}<\infty$
(1)

for any ${\boldsymbol{y}}^{*}$ in the support of $m^{*}$. However, large the
values of $n^{*}$ will produce priors that are relatively concentrated. To
balance these two goals, it is common to choose $n^{*}$ as the size of the
minimum training sample required to satisfy (1) across all models. This
complicates the implementation of this prior in situations where data is not
independent and identically distributed, such as in regression models.

In summary, even though the EP prior attempts to ameliorate the effect of the
${\boldsymbol{y}}^{*}$ by averaging over $m^{*}$ and by using training samples
that are as small as possible, in some applications the prior might be quite
concentrated, and therefore highly informative. Power-expected-posterior
priors (Fouskakis et al., 2015) address this drawback by scaling the
likelihood of the imaginary sample,

$\displaystyle\pi_{\boldsymbol{\gamma}}^{PEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\int\frac{\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)}{\int\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}}m^{*}({\boldsymbol{y}}^{*}\mid\delta)f(\delta\mid\boldsymbol{\gamma})d\delta
d{\boldsymbol{y}}^{*}$

where
$\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)=\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)^{\frac{1}{\delta}}}{\int
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)^{\frac{1}{\delta}}d\boldsymbol{\beta}_{\boldsymbol{\gamma}}}$
is the normalized power likelihood for the training sample
${\boldsymbol{y}}^{*}$ based on model $\boldsymbol{\gamma}$, and $\delta$ is
the power parameter. If $\delta=1$, then PEP prior reduces to the EP prior,
while values of $\delta>1$ yield priors with a larger variance (and therefore,
less information) than the EP prior. A particularly appealing choice is
$\delta=n^{*}$ (or, alternatively, a prior on $\delta$ that is concentrated
around $n^{*}$), which leads to a prior that can be considered as being unit
information (Kass & Wasserman, 1995). Note that $\delta$ plays a similar role
to the $g$ parameter involved in the definition of (mixtures of) $g$ priors,
so that treating $\delta$ as random will typically lead to priors that have
heavier tails, and are therefore more robust (in the sense of Dawid, 1973 and
Andrade & O’Hagan, 2011).

Being able to use the parameter $\delta$ to control the amount of information
contained in the prior means that the choice of the size of training sample is
less critical in the case of PEP priors. In the sequel, we work with
$n^{*}=n$, a choice that is particularly convenient when dealing with GLMs and
other regression models. Indeed, taking $n^{*}=n$ allows us to select
${\boldsymbol{X}}^{*}$, the design matrix associated with the training sample
${\boldsymbol{y}}^{*}$, as ${\boldsymbol{X}}^{*}={\boldsymbol{X}}$, the design
matrix associated with the observed data. A further implicit assumption moving
forward is that $n>p$.

The PEP prior was originally derived for model selection in Gaussian linear
model. In that case, computing the normalizing constant $\int
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)^{\frac{1}{\delta}}d\boldsymbol{\beta}_{\boldsymbol{\gamma}}$
associated with
$\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)$
is straightforward. Indeed, for most standard choices of
$\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)$,
the induced PEP can be written as a location-and-scale mixture of Gaussian
distributions, dramatically simplifying computation within a Markov chain
Monte Carlo framework. This property, however, does not extend to other GLMs.
To address this issue, Fouskakis et al. (2018) introduce two slightly
different modifications of the PEP framework that rely on the unnormalized
power likelihood
$f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)^{\frac{1}{\delta}}$
rather than
$\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)$:
the concentrated reference PEP (CRPEP) and the diffuse reference PEP (DRPEP).
However, while the use of the unnormalized power likelihood avoids some of the
computational difficulties associated with the original PEP prior, many of
them remain. In particular, neither
$\pi_{\boldsymbol{\gamma}}^{CRPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
nor
$\pi_{\boldsymbol{\gamma}}^{DRPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
belong to standard families of distributions. This prevents closed-form
integration of the regression coefficients and therefore requires the use
Reversible Jump Markov chain Monte Carlo algorithms. Furthermore, the
definition of the CRPEP and the DRPEP and the computational approach
introduced by the authors (which relies on Laplace approximations to compute
certain normalizing constants needed for the acceptance probabilities of
various Metropolis-Hastings steps) implicitly assume that the maximum
likelihood estimate of $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$ exists for
any training sample ${\boldsymbol{y}}^{*}$ and model $\boldsymbol{\gamma}$.
However, this constraint, which as we discussed in the introduction might be
binding for some classes of GLMs, is not accounted for in the definition of
the $m^{*}({\boldsymbol{y}}^{*}\mid\delta)$ CRPEP and DRPEP.

##  3 The Laplace power-expected-posterior prior for Generalized Linear Models

Instead of working with the unnormalized power likelihood as in Fouskakis et
al. (2018), in this paper we propose replacing the likelihood of the imaginary
samples with its Laplace approximation before raising it to the power
$1/\delta$. Hence, the name Laplace PEP, or LPEP. More concretely, let the
observations ${\boldsymbol{y}}=(y_{1},\ldots,y_{n})^{T}$ be generated from a
likelihood of the form

$f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)=\prod_{i=1}^{n}h(y_{i},\tau)\exp\left\\{\frac{T(y_{i})G\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)+A\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)}{d(\tau)}\right\\}$

where
$E(y_{i})=\eta({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma})=\theta_{i}$
for some appropriate link function $\eta$,
${\boldsymbol{x}}_{i}=(1,x_{i,1},\dots,x_{i,p})$ is the $p+1$ dimensional
vector of regressors associated with observation $y_{i}$,
$\boldsymbol{\beta}\in S$ is the $p+1$ dimensional vector of regression
coefficients (including the intercept), $S$ is a connected open subset of
$\mathbb{R}^{p+1}$,
$\boldsymbol{\gamma}^{\prime}=(\gamma_{0},\gamma_{1},\ldots,\gamma_{p})$ is a
binary vector of length $p+1$ such that for all
$j\in\\{1,\dots,p\\},\gamma_{j}=1$ if the $j$-th variable is included in the
model (i.e., if $\beta_{j}$ is different from zero) and $\gamma_{j}=0$
otherwise and $\gamma_{0}=1$ (i.e., intercept is always included in the
model), ${\boldsymbol{x}}_{\boldsymbol{\gamma},i}$ and
$\boldsymbol{\beta}_{\gamma}$ denote the sub-vectors of ${\boldsymbol{x}}_{i}$
and $\boldsymbol{\beta}$ with length $p_{\boldsymbol{\gamma}}+1$ where
$p_{\boldsymbol{\gamma}}=\sum_{j=1}^{p}\gamma_{j}$ that include only those
components for which the corresponding $\gamma_{j}$ is equal to 1, and $\tau$
is an overdispersion parameter. In order to simplify our exposition, and
following standard practice in the literature, in the sequel we treat $\tau$
as known. The generalization of our approach to situations where $\tau$ is
unknown is relatively straightforward as $\tau$ will typically be a parameter
that is common to all models under consideration, and can therefore be safely
assigned a standard (potentially improper) non-informative prior (Berger et
al., 1998). When $\tau$ is known, the normalizing constant $h(y_{i},\tau)$ can
be dropped, and the scaling function $d(\tau)$ can be absorbed within the
function $T$ and $A$, leading to the somewhat simpler expression:

$\displaystyle
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)$
$\displaystyle=\exp\left\\{\sum_{i=1}^{n}\left[T^{*}(y_{i})G\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)+A^{*}\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)\right]\right\\},$
(2)

where $A^{*}(\cdot)=A(\cdot)/d(\tau)$ and $T^{*}(\cdot)=T(\cdot)/d(\tau)$.

The second order Laplace approximation to (2) is given by

$\displaystyle
f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)$
$\displaystyle\approx
f^{L}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)\propto\exp\left\\{-\frac{1}{2}\left(\boldsymbol{\beta}_{\gamma}-\hat{\boldsymbol{\beta}}_{\gamma}({\boldsymbol{y}})\right)^{\prime}{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}})\left(\boldsymbol{\beta}_{\gamma}-\hat{\boldsymbol{\beta}}_{\gamma}({\boldsymbol{y}})\right)\right\\},$
(3)

where
$\tilde{\ell}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\sum_{i=1}^{n}\left[T^{*}(y_{i})G\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)+A^{*}\left(\eta\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}\boldsymbol{\beta}_{\gamma}\right)\right)\right]$,
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}})$ denotes the
maximum likelihood estimate for $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$
based on sample ${\boldsymbol{y}}$, and
${\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}})$ is the
$(p_{\boldsymbol{\gamma}}+1)\times(p_{\boldsymbol{\gamma}}+1)$ observed
information matrix with entries

$[{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}})]_{j,j^{\prime}}=-\left.\frac{\partial}{\partial\beta_{\boldsymbol{\gamma},j}\partial\beta_{\boldsymbol{\gamma},j^{\prime}}}\tilde{\ell}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})\right|_{\boldsymbol{\beta}_{\boldsymbol{\gamma}}=\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}(\mathbf{y})}.$

In the case of regular exponential families, it is well known that this
approximation is accurate up to an $\mathcal{O}(\frac{1}{n})$ order term
(e.g., see Schwarz, 1978 and Haughton, 1988). With this in mind, we define the
LPEP as

$\displaystyle\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\int\frac{\tilde{f}^{L}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)}{\int\tilde{f}^{L}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}}m^{*}({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}})f(\delta\mid\boldsymbol{\gamma})d\delta
d{\boldsymbol{y}}^{*},$ (4)

where ${\boldsymbol{X}}$ is the $n\times(p+1)$ design matrix whose rows
correspond to the ${\boldsymbol{x}}_{i}^{\prime}$ vectors.

Defining the LPEP using (3) instead of (2) dramatically simplifies
computation. Indeed, using (3) implies that
$\tilde{f}^{L}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)$
is proportional to a Gaussian kernel,

$\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\propto\delta^{-\frac{p_{\gamma}+1}{2}}\exp\left\\{-\frac{1}{2\delta}\left(\boldsymbol{\beta}_{\gamma}-\hat{\boldsymbol{\beta}}_{\gamma}({\boldsymbol{y}}^{*})\right)^{\prime}{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\left(\boldsymbol{\beta}_{\gamma}-\hat{\boldsymbol{\beta}}_{\gamma}({\boldsymbol{y}}^{*})\right)\right\\},$

and, therefore, for standard choices of
$\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)$
(such as the flat prior
$\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)\propto
1$), $\pi^{LPEP}$ corresponds to a location-and-scale mixture of Gaussian
distributions.

Besides using the the Laplace approximation of
$f({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})$, the
definition in (4) differs from that in Fouskakis et al. (2018) in terms of the
structure of the distribution of the imaginary samples $m^{*}$ in two
important ways. First, note that we do not make $m^{*}$ dependent of the
scaling factor $\delta$. This makes intuitive sense (there is no obvious
reason why the power factor used to re-scale the information in the training
sample should also affect how the training sample is generated) and simplifies
both posterior computation and theoretical analysis. On the other hand,
$m^{*}$ is made to depend explicitly on the design matrix ${\boldsymbol{X}}$.
In particular, we take

$\displaystyle m^{*}({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}})$
$\displaystyle\propto\tilde{m}^{*}({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}})\mathbf{1}\left(\mathbf{y}^{*}\in
A({\boldsymbol{X}})\right),$ (5)

where $\mathbf{1}\left(\mathbf{y}^{*}\in A({\boldsymbol{X}})\right)$ is the
indicator function on the set $A({\boldsymbol{X}})$ and

$\displaystyle A({\boldsymbol{X}})$
$\displaystyle=\left\\{\tilde{{\boldsymbol{y}}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\left(\tilde{{\boldsymbol{y}}}\right)\mbox{
exists and is finite for all }\boldsymbol{\gamma}\right\\}.$ (6)

At first sight, the computational implementation of (6) might seem daunting,
as it in principle requires that the existence of the maximum likelihood
estimator of the parameters be checked for every possible model under
consideration. However, as the following theorem shows, for a broad class of
GLMs that includes binomial and multinomial models with logistic and probit
links, as well as Poisson models with the canonical logarithmic link, it is
enough to check relatively simple conditions on the loglikelihood funciton for
the full model.

######  Theorem 3.1.

Let $\ell_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
denote the log-likelihood function of the GLM in (2) associated with model
$\boldsymbol{\gamma}$ and $\boldsymbol{\gamma}_{F}=(1,1,\ldots,1)$ denote the
full model (i.e., the model that includes all potential regressors). Assume
that

  1. (i)

$\ell_{\boldsymbol{\gamma}_{F}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}})$
is continuous and striclty concave on $S_{\boldsymbol{\gamma}_{F}}$

  2. (ii)

$\lim_{\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}}\to\partial
S_{\boldsymbol{\gamma}_{F}}}\ell_{\boldsymbol{\gamma}_{F}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}})=-\infty$,
where $\partial S_{\boldsymbol{\gamma}_{F}}$ represents the closure of
$S_{\boldsymbol{\gamma}_{F}}$

Then, $\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}})$, the
maximum likelihood estimator, exists under any other model
$\boldsymbol{\gamma}$.

The proof of Theorem 3.1, which is relatively straightforward, can be seen in
Appendix A. Similar theorems covering wider classes of GLMs can be derived
using the results in Wedderburn (1976), but this version suffices for our
purposes.

To conclude, it might be helpful to expand on the relationship between the
LPEP prior as we defined it in this Section and other similar proposals in the
literature within the context of Gaussian linear models. Note that, while the
Laplace approximation in (3) is exact for the Gaussian linear model, our
definition of the LPEP is not equivalent to the PEP prior in Fouskakis et al.
(2015) because the predictive distribution $m^{*}({\boldsymbol{y}}^{*})$ used
the generate the training samples does not depend on $\delta$. On the other
hand, the LPEP for the Gaussian linear model looks similar to a mixture of $g$
priors Liang et al. (2008). In particular, while the covariance matrix of the
Gaussian kernels involved is given by
$\delta\left\\{{\boldsymbol{X}}^{T}{\boldsymbol{X}}\right\\}^{-1}$ in both
priors, the conditional mean of the LPEP depends on the training sample and is
not zero. However, depending on the choice of $m^{*}(\mathbf{y}^{*})$, the
intrinsic prior associated with the LPEP does correspond to a mixture of $g$
priors (please see Section 4.1.3).

##  4 The LPEP prior for logistic regression

In the sequel, we illustrate the LPEP in the context of logistic regression,
where the likelihood function for model $\boldsymbol{\gamma}$ can be written
as

$\displaystyle
f_{\boldsymbol{\gamma}}(\boldsymbol{y}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\exp\left\\{\sum_{i=1}^{n}\left[\boldsymbol{x}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}-\log(1+\exp\left\\{\boldsymbol{x}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right\\})\right]\right\\},$
(7)

and $S_{\boldsymbol{\gamma}}=\mathbb{R}^{p_{\boldsymbol{\gamma}}+1}$. A
natural choice for
$\pi^{N}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ in
this setting is the (improper) flat prior
$\pi^{N}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})\propto
1$, so that

$\displaystyle\pi^{LPEP}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\sum_{{\boldsymbol{y}}^{*}\in\\{0,1\\}^{n}}\Bigg{[}\int\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}),\delta{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\right)\right)f(\delta\mid\boldsymbol{\gamma})d\delta\Bigg{]}m^{*}({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}}),$
(8)

where $\phi_{p}\left(\cdot\mid\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$
denotes the density of the $p$-variate normal distribution with mean vector
$\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$,
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})$ is the
maximum likelihood estimator of $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$
based on the training sample ${\boldsymbol{y}}^{*}$ (which, while not
available in closed form, can be easily evaluated),

${\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})=\sum_{i=1}^{n}\frac{1}{\hat{\theta}_{\boldsymbol{\gamma},i}({\boldsymbol{y}}^{*})(1-\hat{\theta}_{\boldsymbol{\gamma},i}({\boldsymbol{y}}^{*}))}{\boldsymbol{x}}_{\boldsymbol{\gamma},i}{\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{\prime}={\boldsymbol{X}}_{\boldsymbol{\gamma}}^{\prime}{\boldsymbol{W}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}){\boldsymbol{X}}_{\boldsymbol{\gamma}},$

where

$\hat{\theta}_{\boldsymbol{\gamma},i}({\boldsymbol{y}}^{*})=\frac{1}{1+\exp\left\\{\boldsymbol{x}_{\boldsymbol{\gamma},i}^{T}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right\\}},$

and

${\boldsymbol{W}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})=\operatorname{diag}\left\\{\frac{1}{\hat{\theta}_{\boldsymbol{\gamma},1}({\boldsymbol{y}}^{*})(1-\hat{\theta}_{\boldsymbol{\gamma},1}({\boldsymbol{y}}^{*}))},\ldots,\frac{1}{\hat{\theta}_{\boldsymbol{\gamma},n}({\boldsymbol{y}}^{*})(1-\hat{\theta}_{\boldsymbol{\gamma},n}({\boldsymbol{y}}^{*}))}\right\\}.$

Note that the loglikelihood function
$\ell_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\sum_{i=1}^{n}\left[\boldsymbol{x}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}-\log\left(1+\exp\left\\{\boldsymbol{x}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right\\}\right)\right]$
associated with (7) is continuous everywhere for any model
$\boldsymbol{\gamma}$ (and, in particular, for $\boldsymbol{\gamma}_{F}$).
Furthermore, as long as the full design matrix ${\boldsymbol{X}}$ is full
rank, ${\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})$ is
strictly positive definite for any $\boldsymbol{\gamma}$. Hence, the model
satisfies condition (i) in Theorem 3.1. To verify that condition (ii) is
satisfied, it is enough to show that the training sample
${\boldsymbol{y}}^{*}$ is not separable under the full design matrix
${\boldsymbol{X}}$. Konis (2007) discusses an efficient approach to detect
separation in logistic regression models that relies on linear programming.
This approach has been implemented in the R package detectseparation (Kosmidis
& Schumacher, 2020).

We discuss next the choice of $\tilde{m}^{*}({\boldsymbol{y}}^{*})$ in (5). A
common approach is to select $m^{*}({\boldsymbol{y}}^{*})$ as the predictive
under the the null model. In that spirit, for the logistic regression model
(and, more generally, for any binary regression model, independently of the
link function used) we set

$\displaystyle\tilde{m}^{*}({\boldsymbol{y}}^{*})$
$\displaystyle=\frac{\Gamma\left(\sum_{i=1}^{n}y_{i}^{*}+\frac{1}{2}\right)\Gamma\left(n-\sum_{i=1}^{n}y_{i}^{*}+\frac{1}{2}\right)}{\Gamma(n+1)\Gamma\left(\frac{1}{2}\right)\Gamma\left(\frac{1}{2}\right)}$
$\displaystyle=\int
t^{\sum_{i=1}^{n}y_{i}^{*}}(1-t)^{n-\sum_{i=1}^{n}y_{i}^{*}}\,\frac{t^{-\frac{1}{2}}(1-t)^{-\frac{1}{2}}}{\Gamma\left(\frac{1}{2}\right)\Gamma\left(\frac{1}{2}\right)}dt,$

a Beta-Binomial distribution with both parameters equal to $\frac{1}{2}$. This
choice is particularly appealing because it corresponds to the predictive
distribution under the null model and its reference/Jeffreys prior.

To complete the specification of the LPEP we must specify the mixing
distribution for the exponent $\delta$. In this manuscript we consider three
alternatives. The first version of the LPEP we investigate is the unit
information LPEP (UI-LPEP) obtained by fixing $\delta=n^{*}$. We also consider
a version of the hyper-g/n prior discussed in Liang et al. (2008) and Li &
Clyde (2018),

$\displaystyle f^{HGN}(\delta)$
$\displaystyle=\left(1+\frac{\delta}{n^{*}}\right)^{-2}.$

We call this the HGN-LPEP The median of the hyper-g/n prior is equal to
$n^{*}$, and the prior places much of its mass around this value. It can
therefore be considered as a relaxation of the unit information version of the
prior. Note that, under this hyperprior, $\delta$ is independent of the model
under consideration.

Finally, we consider a version of the robust prior recommended by Bayarri et
al. (2012),

$\displaystyle f^{R}(\delta|\boldsymbol{\gamma})$
$\displaystyle=\frac{1}{2(p_{\boldsymbol{\gamma}}+1)^{1/2}}\frac{(n^{*}+1)^{1/2}}{(\delta+1)^{3/2}}\mathbf{1}\left(\delta>\frac{n^{*}-p_{\boldsymbol{\gamma}}}{p_{\boldsymbol{\gamma}}+1}\right),$

which we call the R-LPEP The robust prior satisfies a number of compelling
desiderata for Gaussian linear regression models. For example, its expectation
is $\mathcal{O}(n^{*})$. However, note that, unlike our previous two choices,
this prior depends on the model size.

###  4.1 Properties of the LPEP prior for logistic regression

####  4.1.1 Proper prior

Note that our choice of
$m^{*}\left({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}}\right)$ is proper, and
that because of the constraint on the value of the training samples,

$\pi_{\boldsymbol{\gamma}}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*},\delta\right)=\frac{\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)}{\int\tilde{f}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},{\delta}\right)\pi_{\boldsymbol{\gamma}}^{N}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}}$

is also proper for every $\boldsymbol{\gamma}$. Therefore, the LPEP prior in
(8) is also proper for every $\boldsymbol{\gamma}$.

####  4.1.2 Tail behavior

It straightforward to see that the unit information version of the LPEP (where
$\delta=n^{*}$), $\pi^{LPEP-
UI}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ has
Gaussian tails. On the other hand, as the following theorem shows, the
hyper-g/n and the robust versions of the LPEP have heavier (polynomial) tails
in every direction.

######  Theorem 4.1.

Let

$\zeta^{HGN}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma})=\left.\pi_{\boldsymbol{\gamma}}^{LPEP-
HGN}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})\right|_{\boldsymbol{\beta}_{\boldsymbol{\gamma}}=s{\boldsymbol{v}}}$

and

$\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma})=\left.\pi_{\boldsymbol{\gamma}}^{LPEP-R}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})\right|_{\boldsymbol{\beta}_{\boldsymbol{\gamma}}=s{\boldsymbol{v}}}$

for any vector ${\boldsymbol{v}}$ such that $\|{\boldsymbol{v}}\|=1$. Then
there exist functions $c^{HGN}_{\boldsymbol{\gamma}}({\boldsymbol{v}})$ and
$c^{R}_{\boldsymbol{\gamma}}({\boldsymbol{v}})$ such that

$\lim_{s\to\infty}\frac{\zeta^{HGN}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=c^{HGN}_{\boldsymbol{\gamma}}({\boldsymbol{v}})<\infty$

and

$\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=c^{R}_{\boldsymbol{\gamma}}({\boldsymbol{v}})<\infty$

for every direction ${\boldsymbol{v}}$ and model $\boldsymbol{\gamma}$.

The proof, which is presented in Appendix B, extends results originally
presented in Bayarri et al. (2012).

One important consequence of this result is that, from an estimation (rather
than model selection) perspective, $\pi^{LPEP-
HGN}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ and
$\pi^{LPEP-R}_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
are robust, in the sense of having bounded influence in the case of
likelihood-prior conflict (e.g., see Andrade & O’Hagan, 2006 and Andrade &
O’Hagan, 2011).

####  4.1.3 Intrinsic consistency

In addition to being proper, under mild conditions the LPEP converges to a
non-degenerate prior as the size of the training sample increases. Naturally,
the exact form of the intrinsic prior depends on the asymptotic regime for the
covariates associated with new observations, as well as the exact prior used
for $\delta$. Theorem 4.2 below provides a relevant example.

######  Theorem 4.2.

Assume that, as $n$ (and therefore, $n^{*}$) grows, the covariate vectors
${\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\ldots$ satisfy either of the
following two conditions:

  * (i)

If ${\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\ldots$ forms a deterministic
sequence, then
$\frac{1}{n}{\boldsymbol{X}}^{T}{\boldsymbol{X}}\underset{n\to\infty}{\longrightarrow}\boldsymbol{\Sigma}$.

  * (ii)

If ${\boldsymbol{x}}_{1},{\boldsymbol{x}}_{2},\ldots$ are random, then they
are independent and identically distributed from a distribution with mean
$\mathbf{0}$ and finite covariance matrix $\boldsymbol{\Sigma}$.

Then, the unit information ($\delta=n^{*}$), hyper-g/n and robust versions of
the LPEP have proper, non-degenerate intrinsic priors of the form

$\displaystyle\lim_{n^{*}\to\infty}\pi_{\boldsymbol{\gamma}}^{LPEP-
UI}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\mathbf{0},4\left[\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}\right]^{-1}\right),$
$\displaystyle\lim_{n^{*}\to\infty}\pi_{\boldsymbol{\gamma}}^{LPEP-
HGN}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\int\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\mathbf{0},4\delta^{*}\left[\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}\right]^{-1}\right)\left(1+\delta^{*}\right)^{-2}d\delta^{*},$

and

$\displaystyle\lim_{n^{*}\to\infty}\pi_{\boldsymbol{\gamma}}^{LPEP-R}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=\int\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\mathbf{0},4\delta^{*}\left[\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}\right]^{-1}\right)\frac{1}{2(p_{\boldsymbol{\gamma}}+1)^{\frac{1}{2}}}\left(\frac{1}{\delta^{*}}\right)^{\frac{3}{2}}\mathbf{1}\left(\delta^{*}>\frac{1}{p_{\boldsymbol{\gamma}}+1}\right)d\delta^{*},$

where $\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}$ denotes the square submatrix
of $\boldsymbol{\Sigma}$ that includes only the rows and columns for which
$\gamma_{j}=1$.

A proof of this result (versions of which have been discussed in Li & Clyde,
2018) can be seen in Appendix C. Interestingly, we note that these are the
same intrinsic priors associated with the prior in Bové et al. (2011) under
the same asymptotic regime for ${\boldsymbol{X}}$.

####  4.1.4 Model selection consistency

Model selection consistency refers to the ability of the procedure to choose
the correct model as $n\to\infty$. Intuitively, because the amount of
information in $\pi_{\boldsymbol{\gamma}}^{LPEP}$ is kept approximately
constant as $n^{*}$ increases, we would expect that the associated Bayes
factors would behave asymptotically like those computed from the Bayesian
Information Criteria, which have been well studied and are known to be
consistent. The following theorem, a proof of which can be seen in Appendix D,
formalizes that intuition for the unit information prior.

######  Theorem 4.3.

Assume that a sequence of observations $y_{1},y_{2},\ldots$ is generated from
some model $\boldsymbol{\gamma}_{T}\in\\{1\\}\times\\{0,1\\}^{p}$ (i.e., one
of the models considered by our procedure), and that $p$ is fixed. For the
unit information LPEP with $g=n$, and under appropriate regularity conditions
for how the covariate vectors
${\boldsymbol{x}}_{\boldsymbol{\gamma},1},{\boldsymbol{x}}_{\boldsymbol{\gamma},2},\ldots$
are generated, we have

$\lim_{n\to\infty}\mathsf{Pr}^{LPEP}(\boldsymbol{\gamma}=\boldsymbol{\gamma}_{T}\mid{\boldsymbol{y}})=1.$

####  4.1.5 Information consistency

Information consistency refers to the behavior of the model selection criteria
for fixed sample $n$, as the observed sample ${\boldsymbol{y}}$ becomes
increasingly more “extreme”. In the case of logistic regression, because the
sample space for ${\boldsymbol{y}}$ is finite for any $n$, traditional issues
of information inconsistency do not arise.

###  4.2 Markov chain Monte Carlo sampling

The LPEP prior can be easily combined with the Polya-Gamma augmentation of
Polson et al. (2013) to generate an efficient Markov chain Monte Carlo
algorithm for variable selection in logistic regression. For this purpose it
is convenient to re-express (4) as a hierarchical prior where

$\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*},\delta,\boldsymbol{\gamma}\sim\mathsf{N}\left(\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}),\delta\left\\{{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right\\}^{-1}\right),$

with $y^{*}\sim m^{*}\left(y^{*}\mid{\boldsymbol{X}}\right)$ and $\delta\sim
f(\delta|\boldsymbol{\gamma})$.

Now, from Theorem 3.1 of Polson et al. (2013), we can write

$\displaystyle
f_{\boldsymbol{\gamma}}(\boldsymbol{y}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
$\displaystyle=\prod_{i=1}^{n}\frac{\exp\left\\{y_{i}{\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right\\}}{1+\exp{(\boldsymbol{x}_{i,\boldsymbol{\gamma}}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}})}}$
$\displaystyle\propto\prod_{i=1}^{n}\left(\exp\left\\{(y_{i}-1/2){\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right\\}\int_{0}^{\infty}\exp\left\\{-\frac{\omega_{i}}{2}\left({\boldsymbol{x}}_{\boldsymbol{\gamma},i}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}}\right)^{2}\right\\}f(\omega_{i}\mid
1,0)d\omega_{i}\right),$

where $f(\omega\mid a,b)$ denotes the density of a Pòlya-Gamma random variate
with parameters $a$ and $b$. Therefore, after introducing a vector of
auxiliary random variables
$\boldsymbol{\omega}=(\omega_{1},\ldots,\omega_{n})$,

$\displaystyle
f(\boldsymbol{\gamma},\boldsymbol{\beta}_{\boldsymbol{\gamma}},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})$
$\displaystyle\propto
f(\boldsymbol{\gamma})f(\delta|\boldsymbol{\gamma})\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}),\delta{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right)\phi_{n}\left({\boldsymbol{z}}\mid{\boldsymbol{X}}_{\boldsymbol{\gamma}}\boldsymbol{\beta}_{\boldsymbol{\gamma}},\boldsymbol{\Omega}^{-1}\right),$
(9)

where
${\boldsymbol{z}}=((y_{1}-1/2)/\omega_{1},\ldots,(y_{n}-1/2)/\omega_{n})^{T}$,
$\boldsymbol{\Omega}=\operatorname{diag}\left\\{\omega_{1},\ldots,\omega_{n}\right\\}$,
$f(\boldsymbol{\gamma})$ is a prior on $2^{p}$ dimensional model space and
$f(\delta|\boldsymbol{\gamma})$ is the prior on scale parameter $\delta$.

It is straightforward to see that $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$
can be integrated out of (9), yielding

$\displaystyle
f(\boldsymbol{\gamma},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})=\int
f(\boldsymbol{\gamma},\boldsymbol{\beta}_{\boldsymbol{\gamma}},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})d\boldsymbol{\beta}_{\boldsymbol{\gamma}}\propto
f(\boldsymbol{\gamma})f(\delta|\boldsymbol{\gamma})\phi_{n}({\boldsymbol{z}}|{\boldsymbol{m}}_{{\boldsymbol{z}}}^{\boldsymbol{\gamma}},{\boldsymbol{V}}_{{\boldsymbol{z}}}^{\boldsymbol{\gamma}}),$
(10)

where
${\boldsymbol{m}}_{{\boldsymbol{z}}}^{\boldsymbol{\gamma}}={\boldsymbol{X}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})$,
${\boldsymbol{V}}_{{\boldsymbol{z}}}^{\boldsymbol{\gamma}}=\boldsymbol{\Omega}^{-1}+\delta{\boldsymbol{X}}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}){\boldsymbol{X}}_{\boldsymbol{\gamma}}^{T}$
and, as before, $\phi_{p}(\cdot\mid\boldsymbol{\mu},\boldsymbol{\Sigma})$
denotes the density of the $p$-variate normal distribution with mean vector
$\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$. Various
versions of Metropolis-Hastings algorithms can be implemented to explore the
space of model (e.g., see section 4.5 of George & McCulloch, 1997).

Once the model $\boldsymbol{\gamma}$ and the exponent $\delta$ have been
updated, the regression coefficients can be sampled using the fact that
$\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\boldsymbol{\gamma},\delta,{\boldsymbol{y}}^{*},\boldsymbol{\omega}\sim\mathsf{N}\left({\boldsymbol{m}}_{\boldsymbol{\omega}},{\boldsymbol{V}}_{\boldsymbol{\omega}}\right)$,
where

$\displaystyle{\boldsymbol{m}}_{\boldsymbol{\omega}}$
$\displaystyle={\boldsymbol{V}}_{\boldsymbol{\omega}}\left({\boldsymbol{X}}_{\boldsymbol{\gamma}}\boldsymbol{\Omega}{\boldsymbol{z}}+\frac{1}{\delta}{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right),$
$\displaystyle{\boldsymbol{V}}_{\boldsymbol{\omega}}$
$\displaystyle=\left({\boldsymbol{X}}_{\boldsymbol{\gamma}}^{T}\boldsymbol{\Omega}{\boldsymbol{X}}_{\boldsymbol{\gamma}}+\frac{1}{\delta}{\boldsymbol{H}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right)^{-1}$
(11)

Conditional on $\boldsymbol{\beta}_{\boldsymbol{\gamma}}$, the remaining
parameters ${\boldsymbol{y}}^{*}$ and $\boldsymbol{\omega}$ can be easily
sampled using either Gibbs sampling or random-walk Metropolis-Hastings steps.
Further details of the computational algorithm can be see in Appendix E.

##  5 Simulation studies

We conducted two simulation studies to compare the estimation and model
selection performance of Laplace PEP priors with other existing model
selection techniques. The setup for our simulation study is motivated by that
in Li & Clyde (2018). This section discusses the results from the first
simulation study, the results for the second one can be seen in the
supplementary materials.

The simulation study described in this section uses a sample size of $n=500$
and a total number of covariates $p=p_{\boldsymbol{\gamma}_{F}}=100$, with the
vectors of predictors being drawn independently from a zero-mean, unit-scale
multivariate normal distribution with pairwise correlations given by
$cor(x_{i,j},x_{i,j^{\prime}})=r^{|j-j^{\prime}|}$ for $1\leq j<j^{\prime}\leq
p$. It consists of eight scenarios, which differ in terms of both the sparsity
level in the vector of regression coefficients and the the correlation
structure among predictors. More specifically, we consider all combinations of
four different levels of sparsity
($p_{\boldsymbol{\gamma}_{T}}\in\\{0,5,10,20\\}$, please see Table 1) and two
different correlation coefficients ($r\in\\{0,0.75\\}$).

$p_{\boldsymbol{\gamma}_{T}}$ | $\beta_{\boldsymbol{\gamma}_{T},0}$ | $\beta_{\boldsymbol{\gamma}_{T},1:5}$ | $\beta_{\boldsymbol{\gamma}_{T},6:10}$ | $\beta_{\boldsymbol{\gamma}_{T},11:15}$ | $\beta_{\boldsymbol{\gamma}_{T},16:20}$
---|---|---|---|---|---
0 | $-0.5$ | ${\boldsymbol{0}}$ | ${\boldsymbol{0}}$ | ${\boldsymbol{0}}$ | ${\boldsymbol{0}}$
5 | $-0.5$ | ${\boldsymbol{b}}$ | ${\boldsymbol{0}}$ | ${\boldsymbol{0}}$ | ${\boldsymbol{0}}$
10 | $-0.5$ | ${\boldsymbol{b}}$ | ${\boldsymbol{0}}$ | ${\boldsymbol{b}}$ | ${\boldsymbol{0}}$
20 | $-0.5$ | ${\boldsymbol{b}}$ | $0.5{\boldsymbol{b}}$ | ${\boldsymbol{b}}$ | $0.5{\boldsymbol{b}}$

Table 1: Value of intercept and coefficients in the true logistic regression
model where ${\boldsymbol{b}}=(2,-1,-1,0.5,-0.5)^{T}$

.

A total of 100 datasets were generated for each of our 8 scenarios. We apply
both Bayesian procedures and various penalized likelihood approaches to each
dataset. In terms of Bayesian procedures, in addition to the LPEP prior, we
also consider the methodology of (Li & Clyde, 2018), which relies of a mixture
of g-priors along with a Laplace approximation to compute the associated
marginal likelihood (denoted LCL in the sequel), as well as an “exact” version
of their procedure that relies on the same mixture of g-priors but avoids the
Laplace approximation by implementing a latent-variable augmentation similar
to the one described in Section 4.2 (denoted LCE in the sequel). Comparing LCL
and LCE allows us to disentangle the effect of the Laplace approximation from
that of the prior choice on the performance of these techniques. For each of
these three approaches, we consider three different settings for the
hyperparameter $\delta$: the unit information prior with $\delta=n$, the
hyper-g/n and robust priors (recall Section 4). We use the R package BAS
(Clyde, 2020) to implement LCL, and a slight modification of our own code to
implement LCE. In all cases we assume a Beta-Binomial(1,1) prior over the
model space, and run the MCMC chain for $2^{17}\approx 131,000$ iterations
after a burn-in of $10,000$ iterations. Note that we do not include the CRPEP
and DRPEP priors from Fouskakis et al. (2018) in this simulation study. We do
this for two main reasons. First, the computational complexity of the code
provided by the authors makes a simulation study like this prohibitive. Not
only is each iteration of the algorithm much more expensive than those of the
other approaches, but the algorithm mixes much more slowly, which means that a
much larger number of iterations are required to get accurate results.
Secondly, and more important, the algorithm we had access to broke down for a
number of our simulated datasets. We incorporate comparisons with CRPEP and
DRPEP in our real data examples in Section 6. In terms of penalized likelihood
methods, we compare against LASSO (Tibshirani, 1996), smoothly clipped
absolute deviation (SCAD) (Fan & Li, 2001) and minimax concave penalty (MCP)
(Zhang, 2010). We use the R package glmnet (Friedman et al., 2010) to
implement LASSO, and the package ncvreg (Breheny & Huang, 2011) for SCAD and
MCP.

We first evaluate the performance of these various methods in terms of model
selection performance using three metrics. First, we report the frequency
(over the 100 datasets) with which the MAP model matches the true model
$\boldsymbol{\gamma}_{T}$ (see Table 2). For the penalized likelihood
approaches (for which a single model is reported for each dataset) the
equivalent metric is simply the number of datasets for which the technique
reported the correct model. This metric was used in both Li & Clyde (2018) and
Fouskakis et al. (2018).

$p$ | 100
---|---
$p(\boldsymbol{\gamma})$ | Beta-Binomial(1,1)
$p_{\boldsymbol{\gamma}_{T}}$ | 0 | 5 | 10 | 20
$r$ | 0 | 0.75 | 0 | 0.75 | 0 | 0.75 | 0 | 0.75
$\delta=n$ | LPEP | 99 |  100* | 45 | 4 |  18* | 0 | 0 | 0
LCE |  100* |  100* | 45 | 5 | 8 | 0 | 0 | 0
LCL |  100* |  100* | 46 | 4 | 11 | 0 | 0 | 0
$\delta\sim\text{robust}$ | LPEP | 99 |  100* |  53* |  6* | 15 | 0 | 0 | 0
LCE | 99 |  100* | 45 |  6* | 0 | 0 | 0 | 0
LCL |  100* |  100* | 46 |  6* | 2 | 0 | 0 | 0
$\delta\sim\text{hyper g/n}$ | LPEP | 98 |  100* | 50 | 5 | 17 | 0 | 0 | 0
LCE | 97 | 99 | 25 | 4 | 0 | 0 | 0 | 0
LCL | 65 | 78 | 3 | 0 | 0 | 0 | 0 | 0
| LASSO | 59 | 65 | 0 | 0 | 0 | 0 | 0 | 0
| SCAD | 57 | 59 | 0 | 0 | 0 | 0 | 0 | 0
| MCP | 73 | 66 | 8 | 0 | 3 | 0 | 0 | 0

Table 2: Number of times (over 100 replications) that the MAP model coincides
with the true model in the logistic regression ; BOLD represent group maximum;
* represent overall maximum.

First note that Bayesian methods tend to clearly outperform penalized
likelihood approaches, in some cases quite dramatically. Focusing now on the
Bayesian approaches, we observe that most of them perform very well when the
data is generated from the null model. This is true both for uncorrelated and
highly correlated covariates. The main exception is LCL under the hyper-g/n
hyperprior, where the MAP algorithm matches the true model in only 65 ($r=0$)
and 78 ($r=0.75$) of the datasets. On the other hand, as the number of non-
zero coefficients in the true model increases, we observe that all approaches
struggle to identify the true model, particularly when the covariates are
highly correlated. In particular, when $p_{\boldsymbol{\gamma}_{T}}=20$, none
of the procedures is able to identify the true model. Nonetheless, it appears
that, overall, LPEP (and, specially, the robust and the hyper-g/n versions of
LPEP) perform the best.

While the MAP metric we discussed above provides some insights into model
performance, it tends to be less informative when there is substantial
uncertainty on the posterior distribution over the model space. Therefore, we
also compute for each dataset the $F_{1}$ score for the MAP (Bayesian
procedures) or selected (penalized likelihood procedures) model, see Figure 1.
In this setting, the $F_{1}$ score is defined as the harmonic mean of
proportion of true positives among “selected” covariates (the precision) and
the proportion of “selected” covariates among true positive covariates (the
recall). The $F_{1}$ score ranges between 0 and 1, with a higher value
indicating better model selection performance. Note that results are not
presented for the null model since the $F_{1}$ score is not well defined in
that case.

Figure 1: F1 score for the MAP model estimated by various methods and prior
combinations for 100 simulated datasets ($n=500,p=100$) under different
scenarios of correlation ($r=0$: left; $r=0.75$: right) and true number of
non-zero coefficients specified in rows
($p-true=p_{\boldsymbol{\gamma}_{T}}$); Red dots represent the average F1
score across 100 simulated datasets.

The $F_{1}$ score provides a much more informative picture of the performance
of these models. In all cases, the methods based on LPEP priors tend to
perform the best, with the robust and hyper-g/n versions being slightly better
than that of the unit information prior. We also see that, while all Bayesian
procedures have very similar performance under the unit information prior, LCE
tends to outperform LCL under the robust and hyper-g/n priors (in some cases,
quite dramatically).

Next, we also report the average size of the sampled/selected models for each
data set (see Figure 2). Under the robust and hyper-g/n priors, LCE and,
especially, LCL tend to favor very large models. Interestingly, all Bayesian
procedures under the unit information prior seem to underestimate the model
size when $p_{\boldsymbol{\gamma}_{T}}=20$. The best performing approaches
with average model size close to true model size are again the robust and
hyper-g/n versions of LPEP.

Figure 2: Average size of models selected by various methods and prior
combinations for 100 simulated datasets ($n=500,p=100$) under different
scenarios of correlation ($r=0$: left; $r=0.75$: right) and true number of
non-zero coefficients specified in rows
($p-true=p_{\boldsymbol{\gamma}_{T}}$); Dotted blue line indicates the true
model size $p-true=p_{\boldsymbol{\gamma}_{T}}$ and red dots represent the
average model size over 100 simulated datasets.

Finally, we compare the procedures in terms of parameter estimation
performance using the average mean squared error (AMSE) of the estimated
coefficients

$AMSE(\boldsymbol{\beta})=\frac{1}{p}\sum_{j=1}^{p}(\hat{\beta}_{j}-\beta_{j,\boldsymbol{\gamma}_{T}})^{2}$

where $\hat{\beta}_{j}$ and $\beta_{j,\boldsymbol{\gamma}_{T}}$ are the
estimated and true values of $j^{th}$ covariate, respectively. For the
Bayesian procedures, model-averaged posterior mean estimates are used for this
calculation. For penalized likelihood methods, the sparse point estimates of
the coefficients are used. The results can be seen in Table 3.

$p$ | 100
---|---
$p(\boldsymbol{\gamma})$ | Beta-Binomial(1,1)
$p_{\boldsymbol{\gamma}_{T}}$ | 0 | 5 | 10 | 20
$r$ | 0 | 0.75 | 0 | 0.75 | 0 | 0.75 | 0 | 0.75
$\delta=n$ | LPEP | 0.11 |  0.10* | 2.91 | 7.67 | 7.09 | 17.67 | 14.70 | 33.90
LCE | 0.11 |  0.10* | 3.06 | 7.78 | 7.64 | 18.44 | 16.11 | 36.47
LCL |  0.10* |  0.10* | 2.87 | 7.68 | 6.78 | 18.17 | 16.22 | 36.43
$\delta\sim\text{robust}$ | LPEP | 0.12 |  0.10* |  2.62* |  6.87* |  6.04* |  14.07* | 13.38 |  24.03*
LCE | 0.12 | 0.11 | 4.83 | 7.80 | 47.30 | 23.30 | 96.14 | 52.80
LCL |  0.10* |  0.10* | 8.86 | 8.44 | 214.63 | 60.56 | 275.93 | 115.58
$\delta\sim\text{hyper g/n}$ | LPEP | 0.16 | 0.14 | 2.70 | 6.89 | 6.12 | 14.76 |  13.03* | 24.86
LCE | 0.23 | 0.13 | 6.71 | 8.90 | 38.54 | 26.25 | 51.48 | 44.29
LCL | 0.29 | 0.31 | 34.28 | 22.93 | 104.10 | 72.95 | 130.80 | 94.98
| LASSO | 0.25 | 0.20 | 7.08 | 11.91 | 17.15 | 25.04 | 29.44 | 36.69
| SCAD | 0.21 | 0.16 | 3.07 | 9.02 | 6.62 | 18.80 | 14.88 | 33.00
| MCP | 0.22 | 0.18 | 2.82 | 8.92 | 6.35 | 19.38 | 15.13 | 33.52

Table 3: 1000 times the AMSE for estimated coefficients over 100 replications;
BOLD represent group minimum; * represent overall minimum.

We observe that as the true model size $p_{\boldsymbol{\gamma}_{T}}$ and the
true correlation between covariates increases, the AMSE increases for all
techniques. However, similar to model selection performance, LPEP versions
significantly outperforms all other techniques in terms of estimation
performance under non-null true model scenarios and is comparable to other
techniques when the true model is the null model.

##  6 Real data applications

###  6.1 URINARY: Determinants of urinary incontinence

The URINARY data set describes the results from a small drug study with 21
subjects. The response corresponds to whether the subject developed urinary
incontinence after receiving the drug. The explanatory variables capture drug-
induced physiological changes, which were in the same direction for most
subjects. This data set was first presented in Potter (2005), and is further
discussed in Mansournia et al. (2018). While very small, the data set is
challenging to analyze because it exhibits full separation. In particular, the
maximum likelihood estimates of the regression coefficients are all infinite
(please see the top row of Table 4), indicating that separation is not induced
by any of the variables on its own, but by a non-trivial linear combination of
them.

In addition to the maximum likelihood estimates reported by the R function
glm, Table 4 presents estimates for the regression coefficients for various
Bayesian and penalized likelihood methods. The results for LPEP, CRPEP and
DRPEP are based on 10,000 iterations of the MCMC algorithm obtained after a
burn-in period of 10,000 iterations. On the other hand, for LCL we use full
model enumeration procedure in the R package BAS. In the case of Bayesian
procedures, Table 4 presents model-averaged posterior means, as well as 95%
model-averaged credible intervals for the coefficients. Note that confidence
intervals for the penalized likelihood procedures are not presented since they
are not straightforward to obtain and the R packages we used to fit these
models do not readily provide them. Furthermore, results for CRPEP and DRPEP
are not included under the robust hyper-prior because such a procedure is not
implemented in Fouskakis et al. (2018).

Note that LCL produces large point estimates and very wide credible for the
model coefficients under all hyperpriors. This is no surprise; the prior
proposed by Li & Clyde (2018) is proper only for models for which the maximum
likelihood estimates are finite. This means that, for a data set like URINARY,
some of the Bayes factors associated with LCL are ill-defined. This is also
the reason why we do not show results for LCE; the target posterior
distribution for the associated Markov chain Monte Carlo algorithm is improper
if the full model is included in the analysis. Furthermore, note that the
point estimates associated with CRPEP appear to be different from those
generated by LPEP, DRPEP, and the penalized likelihood methods. This is
clearer when looking at the intercept of the model, which is negative with
high probability under CRPEP but positive with high probability under LPEP and
DRPEP under all hyperpriors.

| $\beta_{0}$ | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$
---|---|---|---|---
| MLE | -83.84 | -2445.04 | -1653.76 | 310.27
| (-1969.33 , 484.67) | (-53259.42 , 10866.11) | (-34653.95 , 5062.34) | (-1166.03 , 6803.89)
$\delta=n$ | LPEP | 0.56 | -0.70 | -0.39 | 0.15
(-1.66 , 2.85) | (-2.32 , 0.10) | (-0.81 , -0.13) | (0.00 , 0.37)
LCL | -83.84 | -2333.88 | -1578.58 | 296.17
(-6009.26 , 5897.13) | (-161488.87 , 158312.76) | (-109266.01 , 107118.14) | (-19896.81 , 20678.41)
CRPEP | -1.15 | -0.70 | -0.34 | 0.00
(-3.21 , 0.52) | (-1.88 , 0.30) | (-0.63 , -0.10) | (-0.00 , 0.00)
DRPEP | 0.69 | -1.00 | 0.00 | 0.06
(-0.55 , 2.13) | (-2.09 , -0.19) | (0.00 , 0.00) | (-0.03 , 0.16)
$\delta\sim\text{robust}$ | LPEP | 0.71 | -0.98 | -0.52 | 0.19
(-1.74 , 3.55) | (-3.82 , 0.05) | (-1.89 , -0.12) | (0.00 , 0.54)
LCL | -83.84 | -2148.67 | -1453.30 | 272.66
(-6250.14 , 5980.74) | (-161065.38 , 154146.29) | (-108979.51 , 104298.98) | (-19890.08 , 20102.78)
$\delta\sim\text{hyper g/n}$ | LPEP | 0.61 | -0.75 | -0.41 | 0.15
(-1.65 , 3.07) | (-2.74 , 0.10) | (-1.02 , -0.09) | (0.00 , 0.39)
LCL | -83.84 | -1288.82 | -871.72 | 163.55
(-6252.44 , 5650.45) | (-124412.67 , 113166.12) | (-84179.77 , 76570.77) | (-15457.94 , 14685.16)
CRPEP | -1.04 | -0.66 | -0.33 | 0.00
(-3.04 , 0.62) | (-1.78 , 0.30) | (-0.65 , -0.08) | (-0.00 , 0.00)
DRPEP | -0.89 | 0.00 | -0.36 | 0.00
(-3.11 , 0.69) | (0.00 , -0.00) | (-0.76 , -0.11) | (-0.00 , 0.00)
| LASSO | 0.36 | -0.70 | -0.31 | 0.11
| SCAD | 0.41 | -0.23 | -0.20 | 0.07
| MCP | 0.40 | -0.17 | -0.20 | 0.07

Table 4: Estimated BMA coefficients and 95% credible intervals for different
Bayesian techniques for urinary dataset; For frequentist techniques, estimated
coefficient is displayed.

Next, we present in Table 5 the posterior probabilities associated with each
of the eight models under consideration under each one of the Bayesian
approaches, along with the model selected by each of the penalized likelihood
methods. In all cases, LCL consistently places probability one on the full
model, which is also the model selected by all the penalized likelihood
methods. The full model is also consistently preferred by LPEP, but there is
more uncertainty. Indeed, under LPEP the model that excludes the first
covariate receives between 0.19 and 0.22 probability, and the model that
excludes the third covariate is assigned between 0.03 and 0.04 posterior
probability. In contrast, CRPEP and DRPEP place zero probability on the full
model. Instead, CRPEP consistently favors the model that excludes the third
covariate, while DRPEP yields contradictory results depending on the
hyperprior: it favors the model that excludes the second variable under the
$g=n$ hyperprior, but the model that only includes the second variable under
the hyper-g/n prior.

| $(\gamma_{1},\gamma_{2},\gamma_{3})$
---|---
| $(0,0,0)$ | $(1,0,0)$ | $(0,1,0)$ | $(1,1,0)$ | $(0,0,1)$ | $(1,0,1)$ | $(0,1,1)$ | $(1,1,1)$
$\delta=n$ | LPEP | 0.00 | 0.00 | 0.06 | 0.03 | 0.00 | 0.00 | 0.22 | 0.69
LCL | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
CRPEP | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00
DRPEP | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00
$\delta\sim\text{robust}$ | LPEP | 0.00 | 0.00 | 0.06 | 0.03 | 0.00 | 0.00 | 0.19 | 0.71
LCL | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
$\delta\sim\text{hyper g/n}$ | LPEP | 0.00 | 0.00 | 0.06 | 0.04 | 0.00 | 0.00 | 0.20 | 0.69
LCL | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
CRPEP | 0.00 | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00
DRPEP | 0.00 | 0.00 | 1.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
| LASSO | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
| SCAD | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00
| MCP | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00

Table 5: Posterior model probabilities for all the models in the model space
for urinary dataset.

###  6.2 GUSTO-I: Survival to treatments for occluded coronary arteries

Next, we consider data from the Global Utilization of Streptokinase and TPA
for Occluded Coronary Arteries (GUSTO-I) trial (Califf et al., 1996), which
has been previously analyzed in Held et al. (2015) and Li & Clyde (2018), and
is publicly available at http://www.clinicalpredictionmodels.org/ (Steyerberg
et al., 2019). Similar to previous analyses, we model the binary endpoint of
30-day survival for a subgroup of $n=2188$ patients using 17 clinical
covariates described in the supplementary materials.

Figure 3 displays the marginal posterior inclusion probabilities (PIPs) for
Bayesian methods and the inferred model under the penalized likelihood
techniques. For techniques related to LCL, we again rely on full model
enumeration. On the other hand, for all other Bayesian techniques, we use
131,000 iterations with a burn-in of 10,000 iterations.

Similar to simulation studies, all penalized likelihood techniques select
denser models than the Bayesian procedures. In line with Li & Clyde (2018) and
Held et al. (2015), we observe that AGE, KILLIP, HYP, HRT and STE have high
PIPs under all methods. However, it is worthwhile noting that the different
versions of LCL perform quite differently. In particular, the version of LCL
that relies on a hyper-g/n hyperprior tends to explore very dense models
leading, to PIPs close to $0.5$ for all variables. Similarly, the hyper-g/n
versions of CRPEP and DRPEP seem to differ from their $g=n$ versions with
respect to PMI and SEX variables. On the other hand, the different versions of
the LPEP prior are roughly in agreement for all variables

Figure 3: Marginal posterior inclusion probabilities (PIPs) for GUSTO-I
dataset (Bayesian procedures) and variables included in the model (penalized
likelihood methods).

As in Li & Clyde (2018), we also compare the different procedures in terms of
their out-of-sample predictive performance. For this purpose, we performed a
10-fold cross-validation study. More specifically, we divided the data into 10
folds, training our model on 9 them and evaluating the predictive accuracy on
the remaining one. Table 6 presents the average value of four different
metrics across all 10 folds. The metrics we rely on are the same ones employed
in Li & Clyde (2018): the area under the ROC curve (AUC), the Calibration
Slope (CS), the Logarithmic Score (LS) and the Brier score (BRIER). AUC and CS
allow us to evaluate the methods in terms of discrimination and calibration.
In both cases, scores closer to 1 indicate better performance. On the other
hand, LS and BRIER measure the predictive accuracy of methods; in both cases
lower scores indicate better performance.

Table 6 suggests that most methods perform similarly. The main exceptions are
both versions of CRPEP and DRPEP, which seem to substantially underperform
across all metrics. LPEP procedures slightly outperforming other methods in
terms of AUC and CS. On the other hand, LASSO seems to slightly outperform the
LPEP procedures in terms of LS and Brier score, but at the cost of selecting
much denser models.

| AUC | CS | LS | BRIER
---|---|---|---|---
$\delta=n$ | LPEP |  0.8324* | 0.9971 | 0.1824 | 0.0496
LCL | 0.8300 | 0.9931 | 0.1831 | 0.0497
CRPEP | 0.7789 | 1.0578 | 0.1965 | 0.0521
DRPEP | 0.7790 | 1.0569 | 0.1963 | 0.0521
$\delta\sim\text{robust}$ | LPEP | 0.8322 | 1.0129 | 0.1822 | 0.0495
LCL | 0.8316 | 0.9804 | 0.1822 | 0.0495
$\delta\sim\text{hyper g/n}$ | LPEP | 0.8319 |  1.0074* | 0.1823 | 0.0495
LCL | 0.8311 | 1.0109 | 0.1818 | 0.0493
CRPEP | 0.7956 | 1.1677 | 0.1951 | 0.0522
DRPEP | 0.7800 | 1.0571 | 0.1961 | 0.0520
| LASSO | 0.8305 | 1.0369 |  0.1816* |  0.0492*
| SCAD | 0.8243 | 0.9135 | 0.1838 | 0.0496
| MCP | 0.8250 | 0.9196 | 0.1838 | 0.0496

Table 6: Average prediction accuracy measures in a 10-fold cross validation
study for GUSTO-I dataset; Bold represents group maximum for AUC, for CS
closest to one, and group minimum for LS and Brier score; * represents the
corresponding best score among all methods.

###  6.3 HOUSE107: Determinants of legislator behavior in 107th U.S. House of
Representatives

Most analyses of congressional voting treat all roll-call votes in the same
way, independently of the type of vote. This might mask considerable variation
in voting behavior across different types of votes. For example, Jessee &
Theriault (2012) provide empirical evidence that the forces affecting
legislators’ voting on procedural and final passage matters have exhibited
important changes over time, with differences between these two vote types
becoming larger, particularly in recent congresses.

One shortcoming of the methodology presented in Jessee & Theriault (2012) is
that it provides legislature-wide measures of agreement among vote types, but
cannot ascribe observed differences to individual legislators. Recently,
Lofland et al. (2017) and Moser et al. (2021) developed methodology that
enables the identification of differences in voting behavior across votes
types for individual legislators. In this section we analyze a dataset where
the response variable corresponds to estimates of whether each legislator in
the 107th U.S. House of Representatives share the same voting behavior across
final passage, amendment and procedural votes. These estimates are obtained
using the model introduced in Moser et al. (2021). Hence, in this case,
$y_{i}=1$ if the $i$-th legislator voting preferences remain unchanged across
all three vote types, and $y_{i}=0$ otherwise. The goal of this analysis is to
understand whether a group of 26 characteristics of the legislator or its
constituency affect the likelihood of such changes. Linking voting behavior
with these underlying characteristics can provide important insights into the
workings of a political system (e.g., see Facchini & Steinhardt, 2011, O’Roark
& Wood, 2011 and Cragg et al., 2013).

Figure 4: Marginal posterior inclusion probabilities (PIPs) for HOUSE107
dataset (Bayesian procedures) and variables included in the model (penalized
likelihood methods).

Detailed description of the variables is available in supplementary materials.
All the results for the Bayesian procedures in this section are based on the
same settings and number of iterations as those used in the previous Sub-
section for the GUSTO-I dataset. Figure 4 shows the marginal PIPs for all 26
variables, along with a list of variables selected by the penalized likelihood
methods. Note that, as in previous illustrations, the penalized likelihood
methods tend to select a superset of the variables selected by the Bayesian
approaches, with LASSO selecting the largest superset. Furthermore, there is
broad agreement among all Bayesian procedures. For example, all techniques,
except both versions of DRPEP, assign high PIP to daysServed. Based on our
experience with this application, including daysServed in the model is
sensible. For example, length of tenure has been previously identified in the
literature as an important predictor of legislator’s effectiveness (e.g., see
Miquel & Snyder Jr, 2006). One place where the different procedures do seem to
disagree is whether medianIncome (the median income in the district
represented by the legislator) explains voting behavior. All LPEP procedures
agree in providing weak to moderate evidence against the inclusion of this
variable, while most other procedures provide weak to moderate evidence in
favor of it. Interestingly, the two versions of the DRPEP seem to disagree,
with the unit information version providing weak evidence in favor of its
inclusion and the hyper-g/n providing weak evidence against it.

##  7 Discussion

The results from our theoretical and empirical studies show that the LPEP
approach to creating non-informative priors for logistic regression is at
least competitive, and in some cases superior, to existing techniques in terms
of model selection and parameter estimation performance. The differences are
particularly striking when comparing the LPEP with the original CRPEP and
DRPEP approaches proposed in Fouskakis et al. (2018). Aside from the increased
computational expense that makes the practical application of CRPEP and DRPEP
prior challenging, our empirical results show that the behavior of these
procedures can be greatly affected by the choice of hyperpriors. Furthermore,
both procedures can lead to unexpected results that substantially deviate from
the consensus of other Bayesian and non-Bayesian methods.

We were surprised by the poor behavior of some of the LCL procedures in some
of our simulation studies. One point to note is that the setup of the
simulations in Section 5 ($n=500$, $p=100$) was only briefly studied in Li &
Clyde (2018). Indeed, most of the simulation studies in Li & Clyde (2018)
focus on settings involving fewer covariates ($p=20$), which are similar to
those in the simulation study that we present in our supplementary materials.
In this lower dimensional setting, LCL behaves quite well. Disentangling the
role of the prior distribution from that of the various approximations used to
speed up computation in the BAS package is one particular challenge in
understanding these negative results. However, by comparing the results for
LCL against its “exact” version (LCE), as well as those under the unit
information ($\delta=n$) with those under the robust and hyper-g/n versions of
the procedure, it would seem like the use of the Laplace approximation for the
Bayes factor is not the culprit. Instead, these results seem to driven by a
combination of some sensitivity to the choice of hyperprior for $\delta$ and
issues with the default approximation procedure implemented in BAS to
integrate over $\delta$. Interestingly, the sensitivity to the hyperprior does
not seem to be present for the LPEP procedures. We believe that this is an
area that is in need of further investigation in the future.

This paper focuses mostly on developing the LPEP for logistic regression.
However, the formulation is very general and can be extended to many other
generalized linear models. Many of the computational advantages of our
procedure extend to binomial, negative binomial and multinomial logic models
where the data augmentation approach of Polson et al. (2013) can be readily
applied. This is also true for probit models in which computation can rely on
the data augmentation approach of Albert & Chib (1993), as well as for
loglinear regression using the approach of Frühwirth-Schnatter et al. (2009).
These extensions will be explored elsewhere.

One final brief note about our theoretical results. Our asymptotic analyses
assume that the number of covariates $p$ is fixed with the sample size $n$. We
believe that our results can be extended to situations in which $p$ grows with
$n$ (as long $n$ remains larger than $p$), potentially along the lines of
Barber et al. (2016). However, this requires the careful introduction of
additional constraints that ensure that the maximum likelihood estimates under
the training sample remains finite as both $p$ and $n$ (and therefore, as
$n^{*}$ grow). We plan to explore this question in our future work.

## Additional materials

The data sets used in the simulation studies in Section 5 along with real data
sets (for e.g. Urinary, endometrial, GUSTO, and House107) and code for Laplace
PEP methodology discussed in the paper with options to implement LCE and Bové
et al. (2011)’s prior are available at https://github.com/Anupreet-
Porwal/LPEP. Code to replicate the results in Section 5 and 6 is available at
https://github.com/Anupreet-Porwal/LPEP-Paper-Analysis. Supplementary
materials, which include a further set of simulation studies along with the
analysis of the endometrial and PIMA data set (Heinze & Schemper, 2002;
Agresti, 2015) is available from the authors.

## Acknowledgements

We would like to thank Dimitris Fouskakis, Ioannis Ntzourfras and Konstantinos
Perrakis for sharing their code for the CRPEP and DRPEP, which we used in our
analyses.

##  Appendix A Proof of Theorem 3.1 (Existence of MLEs)

The two conditions in the theorem together imply that the maximum likelihood
estimator for the full model,

$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{F}}=\operatorname*{arg\,max}_{\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}}}\,\ell_{\boldsymbol{\gamma}_{F}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}}),$

exists, is finite, and, furthermore, is unique (e.g., see Mäkeläinen et al.,
1981).

Now, note that for any other model $\boldsymbol{\gamma}$, maximizing
$\ell_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ is
equivalent to maximizing
$\ell_{\boldsymbol{\gamma}_{F}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}})$
subject to the constraint $\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}}\in
S_{\boldsymbol{\gamma}}$, where
$S_{\boldsymbol{\gamma}}=\\{\boldsymbol{\theta}\in
S_{\boldsymbol{\gamma}_{F}}:\theta_{j}=0\mbox{ if }\gamma_{j}=0\\}$ and that
$S_{\boldsymbol{\gamma}}$ is also an open connected set for all
$\boldsymbol{\gamma}$. Because
$\ell_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}_{F}})$ is
continuous and strongly concave, then its restriction to
$S_{\boldsymbol{\gamma}}$ is also continuous and strongly concave for any
$\boldsymbol{\gamma}$. Furthermore, we also have
$\lim_{\boldsymbol{\beta}_{\boldsymbol{\gamma}}\to\partial
S_{\boldsymbol{\gamma}}}\ell_{\boldsymbol{\gamma}}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})=-\infty$
because $\partial S_{\boldsymbol{\gamma}}\subset\partial
S_{\boldsymbol{\gamma}_{F}}$. Therefore,
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}$ also exists and is finite and
unique for any $\boldsymbol{\gamma}$. $\Box$

##  Appendix B Proof of Theorem 4.1 (Tail behavior)

We develop the argument only for
$\pi_{\boldsymbol{\gamma}}^{LPEP-R}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$.
The proof for $\pi_{\boldsymbol{\gamma}}^{LPEP-
HGN}(\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ follows along almost identical
lines. Note that

$\displaystyle\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=\sum_{{\boldsymbol{y}}^{*}}m^{*}(\mathbf{y}^{*})\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}},$

where

$\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*})=\left.\pi_{\boldsymbol{\gamma}}^{LPEP-R}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*})\right|_{\boldsymbol{\beta}_{\boldsymbol{\gamma}}=s{\boldsymbol{v}}}\\\
=\int_{0}^{\infty}\left(\frac{1}{2\pi\delta}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*})\right|^{1/2}\exp\left\\{-\frac{1}{2\delta}\left(s{\boldsymbol{v}}-\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right)^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*})\left(s{\boldsymbol{v}}-\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\right)\right\\}f^{R}(\delta|\boldsymbol{\gamma})d\delta$

is conditional on a given training sample ${\boldsymbol{y}}^{*}$. (We can
exchange the summation and the limit in this case because, for any $n^{*}$,
the number of potential training samples is finite.) In the sequel, it will
also be important to remember that $m^{*}({\boldsymbol{y}}^{*})$ is defined so
that the maximum likelihood estimators exist for any sample
${\boldsymbol{y}}^{*}$.

To simplify notation, define

$\zeta^{R-C}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*})=\int_{0}^{\infty}\left(\frac{1}{2\pi\delta}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*})\right|^{1/2}\exp\left\\{-\frac{1}{2\delta}s^{2}{\boldsymbol{v}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*}){\boldsymbol{v}}\right\\}f^{R}(\delta|\boldsymbol{\gamma})d\delta.$

Note that

$\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=\\\
\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\zeta^{R-C}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*})}\times\lim_{s\to\infty}\frac{\left(s^{2}\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}\times\lim_{s\to\infty}\frac{\zeta^{R-C}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*})}{\left(s^{2}\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}.$

Clearly, the first two limits converge to finite functions that
${\boldsymbol{v}}$, the training sample ${\boldsymbol{y}}^{*}$ and/or the
model $\boldsymbol{\gamma}$. Hence,

$\displaystyle\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}$
$\displaystyle=c_{\boldsymbol{\gamma},1}({\boldsymbol{v}},y^{*})\lim_{s\to\infty}\left(s^{2}\right)^{\frac{p_{\boldsymbol{\gamma}}+2}{2}}\zeta^{R-C}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*}),$

where $0<c_{\boldsymbol{\gamma},1}({\boldsymbol{v}},y^{*})<\infty$. Plugging
in $\zeta^{R-C}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},y^{*})$ and
$f^{R}(\delta|\boldsymbol{\gamma})$, we can write

$\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=c_{\boldsymbol{\gamma},1}({\boldsymbol{v}},y^{*})\lim_{s\to\infty}\left(s^{2}\right)^{\frac{p_{\boldsymbol{\gamma}}+2}{2}}\times\\\
\int_{0}^{\infty}\left(\frac{1}{2\pi\delta}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*})\right|^{1/2}\exp\left\\{-\frac{1}{2\delta}s^{2}{\boldsymbol{v}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*}){\boldsymbol{v}}\right\\}\times\\\
\frac{1}{2(p_{\boldsymbol{\gamma}}+1)^{1/2}}\frac{(n^{*}+1)^{1/2}}{(\delta+1)^{3/2}}\mathbf{1}\left(\delta>\frac{n^{*}-p_{\boldsymbol{\gamma}}}{p_{\boldsymbol{\gamma}}+1}\right)d\delta.$

Substituting
$\lambda=\left(\frac{n^{*}+1}{p_{\boldsymbol{\gamma}}+1}\right)\frac{1}{\delta+1}$,
we can write above equation as

$\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=c_{\boldsymbol{\gamma},1}({\boldsymbol{v}},y^{*})\left(\frac{1}{2\pi}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*})\right|^{1/2}\times\\\
\lim_{s\to\infty}\left(s^{2}\right)^{\frac{p_{\boldsymbol{\gamma}}+2}{2}}\int_{0}^{1}\left(\frac{\lambda}{m-\lambda}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\lambda^{-\frac{1}{2}}\exp\left\\{-\left(\frac{\lambda}{m-\lambda}\right)qs^{2}\right\\}d\lambda$

where
$q=\frac{1}{2}{\boldsymbol{v}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}(\mathbf{y}^{*}){\boldsymbol{v}}$,
$m=\frac{n^{*}+1}{p_{\boldsymbol{\gamma}}+1}>1$ since
$n^{*}>p_{\boldsymbol{\gamma}}$. Now, from Lemma 2 in Bayarri et al. (2012),

$\lim_{s\to\infty}\left(s^{2}\right)^{\frac{p_{\boldsymbol{\gamma}}+2}{2}}\int_{0}^{1}\left(\frac{\lambda}{m-\lambda}\right)^{\frac{p_{\boldsymbol{\gamma}}+1}{2}}\lambda^{-\frac{1}{2}}\exp\left\\{-\left(\frac{\lambda}{m-\lambda}\right)qs^{2}\right\\}d\lambda=c_{\boldsymbol{\gamma},2}({\boldsymbol{v}},y^{*}),$

where $0<c_{\boldsymbol{\gamma},2}({\boldsymbol{v}},y^{*})<\infty$, and
therefore

$\lim_{s\to\infty}\frac{\zeta^{R}(s\mid{\boldsymbol{v}},\boldsymbol{\gamma},{\boldsymbol{y}}^{*})}{\left(1+s^{2}/(p_{\boldsymbol{\gamma}}+1)\right)^{-\frac{p_{\boldsymbol{\gamma}}+2}{2}}}=c_{\boldsymbol{\gamma},1}({\boldsymbol{v}},y^{*})c_{\boldsymbol{\gamma},2}({\boldsymbol{v}},y^{*})=c_{\boldsymbol{\gamma},3}({\boldsymbol{v}},y^{*}).$

To complete the proof simply define
$c_{\boldsymbol{\gamma}}({\boldsymbol{v}})=\sum_{{\boldsymbol{y}}^{*}}m^{*}(\mathbf{y}^{*})c_{\boldsymbol{\gamma},3}({\boldsymbol{v}},y^{*})$.
Since we have a finite number of terms in the sum and each term is both
positive and finite, so is $c_{\boldsymbol{\gamma}}({\boldsymbol{v}})$.

##  Appendix C Proof of Theorem 4.2 (Intrinsic consistency)

Note that if ${\boldsymbol{y}}^{*}$ leads to separability, then
${\boldsymbol{y}}^{**}=\mathbf{1}-{\boldsymbol{y}}^{*}$ does as well. Hence,
our choice of $m^{*}\left({\boldsymbol{y}}^{*}\mid{\boldsymbol{X}}\right)$ is
symmetric and, as $n^{*}\to\infty$,
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*})\overset{p}{\to}\mathbf{0}$
under $m^{*}$. In turn, this implies that
$\hat{\theta}_{\boldsymbol{\gamma},i}({\boldsymbol{y}}^{*})\overset{p}{\to}1/2$
for all $i$, and therefore
$\frac{1}{n}{\boldsymbol{X}}_{\boldsymbol{\gamma}}^{\prime}{\boldsymbol{W}}_{\boldsymbol{\gamma}}({\boldsymbol{y}}^{*}){\boldsymbol{X}}_{\boldsymbol{\gamma}}\overset{p}{\to}\frac{1}{4}\boldsymbol{\Sigma}_{\boldsymbol{\gamma}}$.
This completes the proof when $g=n$.

In the case where $\delta$ is given a prior distribution define
$\delta=n^{*}\delta^{*}$. Then, under $f^{HGn}(\delta)$, $\delta^{*}$ has
density

$f^{HGn}(\delta^{*})=\left(1+\delta^{*}\right)^{-2}$

which is a proper, non-degenerate prior. The argument for the robust prior
follows along similar lines.

$\Box$

##  Appendix D Proof of Theorem 4.3 (Model selection consistency)

Since

$Pr(\boldsymbol{\gamma}=\boldsymbol{\gamma}_{T}\mid{\boldsymbol{y}})=\frac{1}{1+\sum_{\boldsymbol{\gamma}\neq\boldsymbol{\gamma}_{T}}\frac{f(\boldsymbol{\gamma})}{f(\boldsymbol{\gamma}_{T})}\frac{m^{LPEP}_{\boldsymbol{\gamma}}({\boldsymbol{y}})}{m^{LPEP}_{\boldsymbol{\gamma}_{T}}({\boldsymbol{y}})}}$

it is enough to show that

$\frac{m^{LPEP}_{\boldsymbol{\gamma}}({\boldsymbol{y}})}{m^{LPEP}_{\boldsymbol{\gamma}_{T}}({\boldsymbol{y}})}\xrightarrow[n\to\infty]{P}0.$

The proof follows along similar lines as that in Li & Clyde (2018). We start
by assuming similar regularity conditions:

  * (i)

The true model $\boldsymbol{\gamma}_{T}$ is among the $2^{p}$ models under
consideration, with $p$ fixed.

  * (ii)

For every $i=1,2,\ldots$, the vector ${\boldsymbol{x}}_{i}$ is such that
$\left\|{\boldsymbol{x}}_{i}\right\|_{2}$ is bounded by a constant.

  * (iii)

For all $n$, the smallest eigenvalue of
$\frac{1}{n}{\boldsymbol{X}}^{T}{\boldsymbol{X}}$ is bounded from below by a
positive constant.

Note that these conditions imply the weak consistency and asymptotic normality
of the maximum likelihood estimators for all models under consideration (e.g.,
see Fahrmeir & Kaufmann, 1985).

Consider first
$m_{\boldsymbol{\gamma}}({\boldsymbol{y}}\mid{\boldsymbol{y}}^{*})=\int
f_{\boldsymbol{\gamma}}({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})\phi_{p_{\boldsymbol{\gamma}}+1}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\right),n{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{-1}\left({\boldsymbol{y}}^{*}\right)\right)d\boldsymbol{\beta}_{\boldsymbol{\gamma}}$.
Using a Laplace approximation of
$f_{\boldsymbol{\gamma}}({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})$,

$\displaystyle
f_{\boldsymbol{\gamma}}({\boldsymbol{y}}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}})$
$\displaystyle=f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)\exp\left\\{-\frac{1}{2}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}-\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}-\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)\right\\}\exp\left\\{R\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}},\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)\right\\}$

where
$R\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}},\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}\right)$
is the residual form the second order Taylor expansion of the loglikelihood.
Note that, to simplify notation, we have let
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}=\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\right)$
and
${\boldsymbol{H}}_{\boldsymbol{\gamma}}={\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\right)$.

Hence, following (Tierney & Kadane, 1986),

$m_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid{\boldsymbol{y}}^{*}\right)=f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)n^{-\frac{p_{\boldsymbol{\gamma}}+1}{2}}\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}})\right|}\right)^{1/2}\\\
\exp\left\\{-\frac{1}{2}\left[\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}+\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}_{\boldsymbol{\gamma}}^{-1}{\boldsymbol{d}}_{\boldsymbol{\gamma}}\right]\right\\}\left(1+\mathcal{O}\left(\frac{1}{n}\right)\right)$

where

$\displaystyle{\boldsymbol{d}}_{\boldsymbol{\gamma}}$
$\displaystyle={\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{*},$
$\displaystyle{\boldsymbol{E}}_{\boldsymbol{\gamma}}$
$\displaystyle={\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}.{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}$

As before, we simplify notation by letting
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{*}=\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\right)$,
and
${\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}={\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}^{*}\right)$.
From this,

$\frac{m^{LPEP}_{\boldsymbol{\gamma}}({\boldsymbol{y}})}{m^{LPEP}_{\boldsymbol{\gamma}_{T}}({\boldsymbol{y}})}=\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)}{f_{\boldsymbol{\gamma}_{T}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}\right)}n^{(p_{\boldsymbol{\gamma}_{T}}-p_{\boldsymbol{\gamma}})/2}\\\
\frac{\mathsf{E}_{y^{*}}\left\\{\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}})\right|}\right)^{1/2}\exp\left\\{-\frac{1}{2}\left[\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}+\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}\right]\right\\}\right\\}}{\mathsf{E}_{y^{*}}\left\\{\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}_{T}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}_{T}})\right|}\right)^{1/2}\exp\left\\{-\frac{1}{2}\left[\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}_{T}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}+\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}_{T}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}_{T}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{d}}_{\boldsymbol{\gamma}_{T}}\right]\right\\}\right\\}}\,\mathcal{O}(1)$

Consider first the term

$\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)}{f_{\boldsymbol{\gamma}_{T}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}\right)}n^{(p_{\boldsymbol{\gamma}_{T}}-p_{\boldsymbol{\gamma}})/2}.$

This is just Schwartz criterion, which is well known to be consistent in this
setting. Indeed, under the regularity conditions (i)-(iii), Lemma A.3 in the
supplementary materials of Li & Clyde (2018) can be applied to show that, as
$n$ increases, the likelihood ratio
$\Lambda_{\boldsymbol{\gamma},\boldsymbol{\gamma}_{T}}=\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)}{f_{\boldsymbol{\gamma}_{T}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}\right)}$
has the following behavior.

  * (a)

If $\boldsymbol{\gamma}_{T}\subset\boldsymbol{\gamma}$, then
$\Lambda_{\boldsymbol{\gamma},\boldsymbol{\gamma}_{T}}=\mathcal{O}_{P}(1)$.

  * (b)

If $\boldsymbol{\gamma}_{T}\not\subset\boldsymbol{\gamma}$ then
$\Lambda_{\boldsymbol{\gamma},\boldsymbol{\gamma}_{T}}=\mathcal{O}_{P}\left(e^{-c_{\boldsymbol{\gamma}}n}\right)$
for some positive constant $c_{\boldsymbol{\gamma}}$.

Hence, if $\boldsymbol{\gamma}_{T}\subset\boldsymbol{\gamma}$ then necessarily
$p_{\boldsymbol{\gamma}_{T}}<p_{\boldsymbol{\gamma}}$ and

$\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)}{f_{\boldsymbol{\gamma}_{T}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}\right)}n^{(p_{\boldsymbol{\gamma}_{T}}-p_{\boldsymbol{\gamma}})/2}=\mathcal{O}_{P}(n^{-c})$

for some positive constant $c$ and therefore it tends to zero as $n\to\infty$.
On other hand, $\boldsymbol{\gamma}_{T}\not\subset\boldsymbol{\gamma}$

$\frac{f_{\boldsymbol{\gamma}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right)}{f_{\boldsymbol{\gamma}_{T}}\left({\boldsymbol{y}}\mid\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}\right)}n^{(p_{\boldsymbol{\gamma}_{T}}-p_{\boldsymbol{\gamma}})/2}=\mathcal{O}_{P}\left(e^{-c_{\boldsymbol{\gamma}}n}\right)$

for some positive constant $c_{\boldsymbol{\gamma}}$ no matter whether
$p_{\boldsymbol{\gamma}_{T}}<p_{\boldsymbol{\gamma}}$ or
$p_{\boldsymbol{\gamma}_{T}}\geq p_{\boldsymbol{\gamma}}$, and again it
converges to zero as $n\to\infty$.

All that remains now is to show that

$\frac{\mathsf{E}_{y^{*}}\left\\{\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}})\right|}\right)^{1/2}\exp\left\\{-\frac{1}{2}\left[\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}+\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}\right]\right\\}\right\\}}{\mathsf{E}_{y^{*}}\left\\{\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}_{T}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}_{T}})\right|}\right)^{1/2}\exp\left\\{-\frac{1}{2}\left[\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}_{T}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}_{T}}+\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}_{T}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}_{T}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}_{T}}{\boldsymbol{d}}_{\boldsymbol{\gamma}_{T}}\right]\right\\}\right\\}}$
(12)

is, at most, $\mathcal{O}_{P}(1)$. To do this, we consider the behavior of

$\left(\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}})\right|}\right)^{1/2}\exp\left\\{-\frac{1}{2}\left[\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}+\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}^{T}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}\right]\right\\}$

in probability as both ${\boldsymbol{y}}$ and ${\boldsymbol{y}}^{*}$ grow.
Because of dominated convergence, the behavior in probability with respect to
${\boldsymbol{y}}^{*}$ is the same as that in expectation.

  * •

First, note that
$\frac{\left|{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right|}{\left|{\boldsymbol{E}}_{\boldsymbol{\gamma}}\right|}=\left|\mathbf{I}+\frac{1}{n}{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\right|^{-1}=\left|\mathbf{I}+\frac{1}{n}\left(n{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}\right)\left(\frac{1}{n}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\right)\right|^{-1}$.
Under regularity conditions (ii) and (iii),
$\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\xrightarrow[n\to\infty]{P}{\boldsymbol{J}}_{\boldsymbol{\gamma}}$
and
$\frac{1}{n}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\xrightarrow[n\to\infty]{P}{\boldsymbol{J}}^{*}_{\boldsymbol{\gamma}}$
where ${\boldsymbol{J}}_{\boldsymbol{\gamma}}$ and
${\boldsymbol{J}}^{*}_{\boldsymbol{\gamma}}$ are
$(p_{\boldsymbol{\gamma}}+1)\times(p_{\boldsymbol{\gamma}}+1)$ constant and
strictly positive definite matrices with finite entries (please see below).
Hence,

$\left|\mathbf{I}+\frac{1}{n}{\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\right|^{-1}\xrightarrow[n\to\infty]{P}1.$

We only elaborate on the proof that
$\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\xrightarrow[n\to\infty]{P}{\boldsymbol{J}}_{\boldsymbol{\gamma}}$;
the argument for $\frac{1}{n}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}$ is
analogous. Furthermore, we focus on the diagonal elements of the matrix since
the off-diagonal elements are bounded by the diagonal ones because the matrix
is, by construction at least semi-positive definite. Since
$\left\|{\boldsymbol{x}}_{i}\right\|_{2}<c$ for all $i$ and some constant $c$,
we have $|x_{i,j}|<c$. Therefore
$\left[\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right]_{j,j}=\frac{1}{n}\sum_{i=1}^{n}\hat{w}_{\boldsymbol{\gamma},i}x_{i,j}^{2}<\frac{c^{2}}{n}\sum_{i=1}^{n}\frac{\left(1+\exp\left\\{-{\boldsymbol{x}}_{i}^{T}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right\\}\right)^{2}}{\exp\left\\{-{\boldsymbol{x}}_{i}^{T}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right\\}}$.
Now, the maximum likelihood estimator
$\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}$ is not necessarily consistent
when $\boldsymbol{\gamma}\neq\boldsymbol{\gamma}_{T}$, but it does converge in
probability to a finite constant (e.g., see Fahrmeir, 1990). Combined with the
fact that the ${\boldsymbol{x}}_{i}$s have compact support, this implies that
$\left(1+\exp\left\\{-{\boldsymbol{x}}_{i}^{T}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\right\\}\right)^{-1}$
converges in probability to a constant $\theta^{(0)}_{\boldsymbol{\gamma},i}$
that, for all $i$, is bounded away from 0 and 1, i.e.,
$0<a\leq\theta^{(0)}_{\boldsymbol{\gamma},i}\leq b<1$ for some constants $a$
and $b$. Therefore,
$\left[\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right]_{j,j}<\frac{c^{2}}{\min\\{a(1-a),b(1-b)\\}}<\infty$.

  * •

Next, consider
$\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}$.
Because of our choice for $m^{*}({\boldsymbol{y}}^{*})$, the results in
Fahrmeir (1990) imply that
$\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}$
converges in distribution to a point mass at zero. This is because the null
model is the closest model (in the Kullback–Leibler sense) to
$m^{*}({\boldsymbol{y}}^{*})$. Therefore we also have

$\frac{1}{n}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}\xrightarrow[n\to\infty]{P}0.$

  * •

Finally, consider
$\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}$.
Substituting back the expressions for ${\boldsymbol{d}}_{\boldsymbol{\gamma}}$
and ${\boldsymbol{E}}_{\boldsymbol{\gamma}}$, expanding the quadratic form and
combining similar terms we have:

$\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}=\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}\left[{\boldsymbol{H}}_{\boldsymbol{\gamma}}-{\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right]\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\\\
+\frac{2}{n}\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}-\frac{1}{n^{2}}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}$

Now, from the Woodbury matrix identity we have

$\left[{\boldsymbol{H}}_{\boldsymbol{\gamma}}-{\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\right]=\left({\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}+n{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*-1}\right)^{-1}$

and therefore

$\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}-{\boldsymbol{d}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{E}}^{-1}_{\boldsymbol{\gamma}}{\boldsymbol{d}}_{\boldsymbol{\gamma}}=\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}^{-1}_{\boldsymbol{\gamma}}+n{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*-1}\right)^{-1}\hat{\boldsymbol{\beta}}_{\boldsymbol{\gamma}}\\\
+\frac{2}{n}\hat{\boldsymbol{\beta}}^{T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}-\frac{1}{n^{2}}\hat{\boldsymbol{\beta}}^{*T}_{\boldsymbol{\gamma}}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\left({\boldsymbol{H}}_{\boldsymbol{\gamma}}+\frac{1}{n}{\boldsymbol{H}}_{\boldsymbol{\gamma}}^{*}\right)^{-1}{\boldsymbol{H}}^{*}_{\boldsymbol{\gamma}}\hat{\boldsymbol{\beta}}^{*}_{\boldsymbol{\gamma}}$

Note that the last two terms converge in probability to $0$, while the first
term is $\mathcal{O}_{P}\left(1\right)$.

##  Appendix E Details of the Markov chain Monte Carlo algorithm for logistic
regression

Using the hierarchical representation of the LPEP prior discussed at the start
of Section 4.2, the posterior distribution for the augmented model can be
written as

$\displaystyle\pi(\boldsymbol{\gamma},\boldsymbol{\beta}_{\boldsymbol{\gamma}},\boldsymbol{\omega},\boldsymbol{y}^{*},\delta\mid\boldsymbol{y},\boldsymbol{X})\propto
f_{\boldsymbol{\gamma}}(\boldsymbol{y}\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},\boldsymbol{\omega})\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})f(\delta|\boldsymbol{\gamma})f(\boldsymbol{\gamma})f(\boldsymbol{\omega})m^{*}(\boldsymbol{y}^{*}).$

From this, it is easy to devise samplers for the full conditional posterior
distributions of various blocks of parameters. We focus below on the more
general setting where $\delta$ has been assigned a hyperprior. The
simplifications for the case where $\delta$ is fixed are straightforward and
we do not discuss them explicitly.

  1. 1.

Since the prior support of $\delta$ maybe dependent on model indicator
$\boldsymbol{\gamma}$, the parameters
$(\boldsymbol{\gamma},\delta,\boldsymbol{\beta}_{\boldsymbol{\gamma}})$ are
updated jointly by sampling from
$f(\boldsymbol{\gamma},\delta,\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})$
given in (9). To do this, we write

$f(\boldsymbol{\gamma},\delta,\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})=f(\boldsymbol{\gamma},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})f(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\boldsymbol{\gamma},\delta,{\boldsymbol{y}}^{*},\boldsymbol{\omega}),$

where
$f(\boldsymbol{\gamma},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})\propto\int
f(\boldsymbol{\gamma},\delta,\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})d\boldsymbol{\beta}_{\boldsymbol{\gamma}}$.

The expression for the conditional posterior
$f(\boldsymbol{\gamma},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})$,
up to a proportionality constant, is given by (10). To generate samples from
it, we generate proposal by combining a random walk Metropolis-Hastings
algorithm for $\boldsymbol{\gamma}$ (George & McCulloch, 1997) and a
reflective Gaussian random walk for $\delta$ (similar to section 2.1 of
Thawornwattana et al., 2018). More specifically, we factorize the joint
proposal for $(\boldsymbol{\gamma},\delta)$ as:

$\displaystyle
q\left(\delta^{(prop)},{\boldsymbol{\gamma}}^{(prop)}\mid\delta,\boldsymbol{\gamma}\right)=q\left({\boldsymbol{\gamma}}^{(prop)}\mid\boldsymbol{\gamma}\right)q\left(\delta^{(prop)}\mid\delta,{\boldsymbol{\gamma}}^{(prop)}\right).$

For $q({\boldsymbol{\gamma}}^{(prop)}\mid\boldsymbol{\gamma})$, we use a
symmetric random walk proposal similar to equation (46) of George & McCulloch
(1997) as follows:

     * •

We define two probability vectors $p_{1}=(0.9,0.1)$ and
$p_{2}=(0.6,0.2,0.15,0.05)$.

     * •

Each time, we decide on one of two types of moves according to the probability
vector $p_{1}$.

       * –

If a move type 1 is selected, then the proposed new model
$\boldsymbol{\gamma}^{(prop)}$ is generated by randomly flipping
$d\in\\{1,2,3,4\\}$ components of $\boldsymbol{\gamma}$ with probability
$p_{2,d}$. The components of $\boldsymbol{\gamma}$ to be flipped are selected
uniformly at random given $d$.

       * –

If a move type 2 is selected, then the proposed model
$\boldsymbol{\gamma}^{(prop)}$ is generated by removing one variable currently
included in the model and replacing it with a variable that is currently
excluded, leaving the dimensionality of the model unchanged. The variables to
add and remove are chosen uniformly at random within each set.

Next, given $\boldsymbol{\gamma}^{(prop)}$, we propose $\delta$ using a
reflective Gaussian random walk with a left reflection boundary
$a_{\boldsymbol{\gamma}^{(prop)}}$. More specifically, we define
$\delta^{(prop)}=a_{\boldsymbol{\gamma}^{(prop)}}+|\epsilon-
a_{\boldsymbol{\gamma}^{(prop)}}|$ where $\epsilon\sim N(\delta,\tau^{2})$ and
$a_{\boldsymbol{\gamma}}=0$ under the hyper-g/n prior and
$a_{\boldsymbol{\gamma}}=\frac{n-p_{\boldsymbol{\gamma}}}{p_{\boldsymbol{\gamma}}+1}$
under the robust prior. In both cases we found $\tau=n/2$ to be an efficient
tuning parameter in our studies.

Since the proposal distribution of $\boldsymbol{\gamma}$, given by
$q({\boldsymbol{\gamma}}^{(prop)}\mid\boldsymbol{\gamma})$ is symmetric, the
proposed model $(\boldsymbol{\gamma}^{(prop)},\delta^{(prop)})$ is then
accepted with probability

$\min\left\\{\frac{f(\boldsymbol{\gamma}^{(prop)},\delta^{(prop)}\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})}{f(\boldsymbol{\gamma},\delta\mid{\boldsymbol{y}}^{*},\boldsymbol{\omega})}\frac{q(\delta\mid\delta^{(prop)},{\boldsymbol{\gamma}})}{q(\delta^{(prop)}\mid\delta,{\boldsymbol{\gamma}}^{(prop)})},1\right\\},$

where

$q\left(\delta^{(prop)}\mid\delta,{\boldsymbol{\gamma}}^{(prop)}\right)=\frac{1}{\sqrt{2\pi}\tau}\Bigg{[}\exp\Bigg{\\{}-\frac{1}{2\tau^{2}}(\delta^{(prop)}-\delta)^{2}\Bigg{\\}}+\\\
\exp\Bigg{\\{}-\frac{1}{2\tau^{2}}(2a_{{\boldsymbol{\gamma}}^{(prop)}}-\delta^{(prop)}-\delta)^{2}\Bigg{\\}}\Bigg{]},\quad\quad\delta^{(prop)}\geq
a_{\boldsymbol{\gamma}^{(prop)}}.$

Note that, when move type 2 is selected, since $p_{\boldsymbol{\gamma}}$
remains unchanged,
$q(\delta\mid\delta^{(prop)},{\boldsymbol{\gamma}})=q(\delta^{(prop)}\mid\delta,{\boldsymbol{\gamma}}^{(prop)})$
and the acceptance probability simplifies further. Similar simplification is
observed under hyper-g/n prior since proposal reflection boundary,
$a_{\boldsymbol{\gamma}}=0$.

Once the model $(\boldsymbol{\gamma},\delta)$ has been sampled, the regression
coefficients can be updated using the fact that
$\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\boldsymbol{\gamma},\delta,{\boldsymbol{y}}^{*},\boldsymbol{\omega}\sim\mathsf{N}\left({\boldsymbol{m}}_{\boldsymbol{\omega}},{\boldsymbol{V}}_{\boldsymbol{\omega}}\right)$,
where ${\boldsymbol{m}}_{\boldsymbol{\omega}}$ and
${\boldsymbol{V}}_{\boldsymbol{\omega}}$ are given in (11).

  2. 2.

A posteriori, the entries of $\boldsymbol{\omega}$ are conditionally
independent from each other. Following Polson et al. (2013), it is
straightforward to see that
$\omega_{i}\mid\boldsymbol{\gamma},\boldsymbol{\beta}_{\boldsymbol{\gamma}},\delta,{\boldsymbol{y}}^{*}\sim
PG(1,\boldsymbol{x}_{i,\boldsymbol{\gamma}}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}})$.
Implementations of the samplers for the Pòlya-Gamma distribution are
available, for example, in the R package BayesLogit.

  3. 3.

The conditional distribution of $\boldsymbol{y}^{*}$ is proportional to

$\displaystyle\pi(\boldsymbol{y}^{*}\mid\delta,\boldsymbol{\beta}_{\boldsymbol{\gamma}})\propto\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})m^{*}(\boldsymbol{y}^{*})$

While this distribution is supported over the finite set $\\{0,1\\}^{n}$, a
direct sampler is difficult to construct in part because of its (typically)
large size of the support. Hence, we rely again on Metropolis-Hastings steps.

In order to ensure adequate mixing of the algorithm, we consider both local
and global proposals. At each iteration, the algorithm selects local moves
with probability $0.7$ and global moves with probability $0.3$.

     * •

For the local moves, we propose new ${\boldsymbol{y}}^{*(prop)}$ by randomly
flipping $d\in\\{1,2,3,4,5\\}$ components of ${\boldsymbol{y}}^{*}$ with
probability $(0.5,0.2,0.15,0.10,0.05)$. The components of
${\boldsymbol{y}}^{*}$ to be flipped are selected uniformly at random given
$d$. Because this proposal is symmetric, the acceptance probability for this
move is simply

$\min\left\\{1,\frac{\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*(prop)})m^{*}(\boldsymbol{y}^{*(prop)})}{\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})m^{*}(\boldsymbol{y}^{*})}\right\\}$

     * •

For the global moves, we use an independent proposal similar to that used by
Fouskakis et al. (2018),
$q(\boldsymbol{y}^{*})=\prod_{i=1}^{n}\text{Berl}({\pi_{i}}^{1/\delta})$ where

$\displaystyle\pi_{i}^{*}$
$\displaystyle=\frac{\pi_{0}^{1/n}\pi_{i,\boldsymbol{\gamma}_{-1}}^{1/\delta}}{\pi_{0}^{1/n}\pi_{i,\boldsymbol{\gamma}_{-1}}^{1/\delta}+(1-\pi_{0})^{1/n}(1-\pi_{i,\boldsymbol{\gamma}_{-1}})^{1/\delta}},$

$\pi_{0}=\frac{1}{1+\exp{(-\beta_{0})}}$,
$\pi_{i,\boldsymbol{\gamma}_{-1}}^{*}=\frac{1}{1+\exp{(-\boldsymbol{x}_{i,{\boldsymbol{\gamma}}_{-1}}^{T}\boldsymbol{\beta}_{\boldsymbol{\gamma}_{-1}}})}$,
and $\boldsymbol{\beta}_{\boldsymbol{\gamma}_{-1}}$ and
$\boldsymbol{\gamma}_{-1}$ represent the coefficient vector and indicator
variable excluding the intercept term. The associated acceptance probability
is then

$\displaystyle\min\left\\{1,\frac{\pi_{\boldsymbol{\gamma}}^{PEP}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,{\boldsymbol{y}^{*(prop)}}\right)m^{*}\left({\boldsymbol{y}^{*(prop)}}\right)q(\boldsymbol{y}^{*})}{\pi_{\boldsymbol{\gamma}}^{PEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})m^{*}(\boldsymbol{y}^{*})q\left({\boldsymbol{y}^{*(prop)}}\right)}\right\\}$

  4. 4.

The fact that $\delta$ is jointly sampled with $\boldsymbol{\gamma}$ in step 1
above means that the algorithm might be slow to mix. In order to address this
issue, we incorporate an additional sampler for $\delta$ alone. The target
full conditional distribution is given by:

$\displaystyle\pi(\delta\mid\boldsymbol{\beta}_{\boldsymbol{\gamma}},\boldsymbol{y}^{*},\boldsymbol{\gamma})\propto\pi_{\boldsymbol{\gamma}}^{PEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})f(\delta|\boldsymbol{\gamma})$

For the prior distributions we discuss in this paper, this full posterior
conditional distribution does not belong to a known family. Hence, we again
use a Metropolis-Hastings algorithm to sample $\delta$ that mimics what we did
in step 1. In particular, we propose new values for $\delta$ from a reflective
Gaussian distribution centered around the current value of $\delta$ and with
scale $\tau=n/2$ and a left reflective boundary $a_{\boldsymbol{\gamma}}=0$ fr
the hyper-g/n prior and
$a_{\boldsymbol{\gamma}}=\frac{n-p_{\boldsymbol{\gamma}}}{p_{\boldsymbol{\gamma}}+1}$
for the robust prior. The proposed values are then accepted with probability:

$\displaystyle\min\left\\{1,\frac{\pi_{\boldsymbol{\gamma}}^{LPEP}\left(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta^{(prop)},\boldsymbol{y}^{*}\right)p\left(\delta^{(prop)}\right)}{\pi_{\boldsymbol{\gamma}}^{LPEP}(\boldsymbol{\beta}_{\boldsymbol{\gamma}}\mid\delta,\boldsymbol{y}^{*})p(\delta)}\frac{q(\delta\mid\delta^{(prop)},{\boldsymbol{\gamma}})}{q(\delta^{(prop)}\mid\delta,{\boldsymbol{\gamma}})}\right\\}.$

## References

  * Agresti (2015) Agresti, A. (2015).  _Foundations of linear and generalized linear models_.  John Wiley & Sons.
  * Albert & Anderson (1984) Albert, A. & Anderson, J. A. (1984).  On the existence of maximum likelihood estimates in logistic regression models.  _Biometrika_ 71, 1–10.
  * Albert & Chib (1993) Albert, J. H. & Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  _Journal of the American statistical Association_ 88, 669–679.
  * Andrade & O’Hagan (2006) Andrade, J. A. A. & O’Hagan, A. (2006).  Bayesian robustness modeling using regularly varying distributions.  _Bayesian Analysis_ 1, 169–188.
  * Andrade & O’Hagan (2011) Andrade, J. A. A. & O’Hagan, A. (2011).  Bayesian robustness modelling of location and scale parameters.  _Scandinavian Journal of Statistics_ 38, 691–711.
  * Bai & Ghosh (2018) Bai, R. & Ghosh, M. (2018).  On the beta prime prior for scale parameters in high-dimensional Bayesian regression models.  _arXiv preprint arXiv:1807.06539_ .
  * Barber et al. (2016) Barber, R. F., Drton, M. & Tan, K. M. (2016).  Laplace approximation in high-dimensional Bayesian regression.  In _Statistical Analysis for High-Dimensional Data_ , pp. 15–36. Springer.
  * Bayarri et al. (2012) Bayarri, M. J., Berger, J. O., Forte, A., García-Donato, G. et al. (2012).  Criteria for Bayesian model choice with application to variable selection.  _The Annals of statistics_ 40, 1550–1577.
  * Berger & Pericchi (1996a) Berger, J. O. & Pericchi, L. R. (1996a).  The intrinsic Bayes factor for linear models.  In _Bayesian Statistics 5_ , Eds. A. P. D. J. M. Bernardo, J. O. Berger & A. F. M. Smith, pp. 25–44. Oxford Univ. Press.
  * Berger & Pericchi (1996b) Berger, J. O. & Pericchi, L. R. (1996b).  The intrinsic Bayes factor for model selection and prediction.  _Journal of the American Statistical Association_ 91, 109–122.
  * Berger et al. (2001) Berger, J. O., Pericchi, L. R., Ghosh, J., Samanta, T., De Santis, F., Berger, J. & Pericchi, L. (2001).  Objective Bayesian methods for model selection: Introduction and comparison.  _Lecture Notes-Monograph Series_ pp. 135–207.
  * Berger et al. (1998) Berger, J. O., Pericchi, L. R. & Varshavsky, J. A. (1998).  Bayes factors and marginal distributions in invariant situations.  _Sankhyā: The Indian Journal of Statistics, Series A_ pp. 307–321.
  * Bhattacharya et al. (2015) Bhattacharya, A., Pati, D., Pillai, N. S. & Dunson, D. B. (2015).  Dirichlet–Laplace priors for optimal shrinkage.  _Journal of the American Statistical Association_ 110, 1479–1490.
  * Bové et al. (2011) Bové, D. S., Held, L. et al. (2011).  Hyper-$g$ priors for generalized linear models.  _Bayesian Analysis_ 6, 387–410.
  * Breheny & Huang (2011) Breheny, P. & Huang, J. (2011).  Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection.  _Annals of Applied Statistics_ 5, 232–253.
  * Brown & Griffin (2010) Brown, P. J. & Griffin, J. E. (2010).  Inference with normal-gamma prior distributions in regression problems.  _Bayesian analysis_ 5, 171–188.
  * Califf et al. (1996) Califf, R. M., White, H. D., Van de Werf, F., Sadowski, Z., Armstrong, P. W., Vahanian, A., Simoons, M. L., Simes, R. J., Lee, K. L. & Topol, E. J. (1996).  One-year results from the Global Utilization of Streptokinase and TPA for Occluded Coronary Arteries (GUSTO-I) trial.  _Circulation_ 94, 1233–1238.
  * Carvalho et al. (2010) Carvalho, C. M., Polson, N. G. & Scott, J. G. (2010).  The horseshoe estimator for sparse signals.  _Biometrika_ 97, 465–480.
  * Chen et al. (2000) Chen, M.-H., Ibrahim, J. G. & Shao, Q.-M. (2000).  Power prior distributions for generalized linear models.  _Journal of Statistical Planning and Inference_ 84, 121–137.
  * Clyde (2020) Clyde, M. (2020).  _BAS: Bayesian Variable Selection and Model Averaging using Bayesian Adaptive Sampling_.  R package version 1.5.5.
  * Cragg et al. (2013) Cragg, M. I., Zhou, Y., Gurney, K. & Kahn, M. E. (2013).  Carbon geography: the political economy of congressional support for legislation intended to mitigate greenhouse gas production.  _Economic Inquiry_ 51, 1640–1650.
  * Dawid (1973) Dawid, A. P. (1973).  Posterior expectations for large observations.  _Biometrika_ 60, 664–667.
  * De Santis & Spezzaferri (2001) De Santis, F. & Spezzaferri, F. (2001).  Consistent fractional bayes factor for nested normal linear models.  _Journal of statistical planning and inference_ 97, 305–321.
  * Dellaportas et al. (2002) Dellaportas, P., Forster, J. J. & Ntzoufras, I. (2002).  On Bayesian model and variable selection using MCMC.  _Statistics and Computing_ 12, 27–36.
  * Denti et al. (2021) Denti, F., Azevedo, R., Lo, C., Wheeler, D., Gandhi, S. P., Guindani, M. & Shahbaba, B. (2021).  A horseshoe pit mixture model for Bayesian screening with an application to light sheet fluorescence microscopy in brain imaging.  _arXiv preprint arXiv:2106.08281_ .
  * Facchini & Steinhardt (2011) Facchini, G. & Steinhardt, M. F. (2011).  What drives us immigration policy? Evidence from congressional roll call votes.  _Journal of Public Economics_ 95, 734–743.
  * Fahrmeir (1990) Fahrmeir, L. (1990).  Maximum likelihood estimation in misspecified generalized linear models.  _Statistics_ 21, 487–502.
  * Fahrmeir & Kaufmann (1985) Fahrmeir, L. & Kaufmann, H. (1985).  Consistency and asymptotic normality of the maximum likelihood estimator in generalized linear models.  _The Annals of Statistics_ 13, 342–368.
  * Fan & Li (2001) Fan, J. & Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  _Journal of the American statistical Association_ 96, 1348–1360.
  * Fouskakis et al. (2015) Fouskakis, D., Ntzoufras, I., Draper, D. et al. (2015).  Power-expected-posterior priors for variable selection in Gaussian linear models.  _Bayesian Analysis_ 10, 75–107.
  * Fouskakis et al. (2018) Fouskakis, D., Ntzoufras, I., Perrakis, K. et al. (2018).  Power-expected-posterior priors for generalized linear models.  _Bayesian Analysis_ 13, 721–748.
  * Friedman et al. (2010) Friedman, J., Hastie, T. & Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  _Journal of Statistical Software_ 33, 1.
  * Frühwirth-Schnatter et al. (2009) Frühwirth-Schnatter, S., Frühwirth, R., Held, L. & Rue, H. (2009).  Improved auxiliary mixture sampling for hierarchical models of non-gaussian data.  _Statistics and Computing_ 19, 479–492.
  * George & McCulloch (1997) George, E. I. & McCulloch, R. E. (1997).  Approaches for Bayesian variable selection.  _Statistica sinica_ pp. 339–373.
  * Geweke (1996) Geweke, J. (1996).  Variable selection and model comparison in regression.  _In Bayesian Statistics 5_ .
  * Green (1995) Green, P. J. (1995).  Reversible jump Markov chain Monte Carlo computation and Bayesian model determination.  _Biometrika_ 82, 711–732.
  * Haughton (1988) Haughton, D. M. (1988).  On the choice of a model to fit data from an exponential family.  _The annals of statistics_ pp. 342–355.
  * Heinze & Schemper (2002) Heinze, G. & Schemper, M. (2002).  A solution to the problem of separation in logistic regression.  _Statistics in medicine_ 21, 2409–2419.
  * Held et al. (2015) Held, L., Bové, D. S. & Gravestock, I. (2015).  Approximate Bayesian model selection with the deviance statistic.  _Statistical Science_ pp. 242–257.
  * Ibrahim et al. (2000) Ibrahim, J. G., Chen, M.-H. et al. (2000).  Power prior distributions for regression models.  _Statistical Science_ 15, 46–60.
  * Jessee & Theriault (2012) Jessee, S. A. & Theriault, S. M. (2012).  The two faces of congressional roll-call voting.  _Party Politics_ p. 1354068812458612.
  * Johnson & Rossell (2010) Johnson, V. E. & Rossell, D. (2010).  On the use of non-local prior densities in Bayesian hypothesis tests.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 72, 143–170.
  * Johnson & Rossell (2012) Johnson, V. E. & Rossell, D. (2012).  Bayesian model selection in high-dimensional settings.  _Journal of the American Statistical Association_ 107, 649–660.
  * Kass & Wasserman (1995) Kass, R. E. & Wasserman, L. (1995).  A reference Bayesian test for nested hypotheses and its relationship to the schwarz criterion.  _Journal of the american statistical association_ 90, 928–934.
  * Konis (2007) Konis, K. (2007).  _Linear Programming Algorithms for Detecting Separated Data in Binary Logistic Regression Models_.  Ph.D. thesis, University of Oxford.
  * Kosmidis & Schumacher (2020) Kosmidis, I. & Schumacher, D. (2020).  _detectseparation: Detect and Check for Separation and Infinite Maximum Likelihood Estimates_.  R package version 0.1.
  * Lee et al. (2020) Lee, S. Y., Pati, D. & Mallick, B. K. (2020).  Continuous shrinkage prior revisited: a collapsing behavior and remedy.  _arXiv preprint arXiv:2007.02192_ .
  * Lesaffre & Albert (1989) Lesaffre, E. & Albert, A. (1989).  Partial separation in logistic discrimination.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 51, 109–116.
  * Li & Pati (2017) Li, H. & Pati, D. (2017).  Variable selection using shrinkage priors.  _Computational Statistics & Data Analysis_ 107, 107–119.
  * Li & Clyde (2018) Li, Y. & Clyde, M. A. (2018).  Mixtures of g-priors in generalized linear models.  _Journal of the American Statistical Association_ 113, 1828–1845.
  * Liang et al. (2008) Liang, F., Paulo, R., Molina, G., Clyde, M. A. & Berger, J. O. (2008).  Mixtures of g-priors for Bayesian variable selection.  _Journal of the American Statistical Association_ 103, 410–423.
  * Lofland et al. (2017) Lofland, C. L., Rodríguez, A. & Moser, S. (2017).  Assessing differences in legislators’ revealed preferences: A case study on the 107th U.S. Senate.  _The Annals of Applied Statistics_ 11, 456–479.
  * Mäkeläinen et al. (1981) Mäkeläinen, T., Schmidt, K. & Styan, G. P. (1981).  On the existence and uniqueness of the maximum likelihood estimate of a vector-valued parameter in fixed-size samples.  _The Annals of Statistics_ pp. 758–767.
  * Mansournia et al. (2018) Mansournia, M. A., Geroldinger, A., Greenland, S. & Heinze, G. (2018).  Separation in logistic regression: causes, consequences, and control.  _American journal of epidemiology_ 187, 864–870.
  * McCullagh & Nelder (2019) McCullagh, P. & Nelder, J. A. (2019).  _Generalized linear models_.  Routledge.
  * Miquel & Snyder Jr (2006) Miquel, G. P. I. & Snyder Jr, J. M. (2006).  Legislative effectiveness and legislative careers.  _Legislative Studies Quarterly_ 31, 347–381.
  * Mitchell & Beauchamp (1988) Mitchell, T. J. & Beauchamp, J. J. (1988).  Bayesian variable selection in linear regression.  _Journal of the american statistical association_ 83, 1023–1032.
  * Moser et al. (2021) Moser, S., Rodríguez, A. & Lofland, C. L. (2021).  Multiple ideal points: Revealed preferences in different domains.  _Political Analysis_ 29, 139–166.
  * O’Hagan (1995) O’Hagan, A. (1995).  Fractional Bayes factors for model comparison.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 57, 99–118.
  * O’Roark & Wood (2011) O’Roark, J. B. & Wood, W. C. (2011).  Determinants of congressional minimum wage support: the role of economic education.  _Public Choice_ 147, 209–225.
  * Park & Casella (2008) Park, T. & Casella, G. (2008).  The Bayesian Lasso.  _Journal of the American Statistical Association_ 103, 681–686.
  * Pérez & Berger (2002) Pérez, J. M. & Berger, J. O. (2002).  Expected-posterior prior distributions for model selection.  _Biometrika_ 89, 491–512.
  * Polson & Scott (2012) Polson, N. G. & Scott, J. G. (2012).  Local shrinkage rules, Lévy processes and regularized regression.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 74, 287–311.
  * Polson et al. (2013) Polson, N. G., Scott, J. G. & Windle, J. (2013).  Bayesian inference for logistic models using Pólya–Gamma latent variables.  _Journal of the American statistical Association_ 108, 1339–1349.
  * Potter (2005) Potter, D. M. (2005).  A permutation test for inference in logistic regression with small-and moderate-sized data sets.  _Statistics in medicine_ 24, 693–708.
  * Schwarz (1978) Schwarz, G. (1978).  Estimating the dimension of a model.  _The annals of statistics_ pp. 461–464.
  * Steyerberg et al. (2019) Steyerberg, E. W. et al. (2019).  _Clinical prediction models_.  Springer.
  * Thawornwattana et al. (2018) Thawornwattana, Y., Dalquen, D. & Yang, Z. (2018).  Designing simple and efficient Markov chain Monte Carlo proposal kernels.  _Bayesian Analysis_ 13, 1037–1063.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the Lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 58, 267–288.