# A hierarchical prior for generalized linear models based on predictions for
the mean response

Ethan M. Alt
Department of Biostatistics
University of North Carolina
Chapel Hill, NC 27516
<EMAIL_ADDRESS>
&Matthew A. Psioda
Department of Biostatistics
University of North Carolina
Chapel Hill, NC 27516
<EMAIL_ADDRESS>
Joseph G. Ibrahim
Department of Biostatistics
University of North Carolina
Chapel Hill, NC 27516
<EMAIL_ADDRESS>

###### Abstract

There has been increased interest in using prior information in statistical
analyses. For example, in rare diseases, it can be difficult to establish
treatment efficacy based solely on data from a prospective study due to low
sample sizes. To overcome this issue, an informative prior for the treatment
effect may be elicited. We develop a novel extension of the conjugate prior of
Chen and Ibrahim (2003) that enables practitioners to elicit a prior
prediction for the mean response for generalized linear models, treating the
prediction as random. We refer to the hierarchical prior as the hierarchical
prediction prior. For i.i.d. settings and the normal linear model, we derive
cases for which the hyperprior is a conjugate prior. We also develop an
extension of the HPP in situations where summary statistics from a previous
study are available, drawing comparisons with the power prior. The HPP allows
for discounting based on the quality of individual level predictions, having
the potential to provide efficiency gains (e.g., lower MSE) where predictions
are incompatible with the data. An efficient Markov chain Monte Carlo
algorithm is developed. Applications illustrate that inferences under the HPP
are more robust to prior-data conflict compared to selected non-hierarchical
priors.

_Keywords_ hierarchical models $\cdot$ generalized linear models $\cdot$
bayesian inference $\cdot$ hyperprior

##  1 Introduction

Exponential family models, which include distributions for binary, count, and
continuous data, are among the most utilized models for statistical analysis.
In many application areas, it is desirable to incorporate prior information in
the analysis of data. In the Bayesian paradigm, such information may be
incorporated through an informative prior distribution on the parameters of
interest and, possibly, nuisance parameters. For example, in rare disease
clinical trials, it can be difficult to establish treatment efficacy based
solely on data from a prospective study. In ecological models, informative
priors on the basis of similar studies, historical data, or expert opinion are
often utilized (see, e.g., Kuhnert et al. (2010) and Banner et al. (2020) and
references therein). In applied econometric models, informative priors for
generalized linear models (GLMs) have been, for example, utilized to predict
unemployment rates (Grzenda, 2016) and to model the probability of credit
default (Wang et al., 2018).

When previous studies have been conducted, it is often desirable to construct
a prior based on these data. Three popular priors that have been proposed for
this setting are the power prior (Ibrahim and Chen, 2000), commensurate priors
(Hobbs et al., 2012), and meta-analytic-predictive (MAP) priors (e.g.,
Schmidli et al. (2014)). The primary limitation of the former two priors is
that they require access to an entire historical set in order to implement a
joint prior for the regression coefficients that are correlated a priori. MAP
priors are limited in that they have only been developed under the context
where there is a single parameter of interest. In many data applications, it
is desirable to elicit a joint prior for regression coefficients. For example,
in clinical trials and observational studies, it is often of interest to
determine whether an intervention is more/less efficacious for certain groups
(i.e., effect modification). Another example arises when prediction is of
paramount interest, such as in applied econometrics.

A multivariate normal prior on the parameters of interest and nuisance
parameters may be utilized using summary statistics from publications for
previously completed studies (e.g., as done in Rodhouse et al. (2019)).
However, in regression models, this induces an a priori independent prior on
the regression coefficients, which is primarily a function of constraints on
the information source rather than any a priori belief supporting
independence. In contrast, a posteriori, the posterior distribution for
regression coefficients almost always exhibits dependence.

Moreover, when a multivariate normal prior is elicited for a parameter on the
basis of expert opinion, it may be difficult to translate the information
provided to the regression parameter scale, and also to appropriately
translate uncertainty. For example, in logistic regression models, a prior on
a treatment effect would constitute a prior on the logarithm of the odds
ratio. While one may transform an expert’s opinion about the response
probability to a log odds ratio, it is less clear how to elicit the variance
on this prior. A second type of prior for exponential family models involves
the elicitation of a prior prediction for the mean response (e.g., Diaconis
and Ylvisaker (1979) for i.i.d. settings and Chen and Ibrahim (2003) for
regression settings). This type of prior may be more natural for expert
elicitation than the previous, as the prior is elicited on the same scale as
the data. However, these priors do not directly account for uncertainty in the
prediction. Specifically, uncertainty may only be elicited through discounting
of the informativeness of the prior.

In this paper, we develop a novel extension of the conjugate prior of Chen and
Ibrahim (2003), where the prediction of the mean response is treated as
random. Treatment of the prior prediction as random is intuitive because the
prior prediction will typically made on the basis of summary statistics or
expert opinion, both of which have some degree of uncertainty. We refer to the
hierarchical prior as the hierarchical prediction prior (HPP). In regression
models, the HPP induces a correlation structure on the regression coefficients
a priori. The conjugate priors of Diaconis and Ylvisaker (1979) and Chen and
Ibrahim (2003) may be cast as special cases of the HPP for i.i.d. and
regression models, respectively. Moreover, the HPP is quite flexible, enabling
practitioners to elicit predictions on the mean response based on one or more
covariates.

The posterior density under the HPP is robust to incompatible prior
predictions for the mean response. In particular, we show for i.i.d. settings
and the normal linear model that the posterior means of the predictions are
between the predicted values based on the maximum likelihood estimator (MLE)
and the elicited prior prediction. While a formal proof cannot be offered for
all GLMs, the result held for a Poisson regression example. Moreover, under
some limiting cases, we show for i.i.d. data settings and the normal linear
model that the HPP is a conjugate prior for both the parameters of interest
and the hierarchical parameters.

We compare the posterior distribution under the HPP for a logistic regression
model using the data of Finney (1947) against that of the conjugate prior of
Chen and Ibrahim (2003), which we henceforth refer to as the CI prior. The
results suggest that the posterior distribution of the regression coefficients
under the CI prior is more sensitive to the prior prediction of the mean
response than that for the HPP. We also illustrate how to utilize the HPP for
settings in which a previous study was conducted, but only the point estimates
and standard errors for the regression coefficients are available, comparing
the posterior distribution of the parameter of interest under the HPP against
that for the CI prior, an approximation to the asymptotic power prior (Ibrahim
et al., 2015), and the power prior (Ibrahim and Chen, 2000), the latter for
which the entire historical data set was utilized. When the historical data
set is incompatible with the current data set (e.g., when the historical data
provide evidence that treatment is beneficial and the current data suggest
treatment is unbeneficial), the posterior of the treatment effect under the
HPP consistently placed more mass at 0, while the posterior distributions for
the treatment effect under the other priors suggested that treatment was all
but certain to be efficacious.

In i.i.d. settings, samples from the posterior density under the HPP is
straightforward utilizing any Markov-chain Monte Carlo (MCMC) algorithm. For
regression settings excluding the normal linear model, sampling is more
involved because there is no general closed-form of the normalizing constant
for the CI prior. We utilize a fast and accurate Laplace approximation to the
normalizing constant, which enables efficient sampling using MCMC methods.

The remainder of the article is organized as follows. In Section 2, we review
the conjugate prior of Diaconis and Ylvisaker (1979) and develop the HPP for
independent and identically distributed (i.i.d.) data settings. In Section 3
we develop the HPP for generalized linear models (GLMs). In Section 4, we
illustrate how the HPP may be utilized for settings in which summary
statistics are available (e.g., in a publication). In Section 5, we compare
results from an analysis of the data of Finney (1947) using the HPP and the CI
prior for settings in which the prior is formulated on the basis of expert
opinion. In Section 6, we compare posterior quantities under the HPP with
other priors based on data from a previously conducted study. In Section 7, we
close with some discussion.

##  2 The i.i.d. case

We begin by considering the independently and identically distributed (i.i.d.)
data setting for exponential family models. We develop the HPP for the i.i.d.
case and establish theoretical connections to the conjugate prior of Diaconis
and Ylvisaker (1979).

###  2.1 Hyperprior motivation and construction

Here, we discuss the motivation for the HPP and some of its properties.
Suppose we observe $\\{y_{i},i=1,\ldots,n\\}$ with likelihood function in the
exponential family, given by

$\displaystyle
f(\bm{y}|\theta,\phi)=\prod_{i=1}^{n}\exp\left\\{\frac{1}{a_{i}(\phi)}\left[y_{i}\theta-b(\theta)\right]+c(y_{i};\phi)\right\\},$
(1)

where $a_{i}$ is a positive function and typically $a_{i}(\phi)=w_{i}\phi$,
$\bm{y}=(y_{1},\ldots,y_{n})^{\prime}$, $\theta\in\Theta$ is referred to as
the canonical parameter, where $\Theta$ is the domain of $\theta$, and the
functions $b$ and $c$ index the density or mass function. We assume
$a_{i}(\phi)=\phi$ is known and fixed, and we may suppose without loss of
generality that $\phi=1$. Diaconis and Ylvisaker (1979) showed that each
distribution in the exponential family admits a conjugate prior of the form

$\displaystyle\operatorname{\pi_{\text{DY}}}(\theta|\lambda,m)=\frac{1}{Z(\lambda,\lambda
m)}\exp\left\\{\lambda\left[\theta
m-b(\theta)\right]\right\\},\theta\in\Theta,$ (2)

where $Z(\lambda,\lambda m)=\int_{\Theta}\exp\left\\{\lambda\left[\theta
m-b(\theta)\right]\right\\}d\theta$ is a normalizing constant, $\lambda$ is a
precision parameter typically chosen so that $\lambda\in(0,n]$, and
$m\in\dot{b}(\Theta)$ is a location parameter, where $\dot{b}$ is the first
derivative of the function $b$ and $\dot{b}(\Theta)$ is the image of the set
$\Theta$ under the function $\dot{b}$. Henceforth, we refer to the prior (2)
simply as the “DY prior.” Diaconis and Ylvisaker (1979) showed that, under the
DY prior, $E(y)=m$, that is, $m$ may be interpreted as a prior prediction (or
“guess”) for the mean of $y$. The hyperparameter $\lambda$ controls for the
level of informativeness in the prior. Let
$\bm{y}=(y_{1},\ldots,y_{n})^{\prime}$. The posterior density utilizing the
likelihood (1) and prior (2) is given by

$\displaystyle p(\theta|\bm{y},\lambda,m)$
$\displaystyle=\operatorname{\pi_{\text{DY}}}\left(\theta\left|n+\lambda,\frac{n\bar{y}+\lambda
m}{n+\lambda}\right)\right.,$

where $\bar{y}=n^{-1}\sum_{i=1}^{n}y_{i}$ is the sample mean. The primary
disadvantage of the DY prior is that the posterior is sensitive to the prior
prediction $m$ when $\lambda$ is large. For example, let $\mu=\dot{b}(\theta)$
denote the mean parameter of the distribution of the response variables. It
can be shown that the posterior mean of $\mu$ under the DY prior is given by
$E(\mu|\bm{y},\lambda,m)=(n\bar{y}+\lambda m)/(n+\lambda),$ i.e., the
posterior mean of $\mu$ under the DY prior is a convex combination of the
observed sample mean $\bar{y}$ and the prior prediction $m$, with higher
values of $\lambda$ putting more weight on the value of $m$ in the posterior
mean of $\mu$. Thus, if $m$ is an inaccurate prediction for $E(y)$, the
posterior under the DY prior could give misleading results.

Moreover, the hyperparameter $\lambda$ in the DY prior cannot simultaneously
control the uncertainty surrounding $m$ and the level of borrowing from the
prior. For example, in a Bernoulli model, suppose an expert believes with
$95\%$ probability that the success probability is between 0.2 and 0.4. This
may be accomplished by eliciting $m=0.3$ and $\lambda=78.8$ in the DY prior
(2). This imposes a restriction on the sample size, namely, $n\geq 79$.
Conversely, for fixed values of $n$, $\lambda$, and $m$, we cannot adjust the
DY prior to account for uncertainty surrounding the value of $m$.

Since $m$ is typically elicited on the basis of summary statistics or expert
opinion, both of which are measured with uncertainty, it is natural to view
$m$ as having a probability distribution. Thus, we propose to treat the
hyperparameter $m$ in the DY prior as random. We refer to the joint prior of
$(\theta,m)$ as the hierarchical prediction prior (HPP), and the prior on $m$
simply as the hyperprior.

We now develop the hyperprior. Let $\nu=\dot{b}(m)\in\Theta$. We elicit

$\displaystyle\pi(\nu|\lambda_{0},\mu_{0})\propto\exp\left\\{\lambda_{0}[\nu\mu_{0}-b(\nu)]\right\\},\nu\in\Theta,$
(3)

which we may write in terms of the shape parameter $m$ as

$\displaystyle\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0})\propto\exp\left\\{\lambda_{0}[\dot{b}(m)\mu_{0}-b(\dot{b}(m))]\right\\}\frac{1}{v(m)},m\in\dot{b}(\Theta),$
(4)

where $v(m)=\ddot{b}\circ\dot{b}^{-1}(m)$ is the variance function associated
with the exponential family model (1). For each exponential family model, the
hyperprior is a recognizable density (e.g., a beta density for binomial
models, a gamma density for Poisson models, an inverse-gamma density for gamma
models, and a normal density for normal models). The hyperparameter
$\lambda_{0}$ is a precision parameter controlling for the level of certainty
surrounding the prior prediction.

The HPP is obtained by combining the DY prior (2) and the hyperprior (4),
giving

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\theta,m|\lambda,\lambda_{0}\mu_{0})=\operatorname{\pi_{\text{CI}}}(\theta|\lambda,m)\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0}),\theta\in\Theta,m\in\dot{b}(\Theta).$
(5)

The HPP is thus the CI prior conditional on $m$ multiplied by a density over
$m$. It can be shown that the prior mean of $m$ in the hyperprior is
$E(m)=\mu_{0}$. Thus,

$E(y)=E_{\theta}[E(y|\theta)]=E_{\theta}[\dot{b}(\theta)]=E_{m}[E_{\theta|m}(\dot{b}(\theta))]=E(m)=\mu_{0},$

so that, similar to the DY prior, $\mu_{0}$ may be interpreted as a prior
prediction for $E(y)$. However, unlike the DY prior, the HPP allows
practitioners to directly elicit uncertainty surrounding the prior prediction
$\mu_{0}$ for any fixed level of $\lambda$.

Note that as $\lambda_{0}\to\infty$,
$\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0})$ converges to a point
mass at $\mu_{0}$, so that the HPP converges to the CI prior as
$\lambda_{0}\to\infty$. The marginal prior for $\theta$ under the HPP is given
by

$\operatorname{\pi_{\text{HPP}}}(\theta|\lambda,\lambda_{0},\mu_{0})=\int\operatorname{\pi_{\text{CI}}}(\theta|\lambda,m)\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0})dm,$

i.e., the prior for $\theta$ is the CI prior conditional on $m$ averaged over
a distribution over $m$. In the Bernoulli example above, we may elicit
$\mu_{0}=0.30$ and $\lambda_{0}=78.8$, so that, a priori, $P(0.2\leq m\leq
0.4)=0.95$ for any value of $\lambda$. In general, unlike the DY prior, the
HPP allows practitioners to elicit their uncertainty surrounding the prior
prediction while fixing the amount of information borrowed from the prior.

Combining (1), (2), and (4) yields the joint posterior density

$\displaystyle p(\theta,m|\bm{y},\lambda,\lambda_{0},\mu_{0})$
$\displaystyle\propto\operatorname{\pi_{\text{DY}}}\left(\theta|n+\lambda,\frac{n\bar{y}+\lambda
m}{n+\lambda}\right)\frac{Z(n+\lambda,n\bar{y}+\lambda m)}{Z(\lambda,\lambda
m)}\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0}),$
$\displaystyle\propto\operatorname{\pi_{\text{DY}}}\left(\theta|n+\lambda,\frac{n\bar{y}+\lambda
m}{n+\lambda}\right)p(m|\bm{y},\lambda,\lambda_{0},\mu_{0}).$ (6)

The relationship (6) indicates that the joint posterior density may be
expressed as the posterior using DY prior conditional on $m$ multiplied by the
posterior density of $m$. Note that as $\lambda_{0}\to\infty$,
$p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})$ is simply a point mass at
$m=\mu_{0}$, and the posterior distribution of $\theta$ thus converges to the
posterior under the DY prior.

We may write the marginal posterior distribution of $\theta$ as

$\displaystyle
p(\theta|\bm{y},\lambda,\lambda_{0},\mu_{0})\propto\int\operatorname{\pi_{\text{DY}}}\left(\theta|n+\lambda,\frac{n\bar{y}+\lambda
m}{n+\lambda}\right)p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})dm.$ (7)

The relationship (7) indicates that the marginal posterior density of $\theta$
under the HPP may be interpreted as the posterior utilizing the DY prior with
a fixed value of $m$ averaged over the posterior distribution of $m$. In
Section 2.2 we prove that the posterior mean of $m$ is between $\bar{y}$ and
$\mu_{0}$ as $\lambda_{0}\to\infty$. Simulation results suggest that the
result holds in general. In effect, the marginal posterior distribution of
$\theta$ depends on the data $\bm{y}$ in two ways: through the shape parameter
in the posterior of the DY prior conditional on $m$, and in the marginal
posterior distribution for $m$. Conversely, the posterior density of the DY
prior treats $m=\mu_{0}$ as fixed. If $\mu_{0}$ is a poor prediction for
$E(y)$, then inference on the posterior distribution of $\theta$ could provide
misleading results.

###  2.2 Limiting posterior distributions

As previously stated, the relationship (7) indicates that the marginal
posterior distribution of $\theta$ under the HPP may be arrived at by taking
posterior under the DY prior conditional on $m$ averaged over the posterior
distribution of $m$, which depends on both the prior prediction $\mu_{0}$ as
well as the observed sample mean $\bar{y}$. We may write the posterior
distribution of $m$ as

$\displaystyle
p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})\propto\frac{Z(\lambda+n,\lambda
m+n\bar{y})}{Z(\lambda,\lambda
m)}\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0}),$ (8)

where $Z(a,c)=\int_{\Theta}\exp\\{a[\theta(c/a)-b(\theta)]\\}d\theta$ is the
normalizing constant of the DY prior with precision parameter $c$ and shape
parameter $a$ and $\operatorname{\pi_{\text{HPP}}}(m|\lambda_{0},\mu_{0})$ is
defined in (4). As $\lambda\to\infty$, the hyperprior becomes a conjugate
density, whose posterior mean is a convex combination of the sample mean and
the prior prediction. We state this formally in the following theorem.

Theorem 2.1: Let $y_{1},\ldots,y_{n}$ be observations from the exponential
family (1). Let the prior for $(\theta,m)$ be given by the HPP (5). Then

$p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})\to\operatorname{\pi_{\text{HPP}}}\left(m\left|n+\lambda_{0},\frac{n\bar{y}+\lambda_{0}\mu_{0}}{n+\lambda_{0}}\right)\right.\text{
as }\lambda\to\infty.$

In words, Theorem 2.1 states that as $\lambda\to\infty$, the marginal prior of
$m$ under the HPP is conjugate density. Hence, even as $\lambda\to\infty$, the
posterior distribution of $\bm{m}$ is non-degenerate. Because the marginal
prior for $m$ in the HPP is conjugate in the limit, Theorem 2.1 provides an
approximation for the posterior of $m$ for large values of $\lambda$ in terms
of recognizable densities (e.g., a gamma distribution for Poisson models).

An important consequence of Theorem 2.1 is

$\displaystyle
E(m|\bm{y},\lambda,\lambda_{0},\mu_{0})\to\frac{n\bar{y}+\lambda_{0}\mu_{0}}{n+\lambda_{0}}\text{
as }\lambda\to\infty,$

so that the posterior mean of $m$ converges to a convex combination of the
observed sample mean, $\bar{y}$, and the elicited prior mean, $\mu_{0}$. This
is closely related to, but quite different than, the posterior distribution of
$\mu$ using the DY prior. The parameter $\mu$ is a model parameter, while the
hierarchical parameter $m$ is a prior prediction for $\mu$. We prove Theorem
2.1 for the Poisson, Bernoulli, gamma, and normal models in Section 1 of the
Supplementary Materials.

Corollary 2.1 Given the same setup as Theorem 2.1,

$\lim_{\lambda\to\infty}p(\mu|\bm{y},\lambda,\lambda_{0},\mu_{0})\to\operatorname{\pi_{\text{HPP}}}\left(\mu\left|n+\lambda_{0},\frac{n\bar{y}+\lambda_{0}\mu_{0}}{n+\lambda_{0}}\right)\right..$

Proof of Corollary 2.1: Note that

$p(\mu|\bm{y},\lambda,\lambda_{0},\mu_{0})=\int
p(\mu|m,\bm{y},\lambda,\lambda_{0},\mu_{0})p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})dm.$

As $\lambda\to\infty$, $p(\mu|m,\bm{y},\lambda,\lambda_{0},\mu_{0})$ converges
to a point mass at $m$ and $p(m|\bm{y},\lambda,\lambda_{0},\mu_{0})$ converges
to
$\operatorname{\pi_{\text{HPP}}}(m|n+\lambda_{0},(n\bar{y}+\lambda_{0}\mu_{0})/(n+\lambda_{0}))$
by Theorem 2.1. Hence, the result follows.

In words, Corollary 2.1 states that, as $\lambda\to\infty$, the posterior
distribution for the mean parameter, $\mu$, under the HPP converges to the
posterior density under the DY prior with precision parameter $n+\lambda_{0}$
and mean parameter $(n\bar{y}+\lambda_{0}\mu_{0})/(n+\lambda_{0})$, which is
precisely the same posterior distribution obtained by eliciting
$\lambda=\lambda_{0}$ and $m=\mu_{0}$ in the DY prior (2). When $\lambda\leq
n$ (as is typically the case), the HPP cannot have more influence on the
posterior than the likelihood. This fact and Corollary 2.1 imply that if
$\lambda\to\infty$ and $\lambda_{0}\leq n$ or $\lambda_{0}\to\infty$ and
$\lambda\leq n$, the posterior distribution of $\mu$ cannot depend more on the
prior than the data.

Remark 2.1: By transforming the mean parameter, $\mu$, to the canonical
parameter $\theta$, it is easy to see from Corollary 2.1 that
$p(\theta|\bm{y},\lambda,\lambda_{0},\mu_{0})\to\operatorname{\pi_{\text{DY}}}(\theta|n+\lambda_{0},(n\bar{y}+\lambda_{0}\mu_{0})/(n+\lambda_{0}))$
as $\lambda\to\infty$, therefore, the posterior density of $\theta$ under the
HPP converges to that under the DY prior, with shape parameter equal to a
convex combination of the sample mean and prior prediction.

Remark 2.2: For the i.i.d. normal case, the HPP induces a conjugate prior on
$\theta=\mu$ for finite $\lambda$. The difference between the DY prior and the
HPP for the normal case is that the prior variance of $\mu$ under the HPP is
larger than the DY prior for $\lambda_{0}<\infty$.

##  3 The regression case

In this section, we develop the HPP and illustrate its properties for GLMs. We
derive the posterior distribution of the regression coefficients for the
normal linear model under the HPP, relating it to the posterior under the CI
and multivariate normal priors. Furthermore, we discuss the posterior
distribution of $\bm{m}$ and how possible efficiency gains may be made
achieved even when some of the components of $\bm{\mu}_{0}$ are inaccurate.
Finally, we discuss how to implement the HPP computationally to obtain
posterior samples.

###  3.1 The HPP for GLMs

Throughout this section, suppose we observe
$\\{(y_{i},\bm{x}_{i}),i=1,\ldots,n\\}$, where $y_{i}$ is a response variable
and $\bm{x}_{i}$ is a $p\times 1$ vector of covariates associated with subject
$i$, which may include an intercept term. Suppose $E(y_{i})=\mu_{i}$, where
$g(\mu_{i})=\bm{x}_{i}^{\prime}\bm{\beta}$, where $\bm{\beta}$ is a $p\times
1$ vector of regression coefficients. The function $g$ is referred to as the
$\mu$-link function. The likelihood function of
$\bm{y}=(y_{1},\ldots,y_{n})^{\prime}$ may be written as

$\displaystyle
f(\bm{y}|\phi,\bm{\beta},\bm{X})=\prod_{i=1}^{n}\exp\left\\{\frac{1}{a_{i}(\phi)}\left[y_{i}\theta(\bm{x}_{i}^{\prime}\bm{\beta})-b(\theta(\bm{x}_{i}^{\prime}\bm{\beta}))\right]+c(y_{i},\phi)\right\\},$
(9)

where $\theta(x)=\dot{b}^{-1}\circ g^{-1}(x)$ is referred to as the
$\theta$-link function and the functions $b$ and $c$ index the density or mass
function. We assume $a_{i}(\phi)=\phi$ for $i=1,\ldots,n$ is known and fixed,
and we may suppose without loss of generality that $\phi=1$. Chen and Ibrahim
(2003) showed that the likelihood (9) admits a conjugate prior of the form

$\displaystyle\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})=\frac{1}{Z(\lambda,\lambda\bm{m})}\exp\left\\{\lambda\left[\bm{m}^{\prime}\theta(\bm{X\beta})-\bm{J}^{\prime}b(\theta(\bm{X\beta}))\right]\right\\},\bm{\beta}\in\mathbb{R}^{p},$
(10)

where
$Z(\lambda,\lambda\bm{m})=\int_{\mathbb{R}^{p}}\exp\left\\{\lambda\left[\bm{m}^{\prime}\theta(\bm{X\beta})-\bm{J}^{\prime}b(\theta(\bm{X\beta}))\right]\right\\}d\bm{\beta}$
is the normalizing constant, which has no closed-form expression in general,
$\lambda>0$ is a precision parameter controlling for the informativeness of
the prior, $\bm{m}$ is a $n$-dimensional shape parameter that may be
interpreted as a prior prediction for $E(\bm{y})$, and the function $b$ is
taken componentwise, i.e.,
$b(\bm{\eta})=(b(\eta_{1}),\ldots,b(\eta_{n}))^{\prime}$ for
$\eta\in\mathbb{R}^{n}$. We refer to the prior (10) as the “CI prior.” When
$\bm{m}$ is fixed, Chen and Ibrahim (2003) showed that the posterior is

$\displaystyle p(\bm{\beta}|\lambda,\bm{m})$
$\displaystyle=\operatorname{\pi_{\text{CI}}}\left(\bm{\beta}\left|1+\lambda,\frac{\bm{y}+\lambda\bm{m}}{1+\lambda}\right)\right..$
(11)

Typically, $\lambda\in(0,1]$ so that the effective sample size contributed by
the prior does not exceed that of the data. Note that, similar to the i.i.d.
case, the posterior distribution of $\bm{\beta}$ for fixed $\bm{m}$ has a
shape parameter that is a convex combination of the observed data $\bm{y}$ and
prior prediction $\bm{m}$. Thus, it is clear that the posterior (11) is
sensitive to the choice of $\bm{m}$.

One of the challenges encountered when using the CI prior (10) is that it is
difficult to quantify uncertainty surrounding the value of $\bm{m}$. That is,
the CI prior treats $\bm{m}$ as observed data. While the precision
hyperparameter $\lambda$ controls for the degree of influence of the prior on
the posterior, it does not explicitly reflect uncertainty in the prior
prediction. For example, similar to the i.i.d. case, $\lambda$ cannot
simultaneously control for the uncertainty surrounding $\bm{m}$ as well as the
discounting of the CI prior. Moreover, elicitation of $\bm{m}$ is typically
based on expert opinion or summary statistics, and, thus, it is natural to
view $\bm{m}$ as having a probability distribution.

To this end, we now derive the hyperprior. Let $\bm{\nu}=\dot{b}(\bm{m})$,
where the function $\dot{b}$ is taken component-wise. We may elicit

$\displaystyle\pi(\bm{\nu}|\lambda_{0},\bm{\mu}_{0})\propto\exp\left\\{\lambda_{0}\left[\bm{\nu}^{\prime}\bm{\mu}_{0}-\bm{J}^{\prime}b(\bm{\nu})\right]\right\\}=\prod_{i=1}^{n}\exp\left\\{\lambda_{0}\left[\nu_{i}\mu_{0i}-b(\nu_{i})\right]\right\\},\bm{\nu}\in\Theta^{n},$
(12)

where $\bm{\nu}=\dot{b}(\bm{m})$, $\lambda_{0}>0$ is a precision parameter,
and $\bm{\mu}_{0}\in[\dot{b}(\Theta)]^{n}$ is a $n$-dimensional vector giving
the prior prediction for $E(\bm{y})$, which may depend on covariates. In
general, one may utilize a separate precision parameter for each component of
$\bm{\nu}$, but we will proceed with a common precision parameter for
notational convenience (allowing each component to have its own precision is
discussed in Section 4). Note that (12) is a product of $n$ independent priors
of the form (3). Thus the prior for regression settings (12) can be thought of
as a multivariate generalization of that for i.i.d. settings (3). For example,
if $\bm{\nu}=\nu$ is a scalar, then (12) is equal to (3).

Using (12), the hyperprior on the mean scale is obtained using the
transformation $\bm{m}=\dot{b}^{-1}(\bm{\nu})$, i.e.,

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0})\propto\prod_{i=1}^{n}\exp\left\\{\lambda_{0}\left[\dot{b}^{-1}(m_{i})^{\prime}\mu_{0i}-b(\dot{b}^{-1}(m_{i}))\right]\right\\}\frac{1}{v(m_{i})},\bm{m}\in[\dot{b}(\Theta)]^{n},$
(13)

where $v(x)=\ddot{b}\circ\dot{b}^{-1}(x)$ is the variance function of the
family. Note that the hyperprior (13) is a product of $n$ independent
densities, each having the same form as the i.i.d. case (4), where each
component of $\bm{m}$ has its own mean. The hyperprior is thus a product of
$n$ independent recognizable densities (e.g., beta for binomial models, gamma
for Poisson models, inverse-gamma for gamma models, and normal for normal
models). An attractive feature of the HPP is that, if covariates and
regression coefficients are ignored, the HPP is a conjugate prior for each
component of $\bm{y}$. For example, if $\bm{y}=(y_{1},\ldots,y_{n})^{\prime}$
is a collection of Bernoulli random variables each with mean $\mu_{i}$,
$\\{m_{i},i=1,\ldots,n\\}$ is a priori a collection of $n$ independent beta
random variables with mean $\mu_{0i}$ and dispersion parameter $\lambda_{0}$.

In some cases, it is desirable to elicit a prior on $\bm{m}$ based on
(possibly a subset of) covariates, which we view as fixed. For example, we may
decompose $\bm{\beta}=(\bm{\alpha}^{\prime},\bm{\gamma}^{\prime})^{\prime}$
and $\bm{x}_{i}=(\bm{x}_{1i}^{\prime},\bm{x}_{2i}^{\prime})^{\prime}$. If
$\bm{\alpha}_{0}$ is a prediction for $\bm{\alpha}$, we may take
$\mu_{0i}=g^{-1}(\bm{x}_{1i}^{\prime}\bm{\alpha}_{0})$, where $g$ is the
$\mu$-link function. In this case, we may write (13) as

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0})\propto\exp\left\\{\lambda_{0}\left[\dot{b}^{-1}(\bm{m})^{\prime}g^{-1}(\bm{X}_{1}\bm{\alpha}_{0})-\bm{J}^{\prime}b(\dot{b}^{-1}(\bm{m}))\right]\right\\}.$
(14)

For example, in a Poisson regression model with canonical link function, we
have $b(x)=\exp\\{x\\}$ and $g(\mu)=\exp\\{\mu\\}$, so that we may write the
hyperprior (14) as

$\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\alpha}_{0})\propto\exp\left\\{\lambda_{0}\left[\log(\bm{m})^{\prime}\exp\\{\bm{X}_{1}\bm{\alpha}_{0}\\}-\bm{J}^{\prime}\bm{m}\right]\right\\}.$

Thus, the hyperprior is proportional to the likelihood of a gamma GLM with
known inverse dispersion parameter $\lambda_{0}$, known regression
coefficients $\bm{\alpha}_{0}$, and $\mu$-link function $g(\mu)=\log(\mu)$.
Hence, the hyperprior (14) may be thought of as a random process for
predicting the mean response based on a prediction $\bm{\alpha}_{0}$ based on
covariates $\bm{X}_{1}$.

The HPP for regression models is obtained by combining the CI prior (10) and
the hyperprior (14), giving

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{\beta},\bm{m}|\lambda,\lambda_{0},\bm{\mu}_{0})$
$\displaystyle=\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0}).$
(15)

Thus, the marginal prior for the regression coefficients is

$\operatorname{\pi_{\text{HPP}}}(\bm{\beta}|\lambda,\lambda_{0},\mu_{0})=\int\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0})d\bm{m}.$

The marginal prior of $\bm{\beta}$ under the HPP is thus the CI prior
conditional on $\bm{m}$ averaged over a distribution on $\bm{m}$. As
uncertainty surrounding the prior prediction decreases, i.e., for larger
values of $\lambda_{0}$, the marginal prior of $\bm{\beta}$ under the HPP will
become more similar to the CI prior. In particular, as $\lambda_{0}\to\infty$,
the HPP and the CI priors coincide.

Using the HPP (15), the joint posterior density may be written as

$\displaystyle p(\bm{\beta},\bm{m}|\bm{y},\lambda,\lambda_{0},\bm{\mu}_{0})$
$\displaystyle\propto\operatorname{\pi_{\text{CI}}}\left(\bm{\beta}\left|1+\lambda,\frac{\bm{y}+\lambda\bm{m}}{1+\lambda}\right)\right.\frac{Z(1+\lambda,\bm{y}+\lambda\bm{m})}{Z(\lambda,\lambda\bm{m})}\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0}),$
$\displaystyle\propto
p(\bm{\beta}|\bm{y},\lambda,\bm{m})p(\bm{m}|\bm{y},\lambda,\lambda_{0},\mu_{0}).$
(16)

Thus, the joint posterior density may be expressed as the product of the
posterior under the CI prior with shape parameter
$\frac{\bm{y}+\lambda\bm{m}}{1+\lambda}$ and a density over $\bm{m}$, which
depends on the observed data $\bm{y}$ and the prior prediction $\bm{\mu}_{0}$.
For intercept-only models, i.e., under the special case where
$\bm{\beta}=\beta_{0}$ and $\bm{X}=\bm{J}$, $\bm{m}=m$ is a scalar and it can
be shown that
$\operatorname{\pi_{\text{CI}}}(\bm{\beta}|1+\lambda,(y+\lambda\bm{m})/(1+\lambda))=\operatorname{\pi_{\text{DY}}}(\beta_{0}|n+\lambda,(n\bar{y}+\lambda
m)/(n+\lambda))$. Thus, the joint posterior (16) may be viewed as a
multivariate extension of the joint posterior for the i.i.d. setting, given by
(6).

The CI prior (10) is a member of a special class of priors called $g$-priors
(Zellner, 1986). The case where the hyperparameter $\lambda$ is treated as
random has been considered by Sabanés Bové and Held (2011) and further
expanded upon in (Held et al., 2015). When $\lambda$ is treated as random and
$\bm{m}$ is treated as fixed, the posterior distribution of $\lambda$ reflects
how compatible the prior prediction $\bm{m}$ is with the observed data. For
example, if many of the components of $\bm{m}$ differ from the observed data,
the posterior distribution of $\lambda$ may be concentrated near 0, so that
the prior has little effect in the posterior. By contrast, the HPP treats the
prior prediction $\bm{m}$ as random and $\lambda$ as fixed, and, thus, the
posterior distribution of the regression coefficients depends on the posterior
distribution of $\bm{m}$, which depends on both the prior prediction
$\bm{\mu}_{0}$ and the observed data. In Section 3.3, we discuss in more
detail the important role that the posterior distribution of $\bm{m}$ has in
the posterior distribution of the regression coefficients.

The HPP for GLMs is similar to a hierarchical conditional means prior (CMP)
for GLMs (Bedrick et al., 1996). For the CMP, $p$ _potential_ response and
covariate pairs $(\tilde{y}_{i},\tilde{\bm{x}}_{i})$ are elicited, where $p$
is the number of regression coefficients. The $\tilde{y}_{i}$’s may be
interpreted as a prior prediction for the mean response based on covariate
$\tilde{x}_{i}$. A prior inducing an a priori correlation structure on the
regression coefficients may be obtained by treating each
$\tilde{y}_{i},i=1,\ldots,p$ as random (for example, by utilizing the
hyperprior of the HPP). However, it may be difficult to justify a choice for
the $p$ potential covariate vectors $\tilde{x}_{i}$. By contrast,
specification of $\bm{\mu}_{0}$ is straightforward because it is simply a
prior prediction of the mean response for each of the $n$ observations,
potentially based on _observed_ (rather than potential) covariates.

###  3.2 The linear model

We now explore the posterior under the HPP for the normal linear model.
Specifically, we show that the marginal prior for the regression coefficients
under the HPP and the corresponding posterior density are multivariate normal.
Furthermore, we draw comparisons of the HPP with the CI prior and a
multivariate normal prior. We generalize Theorem 2.1 for the normal linear
model, showing that the limiting distribution of the posterior distribution of
$\bm{\beta}$ under the HPP as $\lambda\to\infty$ is the posterior density
under the CI prior.

Suppose we observe data $D=\\{(y_{i},\bm{x}_{i}),i=1,\ldots,n\\}$ with
likelihood function

$L(\bm{\beta},\tau|D)\propto\tau^{-n/2}\exp\left\\{-\frac{\tau}{2}(\bm{y}-\bm{X\beta})^{\prime}(\bm{y}-\bm{X\beta})\right\\},$

where $\tau=1/\sigma^{2}$ is a precision parameter, $\bm{\beta}$ is a
$p$-dimensional vector of regression coefficients,
$\bm{y}=(y_{1},\ldots,y_{n})^{\prime}$ is a $n$-dimensional response vector,
and $\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\prime}$ is a $n\times p$ matrix
(which may include an intercept term). We assume that $\tau$ is known and
fixed, and we may assume without loss of generality that $\tau=1$.

For the normal linear model, the CI prior is given by

$\displaystyle\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})\propto\exp\left\\{-\frac{\lambda}{2}\left[\bm{\beta}^{\prime}\bm{X}^{\prime}\bm{X}\bm{\beta}-2\bm{\beta}^{\prime}\bm{X}^{\prime}\bm{X}\hat{\bm{\beta}}_{\bm{m}}\right]\right\\},$
(17)

where
$\hat{\bm{\beta}}_{\bm{m}}=(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{m}$.
Note that, for a fixed value of $\bm{m}$, $\hat{\bm{\beta}}_{\bm{m}}$ is the
MLE for a normal linear model with vector of response variables $\bm{m}$ and
design matrix $\bm{X}$. The normalizing constant $Z(\lambda,\lambda\bm{m})$ of
the CI prior is given by

$\displaystyle Z(\lambda,\lambda\bm{m})$
$\displaystyle=(2\pi)^{p/2}\lvert\lambda^{-1}(\bm{X}^{\prime}\bm{X})^{-1}\rvert^{1/2}\exp\left\\{\frac{\lambda}{2}\left[\hat{\bm{\beta}}_{\bm{m}}^{\prime}\bm{X}^{\prime}\bm{X}\hat{\bm{\beta}}_{\bm{m}}\right]\right\\}.$
(18)

It follows from (17) and (18) that the the CI prior on $\bm{\beta}$ for a
fixed value of $\bm{m}$ is multivariate normal, i.e.,

$\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})=\phi_{p}(\bm{\beta}|\hat{\bm{\beta}}_{\bm{m}},\lambda^{-1}(\bm{X}^{\prime}\bm{X})^{-1}),$

where
$\phi_{p}(\bm{x}|\bm{a},\bm{C})=(2\pi)^{-p/2}\lvert\bm{C}\rvert^{-1/2}\exp\left\\{(\bm{x}-\bm{a})^{\prime}\bm{C}^{-1}(\bm{x}-\bm{a})\right\\}$
is the $p$-dimensional multivariate normal density with mean $\bm{a}$ and
positive definite covariance matrix $\bm{C}$. The hyperprior (13) is also a
multivariate normal density. That is,

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0})\propto\phi_{n}(\bm{m}|\bm{\mu}_{0},\lambda_{0}^{-1}\bm{I}_{n}),$
(19)

where $\bm{I}_{n}$ is a $n\times n$ identity matrix. Hence, a priori,
$\bm{m}=(m_{1},\ldots,m_{n})^{\prime}$ is a vector of $n$ independent normal
random variables whose $i^{th}$ component has mean $\mu_{0i}$ and variance
$\lambda_{0}^{-1}$.

We now derive the joint prior distribution of $\bm{\beta}$ and $\bm{m}$ for
the normal linear model. We show that the joint posterior distribution of
$(\bm{\beta},\bm{m})$ is multivariate normal. We further discuss how the
values of $\lambda_{0}$ and $\bm{\mu}_{0}$ affect the prior on $\bm{\beta}$.
Note that the marginal mean and covariance matrix of $\bm{\beta}$ based on the
HPP are given by

$E(\bm{\beta})=E[E(\bm{\beta}|\bm{m})]=E[(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{m}]=(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{\mu}_{0}$

and

$\operatorname{Cov}(\bm{\beta})=E\left[\operatorname{Cov}(\bm{\beta}|\bm{m})\right]+\operatorname{Cov}[E(\bm{\beta}|\bm{m})]=(\lambda^{-1}+\lambda_{0}^{-1})(\bm{X}^{\prime}\bm{X})^{-1},$

so that the joint prior of $(\bm{\beta},\bm{m})$ is multivariate normal, i.e.,

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{\beta},\bm{m})=\phi_{n+p}\left(\begin{pmatrix}\hat{\bm{\beta}}_{\bm{\mu}_{0}}\\\
\bm{\mu}_{0}\end{pmatrix},\begin{pmatrix}(\lambda^{-1}+\lambda_{0}^{-1})(\bm{X}^{\prime}\bm{X})^{-1}&\lambda_{0}^{-1}(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\\\
\lambda_{0}^{-1}\bm{X}^{\prime}(\bm{X}^{\prime}\bm{X})^{-1}&\lambda_{0}^{-1}\bm{I}_{n}\end{pmatrix}\right),$
(20)

where
$\bm{\hat{\beta}}_{\bm{\mu}_{0}}=(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{\mu}_{0}$
is the MLE for a normal linear model treating $\bm{\mu}_{0}$ as data. Hence,
for the normal linear model, the prior mean of $\bm{\beta}$ induced by the HPP
is precisely the same as the CI prior, treating $\bm{m}=\bm{\mu}_{0}$ as
fixed, and the prior variances of the components of $\bm{\beta}$ under the HPP
are larger than the CI prior. Moreover, (20) indicates that larger values of
$\lambda_{0}$ shrink the correlation between $\bm{\beta}$ and $\bm{m}$ toward
0.

Using (20), we may write the marginal prior of $\bm{\beta}$ as

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{\beta}|\lambda,\lambda_{0},\bm{\mu}_{0})=\int\operatorname{\pi_{\text{HPP}}}(\bm{\beta},\bm{m}|\lambda,\lambda_{0},\bm{\mu}_{0}))d\bm{m}=\phi_{p}(\bm{\beta}|\hat{\bm{\beta}}_{\bm{\mu}_{0}},(\lambda_{\text{H}}(\bm{X}^{\prime}\bm{X}))^{-1}),$
(21)

where $\lambda_{\text{H}}=(\lambda_{0}\lambda)/(\lambda_{0}+\lambda)$. Note
that the marginal prior (21) has the form of the CI prior with
$\bm{m}=\bm{\mu}_{0}$ and $\lambda=\lambda_{\text{H}}$. The connection between
the marginal prior of $\bm{\beta}$ based on the HPP and the CI prior holds
only when each component of $\bm{m}$ in the HPP has the same precision
parameter $\lambda_{0}$. While it is permissible to allow each component of
$\bm{m}$ to have its own degree of uncertainty in the HPP, the CI prior cannot
simultaneously reflect uncertainty about $\bm{m}$ and control the level of
discounting. Note also that the marginal posterior (21) may be obtained by
eliciting a multivariate normal prior on $\bm{\beta}$ with mean
$\hat{\bm{\beta}}_{\bm{\mu}_{0}}$ and covariance
$(\lambda_{\text{H}}\bm{X}^{\prime}\bm{X})^{-1}$ directly.

We now derive the posterior density of the regression coefficients for the
normal linear model. Since, as stated above, the marginal prior of
$\bm{\beta}$ is
$\operatorname{\pi_{\text{HPP}}}(\bm{\beta}|\lambda,\lambda_{0},\bm{\mu}_{0})\propto\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda_{\text{H}},\bm{\mu}_{0})$,
we have via the conjugacy of the CI prior that

$\displaystyle
p(\bm{\beta}|\bm{y},\lambda,\lambda_{0},\bm{\mu}_{0})\propto\operatorname{\pi_{\text{CI}}}\left(\bm{\beta}|1+\lambda_{\text{H}},\frac{\bm{y}+\lambda_{\text{H}}\bm{\mu}_{0}}{1+\lambda_{\text{H}}}\right)\propto\phi_{p}(\bm{\beta}|\bm{\mu}_{\bm{\beta}},(1+\lambda_{\text{H}})^{-1}(\bm{X}^{\prime}\bm{X})^{-1}),$
(22)

where

$\bm{\mu}_{\bm{\beta}}=\frac{1}{1+\lambda_{\text{H}}}\hat{\bm{\beta}}+\left(1-\frac{1}{1+\lambda_{\text{H}}}\right)\hat{\bm{\beta}}_{\bm{\mu}_{0}},$

where $\hat{\bm{\beta}}=(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}\bm{y}$ is
the MLE for $\bm{\beta}$. Hence, the marginal posterior distribution of
$\bm{\beta}$ is a multivariate normal distribution with mean
$\bm{\mu}_{\bm{\beta}}$ and covariance matrix
$(1+\lambda_{\text{H}})^{-1}(\bm{X}^{\prime}\bm{X})^{-1}$. Note that
$\bm{\mu}_{\bm{\beta}}$ is a convex combination of the MLE of $\bm{\beta}$
based on $\bm{y}$ and the “MLE” based on treating $\bm{\mu}_{0}$ as observed
data. Since $\lambda_{\text{H}}$ increases as either $\lambda_{0}$ or
$\lambda$ increases, larger values of $\lambda$ or $\lambda_{0}$ cause the
posterior mean of $\bm{\beta}$ to be closer to $\bm{\mu}_{0}$ than smaller
values. Note that when $\lambda\leq 1$,
$\lambda_{\text{H}}\leq\lambda_{0}/(1+\lambda_{0})<1$. Hence, the posterior
mean of $\bm{\beta}$ will always be closer to the MLE than
$\hat{\bm{\beta}}_{\bm{\mu}_{0}}$ when $\lambda\leq 1$.

When $\lambda_{\text{H}}$ is small, i.e., for small values of $\lambda$ or
$\lambda_{0}$, the posterior covariance matrix of $\bm{\beta}$ becomes similar
to that when using an improper uniform prior for the regression coefficients.
This is intuitive, as small values of $\lambda$ induce less borrowing from the
CI prior, and small values of $\lambda_{0}$ indicate substantial uncertainty
surrounding the prior prediction $\bm{\mu}_{0}$.

An interesting limiting case is when $\lambda\to\infty$. Note that
$\lim_{\lambda\to\infty}\lambda_{\text{H}}=\lambda_{0}$. Hence, it is clear
from (22) that

$\displaystyle
p(\bm{\beta}|\bm{y},\lambda,\lambda_{0},\bm{\mu}_{0})\to\phi_{p}\left(\bm{\beta}\left|\frac{\hat{\bm{\beta}}+\lambda_{0}\hat{\bm{\beta}}_{\bm{\mu}_{0}}}{1+\lambda_{0}},\frac{(\bm{X}^{\prime}\bm{X})^{-1}}{1+\lambda_{0}}\right)\right.\text{
as }\lambda\to\infty.$ (23)

Hence, as $\lambda\to\infty$, the posterior distribution of $\bm{\beta}$ is a
multivariate normal distribution with posterior mean equal to a weighted
average of $\hat{\bm{\beta}}$ and $\hat{\bm{\beta}}_{\bm{\mu}_{0}}$. The
relationship (23) indicates that, when $\lambda$ is allowed to unboundedly
increase, the posterior distribution is non-degenerate, but depends heavily on
the value of $\bm{\mu}_{0}$ when $\lambda_{0}$ is large. In practice, we
typically restrict $\lambda\in(0,1]$, so that the prior does not have more
weight on the posterior than the data. However, the form of the CI prior and
the relationship (23) indicates that the posterior of $\bm{\beta}$ will depend
more on the data than the prior prediction whenever $\lambda<1$ or
$\lambda_{0}<1$.

###  3.3 The posterior distribution of m

We now discuss the posterior distribution of $\bm{m}$. We will see that the
posterior distribution of $\bm{m}$ has important implications on the posterior
distribution of the regression coefficients. Namely, because the posterior
distribution of the regression coefficients under the HPP is the posterior
distribution of the CI prior averaged over the posterior of $\bm{m}$, both of
which depend on the observed data, the HPP tends to put more weight on the
observed data than the prior prediction in the posterior of the regression
coefficients compared to the CI prior.

We may write the posterior distribution of $\bm{m}$ under the HPP as

$\displaystyle
p(\bm{m}|\bm{y},\lambda,\lambda_{0},\bm{\mu}_{0})\propto\frac{Z(1+\lambda,\bm{y}+\lambda\bm{m})}{Z(\lambda,\lambda\bm{m})}\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0}),$
(24)

where
$Z(a,\bm{c})=\int\exp\left\\{a\left[(\bm{c}^{\prime}/a)\theta(\bm{X\beta})-b(\theta(\bm{X\beta}))\right]\right\\}$
is the normalizing constant of the CI prior (10), which has no general closed-
form solution. For the i.i.d. case, it was shown in Section 2.2 that, as
$\lambda\to\infty$, the posterior mean of $m$ was between the prior prediction
$\mu_{0}$ and the MLE $\hat{\mu}=\bar{y}$. Similarly, for GLMs, simulation
studies suggest that the posterior mean of $\bm{m}$ is between
$\hat{\bm{\mu}}=g^{-1}(\bm{X}\hat{\bm{\beta}})$ and $\bm{\mu}_{0}$, even for
finite values of $\lambda$.

Note that we may write the posterior distribution of $\bm{\beta}$ under the
HPP as

$\displaystyle
p(\bm{\beta}|\lambda,\lambda_{0},\bm{\mu}_{0})=\int\operatorname{\pi_{\text{CI}}}\left(\bm{\beta}|1+\lambda,\frac{\bm{y}+\lambda\bm{m}}{1+\lambda}\right)p(\bm{m}|\lambda,\lambda_{0},\bm{\mu}_{0})d\bm{m}.$
(25)

The relationship (25) indicates that the posterior distribution of
$\bm{\beta}$ may be interpreted as the CI density with shape parameter
$(\bm{y}+\lambda\bm{m})/(1+\lambda)$ averaged over the posterior distribution
of $\bm{m}$, which depends on the data $\bm{y}$. Thus, the posterior
distribution of $\bm{\beta}$ depends on the data in two ways: through the
conditional posterior of the CI prior (treating $\bm{m}$ as known) and through
the posterior of $\bm{m}$.

When $\lambda$ is treated as random and $\bm{m}$ is treated as fixed (e.g., as
done in Sabanés Bové and Held (2011)), the posterior distribution of $\lambda$
reflects how accurate the prior prediction is _overall_. For example, if many
of the prior predictions are inaccurate and some are highly accurate, the
posterior distribution of $\lambda$ may be concentrated near 0, so that the
prior has little influence on the posterior. Conversely, in the HPP, $\bm{m}$
is treated as random and $\lambda$ is fixed, so that the posterior always
borrows from the HPP. The HPP discounts at the prediction level, so that it is
possible for efficiency gains (e.g., lower mean square error) to be made even
when some of the prior predictions are inaccurate.

For the normal linear model, it can be shown that the posterior distribution
of $\bm{m}$ under the HPP is multivariate normal with mean
$E(\bm{m}|\bm{y})=\bm{\Lambda}\bm{\mu}_{0}+(\bm{I}_{n}-\bm{\Lambda})\hat{\bm{y}},$
where
$\bm{\Lambda}=\left(\lambda_{0}\bm{I}_{n}+\frac{\lambda}{1+\lambda}\bm{H}\right)^{-1}\lambda_{0}\bm{I}_{n}$,
$\bm{H}=\bm{X}(\bm{X}^{\prime}\bm{X})^{-1}\bm{X}^{\prime}$ is the orthogonal
projection operator onto the space spanned by the columns of $\bm{X}$, and
$\bm{\hat{y}}=\bm{X}\hat{\bm{\beta}}$ is the predicted values of $\bm{y}$,
where $\hat{\bm{\beta}}$ is the MLE of $\bm{\beta}$. Details of the derivation
are given in Section 3 of the Supplementary Materials. Hence, the posterior
mean of $\bm{m}$ for the normal linear model is a convex combination of the
prior prediction $\bm{\mu}_{0}$ and the predicted values $\hat{\bm{y}}$. Note
that for larger values of $\lambda$ (i.e., when more information from the
prior is borrowed) and for fixed values of $\lambda_{0}$, the posterior mean
of $\bm{m}$ depends more on the predicted values than the prior prediction.
This fact and the relationship (24) illustrate the robustness of the HPP,
namely, the posterior distribution of the regression coefficients is averaged
over a distribution depending on the prior prediction $\bm{\mu}_{0}$, but
highly modified by the data $\bm{y}$.

###  3.4 Computational development

We now discuss how to obtain posterior samples under the HPP. We develop a
Laplace approximation to the normalizing constant of the CI prior. We then
propose an efficient Hamiltonian Monte Carlo algorithm to obtain posterior
samples.

The joint posterior may be written as

$\displaystyle
p(\bm{\beta},\bm{m}|\bm{y},\lambda,\lambda_{0},\mu_{0})\propto\operatorname{\pi_{\text{CI}}}\left(\bm{\beta}\left|1+\lambda,\frac{\bm{y}+\lambda\bm{m}}{1+\lambda}\right)\right.\frac{\operatorname{\pi_{\text{HPP}}}(\bm{m}|\lambda_{0},\bm{\mu}_{0})}{Z(\lambda,\lambda\bm{m})}.$
(26)

While posterior inference in the i.i.d. case is analytically tractable, that
for the regression setting is more complicated because the normalizing
constant $Z(\lambda,\lambda\bm{m})$ in (26) does not have a closed form in
general.

We may utilize an integrated Laplace approximation to estimate the normalizing
constant. Sabanés Bové and Held (2011) took a similar approach to estimating
the marginal likelihood when $\lambda$ is treated as random. Note that
$\operatorname{\pi_{\text{CI}}}(\bm{\beta}|a,\bm{b}/a)$ is proportional to the
likelihood of a GLM with response variable $\bm{b}$ and inverse dispersion
parameter $a$ with link function $\theta$. Thus, we may utilize maximum
likelihood methods to efficiently obtain, for each proposed value
$\tilde{\bm{m}}$ of $\bm{m}$, $\hat{\bm{\beta}}_{\tilde{\bm{m}}}$, the value
of $\bm{\beta}$ that maximizes
$\operatorname{\pi_{\text{CI}}}(\bm{\beta}|1+\lambda,(\bm{y}+\lambda\bm{m})/(1+\lambda)$,
and $\mathcal{J}(\hat{\bm{\beta}}_{\tilde{\bm{m}}})$, the observed information
matrix evaluated at the maximizer. A Laplace approximation to the normalizing
constant $Z(\lambda,\lambda\bm{m})$ of the prior (10) is given by

$\displaystyle\hat{Z}_{\text{L}}(\lambda,\lambda\bm{m})\equiv(2\pi)^{1/2}\lvert\lambda\mathcal{J}(\hat{\bm{\beta}}_{\bm{m}})\rvert^{-1/2}\exp\left\\{\lambda\left[\bm{m}^{\prime}\theta(\bm{X}\hat{\bm{\beta}}_{\bm{m}})-\bm{J}^{\prime}b(\theta(\bm{X}\hat{\bm{\beta}}_{\bm{m}}))\right]\right\\}.$

Because $\hat{Z}_{\text{L}}$ is a nonstochastic estimator of the normalizing
constant, we may utilize Hamiltonian MCMC for posterior sampling, such as the
highly efficient No U-Turn Sampler (NUTS) of Hoffman and Gelman (2014).

Under the special case that all covariates are categorical, there are closed-
form expressions for the normalizing constants, and sampling from the
posterior (26) via MCMC is straightforward. Suppose without loss of generality
that there is a single categorical covariate with $J$ levels. Let
$\bm{y}_{j}=(y_{1j},\ldots,y_{n_{j},j})$ denote the $n_{j}$ response variables
belonging to the $j^{th}$ category, $j=1,\ldots,J$. The likelihood may be
written as

$\displaystyle
L(\bm{\beta}|\bm{y})\propto\prod_{j=1}^{J}\exp\left\\{n_{j}\left[\beta_{j}\bar{y}_{j}-b(\beta_{j})\right]\right\\},$
(27)

where $\bar{y}_{j}=n_{j}^{-1}\sum_{i=1}^{n_{j}}y_{ij}$ is the sample mean for
the $j^{th}$ category. Note that the likelihood (27) is proportional to the
product of $J$ exponential family likelihoods (1), where observations are
i.i.d. within each category. It follows that the joint conjugate prior of
$\bm{\beta}$ is the product of $J$ independent DY priors, i.e.,

$\displaystyle\pi(\bm{\beta}|\lambda,\bm{m})\propto\prod_{j=1}^{J}\frac{1}{Z(\lambda,\lambda
m_{j})}\exp\left\\{\lambda\left[m_{j}\beta_{j}-b(\beta_{j})\right]\right\\},$
(28)

where $\bm{m}=(m_{1},\ldots,m_{J})^{\prime}$ is a $J$-dimensional vector.
Since (28) is a product of independent conjugate priors with a scalar shape
parameter, the normalizing constant for each class has a known closed form
(e.g., it is a beta function for a logistic regression model).

##  4 Utilizing previous studies

In this section, we describe how the HPP may be utilized when we possess
summary statistics from a previous study obtained, for example, from a
publication. We compare and contrast our approach with the power prior of
Ibrahim and Chen (2000).

Suppose that we possess historical data $D_{0}=\\{y_{0i},i=1,\ldots,n_{0}\\}$
with density (or mass function) (1). For notational ease, we assume that
dispersion parameters are known and fixed. The power prior of Ibrahim and Chen
(2000) in i.i.d. settings is given by

$\displaystyle\pi_{\text{PP}}(\theta|a_{0},D_{0})\propto\exp\left\\{a_{0}\left[\theta\bar{y}_{0}-b(\theta)\right]\right\\}\pi_{0}(\theta),\theta\in\Theta,$
(29)

where $a_{0}\in(0,n]$ is a parameter controlling how much the practitioner
wishes to borrow from the historical data and $\pi_{0}(\theta)$ is an initial
prior for $\theta$. If $\pi_{0}(\theta)\propto 1$, then the power prior (29)
is precisely the DY prior (2) with $\lambda=a_{0}$ and $m=\bar{y}_{0}$. Thus,
for i.i.d. exponential family models, one may obtain the same posterior
density using the summary statistic $S_{0}=\bar{y}$ as that using the entire
historical data set $D_{0}$.

Suppose now that we possess historical data with covariates, say,
$D_{0}=\\{(y_{0i},\bm{x}_{0i}),i=1,\ldots,n_{0}\\}$. The power prior for GLMs
is given by

$\displaystyle\pi_{\text{PP}}(\bm{\beta}|a_{0},D_{0})\propto\exp\left\\{a_{0}\left[\bm{y}_{0}^{\prime}\theta(\bm{X}_{0}\bm{\beta})-\bm{J}^{\prime}b(\theta(\bm{X}_{0}\bm{\beta}))\right]\right\\}\pi_{0}(\bm{\beta}),\bm{\beta}\in\mathbb{R}^{p},$
(30)

where $a_{0}\in(0,1]$ is described above and $\pi_{0}$ here is an initial
prior for $\bm{\beta}$, which we may take as $\pi_{0}(\bm{\beta})\propto 1$.
The prior (30) is similar to CI prior with $\lambda=a_{0}$, but $\bm{y}_{0}$
is a $n_{0}$-dimensional vector while $\bm{m}$ in the CI prior is a
$n$-dimensional vector and the covariates from the historical data set,
$\bm{X}_{0}$, are utilized instead of those from the current data, $\bm{X}$.
Note that for GLMs, unlike the i.i.d. case, the power prior requires the full
historical data set.

We may, however, utilize the MLE from the historical data,
$\hat{\bm{\beta}}_{0}$, to obtain a prior prediction for the CI prior as
$\bm{\mu}_{0}=g^{-1}(\bm{X}\hat{\bm{\beta}}_{0})$. However,
$\hat{\bm{\beta}}_{0}$ is a statistic and each of its components has a
variance, and there is no direct mechanism for implementing this uncertainty
into the CI prior.

By contrast, using the HPP, we may elicit uncertainty surrounding the prior
prediction $\bm{\mu}_{0}$ via the precision parameter $\lambda_{0}$. Using the
delta method, we may approximate the variance of
$\mu_{0i}=g^{-1}(\bm{x}_{i}^{\prime}\hat{\bm{\beta}}_{0})$ as

$\displaystyle\text{Var}(g^{-1}(\bm{x}_{i}^{\prime}\hat{\bm{\beta}}_{0}))\approx\frac{\bm{x}_{i}^{\prime}\text{Cov}(\hat{\bm{\beta}}_{0})\bm{x}_{i}}{\left[\dot{g}\circ
g^{-1}(\bm{x}_{i}^{\prime}\hat{\bm{\beta}}_{0})\right]^{2}}\equiv\tau_{0i},i=1,\ldots,n.$
(31)

Typically, previous studies report only the estimated standard errors,
$\\{\hat{\sigma}_{0j},j=1,\ldots,p\\}$, of $\hat{\bm{\beta}}_{0}$ instead of
its estimated covariance matrix. In such cases, we may substitute
$\widehat{\bm{\Sigma}}=\text{diag}\\{\hat{\sigma}_{0j}^{2},j=1,\ldots,p\\}$
for $\text{Cov}(\hat{\bm{\beta}}_{0})$ in (31), providing a reasonable
approximation for $\tau_{0i}$ in (31), say, $\hat{\tau}_{0i}$ as

$\displaystyle\hat{\tau}_{0i}=\frac{\sum_{j=1}^{p}x_{ij}^{2}\hat{\sigma}_{0j}^{2}}{\left[\dot{g}\circ
g^{-1}(\bm{x}_{i}^{\prime}\hat{\bm{\beta}}_{0})\right]^{2}},i=1,\ldots,n.$

Once an estimate $\hat{\tau}_{0i}$ for $\tau_{0i}$ is obtained, we may elicit
$\mu_{0i}$ and $\lambda_{0i}$, such that $\text{Var}(m_{i})=\hat{\tau}_{0i}$,
$i=1,\ldots,n$. Here, we allow each component of $\bm{m}$ to have its own
precision. That is, we augment the hyperprior (13) to

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{m}|\bm{\lambda}_{0},\bm{\mu}_{0})\propto\prod_{i=1}^{n}\exp\left\\{\lambda_{0i}\left[\dot{b}^{-1}(m_{i})\mu_{0i}-(b\circ\dot{b}^{-1})(m_{i})\right]\right\\}\frac{1}{v(m_{i})},$

where $\bm{\lambda}_{0}=(\lambda_{01},\ldots,\lambda_{0n})^{\prime}$ is a
$n$-dimensional vector of precision parameters. We stress that the
hyperparameters $\bm{\lambda}_{0}$ and $\bm{\mu}_{0}$ in the example described
are not elicited based on an individual’s opinion. Rather, the hyperparameters
are deterministic functions of summary statistics from a previous study.

Both the HPP and the power prior of Ibrahim and Chen (2000) induce an a priori
correlation structure on the regression coefficients. However, with the power
prior, the correlation structure depends on the observed historical response
variable $\bm{y}_{0}$ and design matrix $\bm{X}_{0}$. By contrast, the
correlation structure in the HPP is based on the current design matrix
$\bm{X}$, the prior prediction
$\bm{\mu}_{0}=g^{-1}(\bm{X}\hat{\bm{\beta}}_{0})$, and the prior precision
$\bm{\lambda}_{0}$. When one is in possession of an entire historical data
set, the power prior may be used, which is less computationally demanding and
has various desirable properties (see, e.g., Ibrahim et al. (2015)). However,
the HPP may be used more generally, as it only requires a prior prediction for
the mean response and associated uncertainty, which may be obtained either
from a full historical data set or using summary statistics from a
publication.

When the parameter $a_{0}$ in the power prior (30) is treated as fixed, the
prior can be highly influential on the posterior density. For example, in
Section 6, we generate a historical data set with a positive treatment effect
and a current data set with a null treatment effect. For these incompatible
data sets, the power prior suggested that treatment was all but certain to be
efficacious a posteriori. By contrast, the posterior density of the treatment
effect under the HPP suggested that treatment may be unbeneficial, which is
the correct result based on how the current data were generated.

##  5 Data application

In this section, we obtain posterior quantities under the HPP using the data
of Finney (1947) (henceforth referred to as the Finney data), which was
obtained from the R package _robustbase_ (Maechler et al., 2021). We compare
posterior quantities with the CI prior. In particular, we illustrate how the
HPP is less sensitive to the prior prediction than the CI prior.

The Finney data were obtained from a controlled study in which the effects of
the rate and volume of air from a single deep breath on whether a transient
reflex vaso-constriction on the skin of the digits occurred. The response
variable, $y$, is binary, such that $y=1$ if a vaso-constriction occurred and
$y=0$ if one did not occur. The number of observations recorded is $n=39$.
There are two continuous covariates: $x_{1}=\log(\text{volume})$ and
$x_{2}=\log(\text{rate})$, where “volume” is the inhaled volume of the
individual and “rate” is the rate of inhalation. We assume a logistic
regression model for the responses. The likelihood may thus be written as

$\displaystyle
L(\bm{\beta}|\bm{y},\bm{X})=\prod_{i=1}^{n}\exp\left\\{y_{i}\bm{x}_{i}^{\prime}\bm{\beta}-\log(1+\exp\\{\bm{x}_{i}^{\prime}\bm{\beta}\\})\right\\},$

where $\bm{x}_{i}=(1,x_{1i},x_{2i})^{\prime}$, $i=1,\ldots,n$,
$\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{n})^{\prime}$ is the $n\times 3$ design
matrix, and $\bm{\beta}=(\beta_{0},\beta_{1},\beta_{2})$, where $\beta_{0}$ is
an intercept term. The CI prior takes the form

$\operatorname{\pi_{\text{CI}}}(\bm{\beta}|\lambda,\bm{m})\propto\prod_{i=1}^{n}\exp\left\\{\lambda\left[m_{i}\bm{x}_{i}^{\prime}\bm{\beta}-\log(1+\exp\\{\bm{x}_{i}^{\prime}\bm{\beta}\\})\right]\right\\},$

where $\bm{m}=(m_{1},\ldots,m_{n})^{\prime}\in(0,1)^{n}$ is the shape
parameter and $\lambda\in(0,1]$ is the precision parameter. The hierarchical
model is completed by specifying the hyperprior of the HPP, which is given by

$\displaystyle\operatorname{\pi_{\text{HPP}}}(\bm{m}|\bm{\mu_{0}},\lambda_{0})\propto\prod_{i=1}^{n}m_{i}^{\lambda_{0}\mu_{0i}}(1-m_{i})^{\lambda_{0}(1-\mu_{0i})},$
(32)

where $\bm{\mu}\in(0,1)^{n}$ is the prior prediction for the mean of $\bm{y}$
and $\lambda_{0}>0$ is the associated precision parameter. Note that (32) is
proportional to a collection of $n$ independent beta random variables each
with mean $\mu_{0i}$ and shared precision parameter $\lambda_{0}$, where
$\mu_{0i}$ denotes the $i^{th}$ component of $\bm{\mu}_{0}$.

For the purposes of illustration, we suppose that an expert informs us that
lower volume, on average, yields a lower probability for the event. However,
the expert is unsure about the effect of the rate of inhalation, and declines
to make a prediction for the probability of an event based on that covariate.
The expert then gives their prediction based on the level of volume, providing
a value of $\mu_{0i}=\mu_{0}$ when volume is less than 1 and
$\mu_{0i}=\mu_{1}$ when volume is larger than 1, for $i=1,\ldots,n$. The
expert further dictates their uncertainty surrounding these values, which is
captured by $\lambda_{0}>0$.

The following hyperparameters for the exercise were used:
$(\mu_{0},\mu_{1})\in\\{(0.3,0.7),(0.4,0.6)\\}$, $\lambda_{0}\in\\{4,10\\}$,
and $\lambda\in\\{0.75,1.00\\}$. We note that the choices for prior mean were
arbitrary and they need not be symmetric in general. Models were fit using the
_rstan_ package in R (Stan Development Team, 2020), which utilizes the No
U-Turn Sampler (NUTS) of Hoffman and Gelman (2014). The Laplace method for
estimating the normalizing constant in Section 3.4 was utilized. A total of
25,000 samples were obtained, for which 1,000 was used as a burn-in period. No
thinning was performed. The minimum effective sample size (ESS) for all
samples of $\bm{\beta}$ across all scenarios was 9,269.

Table 1, shows the results from the analysis of the Finney data, comparing the
posterior under the HPP against the CI prior. The table shows that posterior
quantities utilizing the HPP are more similar to those using the CI prior for
increasing values of $\lambda_{0}$. For example, for
$(\mu_{0},\mu_{1})=(0.3,0.7)$ and when $\lambda=0.75$, the posterior means for
$\beta_{1}$ are, respectively, $2.72$, $2.63$, and $2.56$ for $\lambda_{0}$
equal to $4$, $10$, and $\infty$. For example, when $\lambda_{0}=4$, the
posterior mean for $\beta_{1}$ is $2.72$ and $2.56$, respectively, when
$\lambda=0.75$, and it is $2.53$ and $2.35$, respectively, when $\lambda=1$. A
heuristic reason for this is because as $\lambda$ increases, the posterior
depends more on the prior. Smaller values of $\lambda_{0}$ suggest more
uncertainty in the prior prediction $\mu_{0}$. Hence, for fixed values of
$\lambda$, lower values of $\lambda_{0}$ will make the posterior distribution
under the HPP more dissimilar to that under the CI prior.

The primary advantage of the HPP over the CI prior is that it enables
practitioners to elicit an appropriate amount of uncertainty surrounding the
prior prediction. As previously stated, in general, each component of the
prior prediction may have its own variance. Moreover, Table 1 clearly shows
that the results become similar when the precision parameter $\lambda_{0}$
increases. However, due to the estimation of the normalizing constant of the
conjugate prior in the MCMC algorithm, the HPP is more computationally
expensive. Despite this, we were able to obtain a large amount of samples in a
matter of a few minutes.

Table 1: Posterior quantities from Finney data

|  |  | $(\mu_{0},\mu_{1})$
---|---|---|---
|  |  | $(0.3,0.7)$ | $(0.4,0.6)$
$\lambda_{0}$ | $\lambda$ | Parameter | Mean | SD | HPD | Mean | SD | HPD
4 | 0.75 | $\beta_{0}$ | $-0.858$ | $\phantom{-}0.448$ | $(-1.837,-0.074)$ | $-0.796$ | $\phantom{-}0.437$ | $(-1.741,-0.019)$
|  | $\beta_{1}$ | $\phantom{-}2.721$ | $\phantom{-}0.798$ | $(~{}~{}1.308,~{}~{}4.439)$ | $\phantom{-}2.230$ | $\phantom{-}0.714$ | $(~{}~{}0.951,~{}~{}3.740)$
|  | $\beta_{2}$ | $\phantom{-}1.403$ | $\phantom{-}0.584$ | $(~{}~{}0.431,~{}~{}2.711)$ | $\phantom{-}1.353$ | $\phantom{-}0.575$ | $(~{}~{}0.378,~{}~{}2.632)$
|  1.00 | $\beta_{0}$ | $-0.736$ | $\phantom{-}0.390$ | $(-1.560,-0.031)$ | $-0.672$ | $\phantom{-}0.383$ | $(-1.489,~{}~{}0.014)$
|  | $\beta_{1}$ | $\phantom{-}2.529$ | $\phantom{-}0.725$ | $(~{}~{}1.211,~{}~{}4.063)$ | $\phantom{-}1.998$ | $\phantom{-}0.650$ | $(~{}~{}0.803,~{}~{}3.359)$
|  | $\beta_{2}$ | $\phantom{-}1.188$ | $\phantom{-}0.499$ | $(~{}~{}0.330,~{}~{}2.285)$ | $\phantom{-}1.134$ | $\phantom{-}0.500$ | $(~{}~{}0.282,~{}~{}2.250)$
10 | 0.75 | $\beta_{0}$ | $-0.793$ | $\phantom{-}0.412$ | $(-1.668,-0.051)$ | $-0.734$ | $\phantom{-}0.408$ | $(-1.603,~{}~{}0.003)$
|  | $\beta_{1}$ | $\phantom{-}2.631$ | $\phantom{-}0.760$ | $(~{}~{}1.262,~{}~{}4.244)$ | $\phantom{-}2.131$ | $\phantom{-}0.692$ | $(~{}~{}0.874,~{}~{}3.587)$
|  | $\beta_{2}$ | $\phantom{-}1.291$ | $\phantom{-}0.525$ | $(~{}~{}0.398,~{}~{}2.440)$ | $\phantom{-}1.242$ | $\phantom{-}0.528$ | $(~{}~{}0.324,~{}~{}2.395)$
|  1.00 | $\beta_{0}$ | $-0.669$ | $\phantom{-}0.368$ | $(-1.455,~{}~{}0.006)$ | $-0.613$ | $\phantom{-}0.360$ | $(-1.370,~{}~{}0.046)$
|  | $\beta_{1}$ | $\phantom{-}2.422$ | $\phantom{-}0.682$ | $(~{}~{}1.170,~{}~{}3.826)$ | $\phantom{-}1.900$ | $\phantom{-}0.622$ | $(~{}~{}0.761,~{}~{}3.193)$
|  | $\beta_{2}$ | $\phantom{-}1.075$ | $\phantom{-}0.453$ | $(~{}~{}0.288,~{}~{}2.068)$ | $\phantom{-}1.032$ | $\phantom{-}0.453$ | $(~{}~{}0.240,~{}~{}2.013)$
$\infty$ | 0.75 | $\beta_{0}$ | $-0.747$ | $\phantom{-}0.396$ | $(-1.592,-0.033)$ | $-0.682$ | $\phantom{-}0.385$ | $(-1.484,~{}~{}0.015)$
|  | $\beta_{1}$ | $\phantom{-}2.563$ | $\phantom{-}0.729$ | $(~{}~{}1.228,~{}~{}4.080)$ | $\phantom{-}2.039$ | $\phantom{-}0.658$ | $(~{}~{}0.837,~{}~{}3.410)$
|  | $\beta_{2}$ | $\phantom{-}1.212$ | $\phantom{-}0.498$ | $(~{}~{}0.352,~{}~{}2.312)$ | $\phantom{-}1.159$ | $\phantom{-}0.492$ | $(~{}~{}0.305,~{}~{}2.223)$
|  1.00 | $\beta_{0}$ | $-0.626$ | $\phantom{-}0.341$ | $(-1.331,-0.003)$ | $-0.560$ | $\phantom{-}0.331$ | $(-1.250,~{}~{}0.050)$
|  | $\beta_{1}$ | $\phantom{-}2.351$ | $\phantom{-}0.653$ | $(~{}~{}1.154,~{}~{}3.707)$ | $\phantom{-}1.801$ | $\phantom{-}0.579$ | $(~{}~{}0.731,~{}~{}2.985)$
|  | $\beta_{2}$ | $\phantom{-}0.998$ | $\phantom{-}0.415$ | $(~{}~{}0.256,~{}~{}1.884)$ | $\phantom{-}0.939$ | $\phantom{-}0.415$ | $(~{}~{}0.204,~{}~{}1.832)$

The HPP is less sensitive to the prior prediction than the CI prior. For
example, suppose now the expert believes that higher levels of volume _reduce_
the probability of an event. This expert elicits
$(\mu_{0},\mu_{1})\in\\{(0.6,0.4),(0.7,0.3)\\}$ with the same levels of
precision as before. The $95\%$ highest posterior density (HPD) regions under
the prior predictions are presented in Figure 1 of the Supplementary
Materials. Under the prior response elicitation least consistent with the
data, e.g., for $(\mu_{0},\mu_{1})=(0.7,0.3)$, the left endpoint of the HPD
interval for $\beta_{1}$ is only positive when $\lambda=0.75$ and
$\lambda_{0}=4$, i.e., when there is more uncertainty surrounding the prior
prediction and when the degree of informativeness of the prior is moderate.

##  6 Data application with a previous study

In this section, we show results applying the proposed HPP against other
priors for two generated historical data sets. The other selected priors are
the CI prior, the power prior (PP) of Ibrahim and Chen (2000), and a “Gaussian
power prior” (GPP), which is an approximation to the asymptotic power prior of
Ibrahim et al. (2015) and is given by

$\pi_{\text{GPP}}(\bm{\beta}|\bm{\mu}_{\bm{\beta}},\bm{\Sigma}_{\bm{\beta}},\lambda)\propto\left[(2\pi)^{-1/2}\lvert\bm{\Sigma}_{\bm{\beta}}\rvert^{-1/2}\exp\left\\{-\frac{1}{2}(\bm{\beta}-\bm{\mu}_{\bm{\beta}})^{\prime}\bm{\Sigma}_{\bm{\beta}}^{-1}(\bm{\beta}-\bm{\mu}_{\bm{\beta}}\right\\}\right]^{\lambda},$

where $\bm{\mu}_{\bm{\beta}}$ is the prior mean of $\bm{\beta}$ and
$\bm{\Sigma}_{\bm{\beta}}$ is the prior covariance matrix of $\bm{\beta}$. We
assume for the GPP that only the maximum likelihood estimates and associated
standard errors are available from the previous study, and we elicit
$\bm{\mu}_{\bm{\beta}}=\hat{\bm{\beta}}$ and
$\bm{\Sigma}_{\bm{\beta}}=\text{diag}\\{\hat{\sigma}_{j},j=1,\ldots,p\\}$,
where $\hat{\bm{\beta}}=(\hat{\beta}_{1},\ldots,\hat{\beta}_{p})^{\prime}$ is
a vector of MLEs for the regression coefficients and $\hat{\sigma}_{j}$ is the
standard error for the $j^{th}$ coefficient, $j=1,\ldots,p$. The GPP is
precisely equivalent to the asymptotic power prior, except the off-diagonal
elements of $\bm{\Sigma}_{\bm{\beta}}$ are set to 0. The reason for this is,
in practical situations where a prior is being formulated on the basis of
results in published studies, the correlations of the regression coefficients
are not typically reported.

The PP will use a full historical data set, while hyperparameter elicitation
for the three remaining priors will utilize the maximum likelihood estimates
(MLEs) and standard errors. For the HPP, $\bm{\lambda}_{0}$ will be elicited
based on the standard errors as described in Section 4. For the GPP, the prior
variance of the regression coefficients will be set equal to the squared
standard errors.

###  6.1 The generated data sets

The ACTG036 study (Merigan et al., 1991) was a clinical trial comparing AZT
with a placebo in asymptomatic patients with hereditary coagulation disorders
(hemophilia). The outcome variable was CD4 count, defined to be the number of
CD4 cells (a type of white blood cell that destroys bacteria and fights off
infection) per cubic millimeter of blood. Covariates included in the model
were treatment ($x_{1i}=1$ if subject $i$ received AZT, 0 otherwise), race
$(x_{2i}=1$ if subject $i$ was white, 0 otherwise), and age ($x_{3i})$,
treated as a continuous variable. We assume a Poisson regression model with
canonical link. The likelihood function for the current data is given by

$\displaystyle
L(\bm{\beta}|\bm{y},\bm{X})=\exp\left\\{\sum_{i=1}^{n}y_{i}\bm{x}_{i}^{\prime}\bm{\beta}-\exp\\{\bm{x}_{i}^{\prime}\bm{\beta}\\}\right\\},$

where $n=75$ is the sample size of the current data, $y_{i}$ is the CD4 count
for subject $i$, $\bm{x}_{i}=(1,x_{1i},x_{2i},x_{3i})^{\prime}$, and
$\bm{\beta}=(\beta_{0},\beta_{1},\beta_{2},\beta_{3})^{\prime}$. Henceforth,
we denote the current data by $D=\\{(y_{i},\bm{x}_{i}),i=1,\ldots,n\\}$.
Historical data of size $n_{0}=50$ was generated from the same model, giving
the data set $D_{0}=\\{(y_{0i},\bm{x}_{0i}),i=1,\ldots,n_{0}\\}$. The
quantities $\beta_{0},\beta_{2},$ and $\beta_{3}$ were obtained directly from
the MLEs of the ACTG036 data set. For the historical data, the treatment
effect was set equal to the MLE of the ACTG036 study, that is,
$\beta_{1}\approx 0.048$. Two current data sets were generated: an
“incompatible” current data set, where $\beta_{1}=0$, and a “compatible” data
set, where $\beta_{1}\approx 0.048$.

###  6.2 Results

Prior and posterior density plots for the treatment effect for the compatible
and incompatible current data sets are depicted in Figure 1. The left panel
depicts prior and posterior densities when the historical and current data are
incompatible. The right panel shows the prior and posterior densities when the
data sets are compatible. Of the four prior densities, the HPP yields the
highest prior variance for the treatment effect $\beta_{1}$ (i.e., it is the
least informative prior). When $\lambda\in\\{0.75,1.00\\}$, the HPP is the
only prior that places mass around the null, while the other three priors
suggest that treatment is all but certain to be efficacious a priori. Since
the CI prior, which is the most informative prior, is a special case of the
HPP, we see that, for a fixed value of $\lambda$, the HPP was the most
flexible prior for expressing uncertainty.

Figure 1: Prior and posterior densities of the treatment effect for the
generated data under various levels of borrowing $(\lambda)$. The historical
data suggest efficacy ($\beta_{1}=0.048)$. The incompatible current data was
generated from $\beta_{1}=0$, while that for the compatible was generated from
$\beta_{1}=0.048$.

When the current data is incompatible with the historical data, i.e., when
$\beta_{1}=0$ for the current data and $\beta_{1}\approx 0.048$ for the
historical data, the posterior distribution of $\beta_{1}$ utilizing the HPP
is generally more robust than its competitors (in the sense of placing the
most mass a posteriori on a nonpositive treatment effect). When enabling the
prior to have as much influence on the posterior as the data, (i.e., when
$\lambda=1.00$), the posterior densities under the CI prior, GPP, and PP
suggest treatment is all but certain to be efficacious a posteriori. In
contrast, the posterior density under the HPP suggests that treatment could be
unbeneficial, which is the correct result.

When the two data sets are compatible, i.e., when $\beta_{1}\approx 0.048$ for
both the current and historical data sets, the four posterior densities are
quite similar. Across the various levels of $\lambda$, the posterior densities
from all four priors suggest that treatment is efficacious a posteriori. The
posterior means across the four priors and all levels of $\lambda$ are
essentially the same, although the GPP yielded a somewhat higher posterior
mean than the other priors. While increasing $\lambda$ reduced the posterior
variance markedly for the three competitor priors, the posterior under the HPP
was virtually unchanged by the level of $\lambda$.

In general, for this application, the effect of $\lambda$ on the prior and the
posterior densities of the treatment effect under the HPP is relatively
marginal. When the data sets are incompatible, increasing the value of
$\lambda$ from $0.50$ to $0.75$ increases the posterior mean of the treatment
effect by $14\%$. By contrast, the same change for the power prior yields
roughly a $29\%$ increase, more than twice as much. This suggests that, when
each component of the prior prediction for the HPP is given its own level of
precision, the posterior density under the HPP is much less sensitive to
increasing values of $\lambda$ than the other priors.

The primary advantage of the HPP over the CI prior is that it allows
practitioners to incorporate uncertainty in the prior guess in a flexible way.
However, the HPP adds computational cost since the normalizing constant of the
CI prior must be estimated. The GPP induces an a priori independent prior on
the regression coefficients, which is not realistic practically. By contrast,
the HPP induces a correlation structure a priori on the regression
coefficients. The primary advantage of the PP over the HPP is that the PP is
more computationally efficient than the GPP. However, we have seen that the PP
can be quite informative when $\lambda$ is fixed, with high posterior
probability that $\beta_{1}$ is positive. By contrast, the posterior under the
HPP suggests a nonzero chance that $\beta_{1}\leq 0$ when the data sets are
compatible a posteriori. Moreover, the power prior is somewhat restrictive in
that it requires access to a full historical data set.

###  6.3 The posterior of m

In Section 3.3, we argued that the posterior mean of $\bm{m}$ is “between” the
prior prediction $\bm{\mu}_{0}$ and the data. We showed that the posterior
mean of $\bm{m}$ for the normal linear model was a convex combination of
$\hat{\bm{y}}$, the predicted values, and $\bm{\mu}_{0}$. While a formal proof
may not be offered, based on a data set analysis, it appears that, for other
GLMs, the posterior mean of $\bm{m}$ is also between $\bm{\mu}_{0}$ and
$\hat{\bm{\mu}}=g^{-1}(\bm{X}\hat{\bm{\beta}})$, the predicted mean based on
the MLE.

Figure 2 depicts $\bm{\mu}_{0}=g^{-1}(\bm{X}\hat{\bm{\beta}}_{0})$, the
prediction for $E(\bm{y})$ based on the MLE of the historical data (where
$\beta_{1}=0.048)$, and $\hat{\bm{\mu}}$, the predicted mean utilizing the MLE
of the incompatible current data (where $\beta_{1}=0$), against the posterior
mean of $\bm{m}$. The black line is a 45-degree line, i.e., points close to
the black line are close to the posterior mean of $\bm{m}$.

The left panel of Figure 2 depicts the predicted means for the historical data
set and the posterior of $\bm{m}$ for observations in the control group. For
this group, the historical data set is compatible since the covariate
corresponding to $\beta_{1}$ (i.e., the treatment indicator) is zero. Note
that the vertical distance between the blue and orange points are relatively
small, and both points are close to the 45-degree line. This indicates that
the posterior mean of $\bm{m}$ for the control group is very close to both
$\hat{\bm{\mu}}$ and $\bm{\mu}_{0}$, i.e., that $\bm{\mu}_{0}$ was an accurate
guess for $E(\bm{y}_{0})$, where $\bm{y}_{0}$ is the vector of response
variables for observations assigned to control.

The right panel of Figure 2 depicts the posterior mean of $\bm{m}$ against the
predicted mean based on the current data (in blue) and historical data (in
orange). In this case, the orange points are always above the blue points.
This is because, in the current data, $\beta_{1}=0$, whereas $\beta_{1}=0.048$
in the historical data set (i.e., the predicted mean under the historical data
is larger than the current data).

Figure 2: Posterior mean of $\bm{m}$ for the Poisson example in Section 6 with
$\lambda=1$. The blue points are the predicted mean utilizing the current data
set, which suggest $\beta_{1}=0$. The orange points are the predicted mean
utilizing the historical data set, which suggest $\beta_{1}=0.048$. The black
line is a 45-degree line.

##  7 Discussion

We have developed Bayesian hierarchical generalized linear models utilizing
conjugate priors, where the mean parameter is treated as random. The HPP is
simple and intuitive (e.g., for logistic regression models, the hyperprior is
simply a product of independent beta priors, which is, ignoring covariates, a
conjugate prior for each component of the response variable). However, we note
that the development is flexible–any (proper) hyperprior with support on the
mean of the responses may be utilized in practice. For example, a normal
distribution truncated to $[0,1]$ may be utilized for logistic regression
models.

The incorporation of uncertainty in the prior prediction of the mean of the
responses has a natural practical application. As the prior prediction for the
mean of the response typically comes in the form of historical data and/or
expert opinion, there is uncertainty surrounding the elicited value. Our
proposed model enables practitioners to capture this uncertainty. For example,
in rare disease settings, the prior prediction may be elicited by a panel of,
say, $J$ experts, where $J\geq 1$. The opinions of these experts may be
aggregated to a single prior (mathematically or behaviorally), to elicit a
predicted response probability $\bm{\mu}_{0}$ and incorporate uncertainty
surrounding this value via $\lambda_{0}$. We have also explored the
relationship of the posterior densities under the HPP and the conjugate prior.

We have shown several advantages of the HPP over other priors. First, unlike
the CI prior, the HPP allows practitioners to fix the degree of influence of
the prior while simultaneously allowing for uncertainty surrounding the prior
prediction. Second, the HPP may be utilized to set a correlated joint prior on
regression coefficients, for example, in settings where summary statistics
from a previous study are available, but the full historical data set is not.
By contrast, the power prior (Ibrahim and Chen, 2000) requires the possession
of a full historical data set. Finally, the HPP may be less informative than
the CI prior, power prior, and GPP, which may be attractive to clinical
regulators. In particular, data applications indicate that posterior
quantities using the HPP are robust against prior predictions that are not
compatible with the observed data.

An extension of this work could be treating both $\bm{m}$ and $\lambda$ in the
CI prior as random. When $\bm{m}$ is fixed and $\lambda$ is random, where
$\lambda$ controls for the level of borrowing from the prior. If the prior
prediction is incompatible with the data set, the posterior of $\lambda$ would
be near zero. However, when $\bm{m}$ and $\lambda$ are both treated as random,
it is not clear what the effect would be on the posterior. Another avenue of
exploration is an HPP for time-to-event data. While the idea of eliciting a
prior prediction for the mean survival time is straightforward, it is not
clear how to handle right censored data in such a prior. In longitudinal data,
the development of an HPP may be particularly useful since $\bm{m}$ could be
elicited as a mean prediction within each cluster.

## References

  * Chen and Ibrahim [2003] Ming-Hui Chen and Joseph G. Ibrahim.  Conjugate priors for generalized linear models.  _Statistica Sinica_ , 13(2):461–476, 2003.  ISSN 1017-0405.  URL http://www.jstor.org/stable/24307144.  Publisher: Institute of Statistical Science, Academia Sinica.
  * Kuhnert et al. [2010] Petra M. Kuhnert, Tara G. Martin, and Shane P. Griffiths.  A guide to eliciting and using expert knowledge in Bayesian ecological models.  _Ecology Letters_ , 13(7):900–914, 2010.  ISSN 1461-0248.  doi:https://doi.org/10.1111/j.1461-0248.2010.01477.x.  URL http://onlinelibrary.wiley.com/doi/abs/10.1111/j.1461-0248.2010.01477.x.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1461-0248.2010.01477.x.
  * Banner et al. [2020] Katharine M. Banner, Kathryn M. Irvine, and Thomas J. Rodhouse.  The use of Bayesian priors in Ecology: The good, the bad and the not great.  _Methods in Ecology and Evolution_ , 11(8):882–889, 2020.  ISSN 2041-210X.  doi:https://doi.org/10.1111/2041-210X.13407.  URL http://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.13407.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/2041-210X.13407.
  * Grzenda [2016] Wioletta Grzenda.  Informative versus non-informative prior distributions and their impact on the accuracy of Bayesian inference.  _Statistics in Transition. New Series_ , 17:763–780, December 2016.  doi:10.21307/stattrans-2016-051.
  * Wang et al. [2018] Zheqi Wang, Jonathan N. Crook, and Galina Andreeva.  Improving Model Predictive Accuracy Using a Bayesian Approach: Application to PD Modelling of Mortgage Loans.  SSRN Scholarly Paper ID 3472200, Social Science Research Network, Rochester, NY, December 2018.  URL http://papers.ssrn.com/abstract=3472200.
  * Ibrahim and Chen [2000] Joseph G. Ibrahim and Ming-Hui Chen.  Power prior distributions for regression models.  _Statistical Science_ , 15(1):46–60, February 2000.  ISSN 0883-4237, 2168-8745.  doi:10.1214/ss/1009212673.  URL http://projecteuclid.org/euclid.ss/1009212673.  Publisher: Institute of Mathematical Statistics.
  * Hobbs et al. [2012] Brian P. Hobbs, Daniel J. Sargent, and Bradley P. Carlin.  Commensurate Priors for Incorporating Historical Information in Clinical Trials Using General and Generalized Linear Models.  _Bayesian analysis (Online)_ , 7(3):639–674, August 2012.  ISSN 1931-6690.  doi:10.1214/12-BA722.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4007051/.
  * Schmidli et al. [2014] Heinz Schmidli, Sandro Gsteiger, Satrajit Roychoudhury, Anthony O’Hagan, David Spiegelhalter, and Beat Neuenschwander.  Robust meta-analytic-predictive priors in clinical trials with historical control information.  _Biometrics_ , 70(4):1023–1032, 2014.  ISSN 1541-0420.  doi:10.1111/biom.12242.  URL http://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12242.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/biom.12242.
  * Rodhouse et al. [2019] Thomas J. Rodhouse, Rogelio M. Rodriguez, Katharine M. Banner, Patricia C. Ormsbee, Jenny Barnett, and Kathryn M. Irvine.  Evidence of region-wide bat population decline from long-term monitoring and Bayesian occupancy models with empirically informed priors.  _Ecology and Evolution_ , 9(19):11078–11088, 2019\.  ISSN 2045-7758.  doi:https://doi.org/10.1002/ece3.5612.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/ece3.5612.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/ece3.5612.
  * Diaconis and Ylvisaker [1979] Persi Diaconis and Donald Ylvisaker.  Conjugate Priors for Exponential Families.  _Annals of Statistics_ , 7(2):269–281, March 1979\.  ISSN 0090-5364, 2168-8966.  doi:10.1214/aos/1176344611.  URL http://projecteuclid.org/euclid.aos/1176344611.  Publisher: Institute of Mathematical Statistics.
  * Finney [1947] D. J. Finney.  The Estimation from Individual Records of the Relationship between Dose and Quantal Response.  _Biometrika_ , 34(3-4):320–334, December 1947\.  ISSN 0006-3444.  doi:10.1093/biomet/34.3-4.320.  URL https://doi.org/10.1093/biomet/34.3-4.320.
  * Ibrahim et al. [2015] Joseph G. Ibrahim, Ming-Hui Chen, Yeongjin Gwon, and Fang Chen.  The power prior: theory and applications.  _Statistics in Medicine_ , 34(28):3724–3749, 2015\.  ISSN 1097-0258.  doi:10.1002/sim.6728.  URL http://onlinelibrary.wiley.com/doi/abs/10.1002/sim.6728.
  * Zellner [1986] Arnold Zellner.  On Assessing Prior Distributions and Bayesian Regression Analysis with g-Prior Distributions.  In Prem Goel and Arnold Zellner, editors, _Bayesian Inference and Decision Techniques: Essays in Honor of Bruno De Finetti_ , volume 6, pages 233–243. Elsevier Science Ltd, June 1986.
  * Sabanés Bové and Held [2011] Daniel Sabanés Bové and Leonhard Held.  Hyper-$g$ priors for generalized linear models.  _Bayesian Analysis_ , 6(3):387–410, September 2011.  ISSN 1936-0975, 1931-6690.  doi:10.1214/11-BA615.  URL http://projecteuclid.org/journals/bayesian-analysis/volume-6/issue-3/Hyper-g-priors-for-generalized-linear-models/10.1214/11-BA615.full.  Publisher: International Society for Bayesian Analysis.
  * Held et al. [2015] Leonhard Held, Daniel Sabanés Bové, and Isaac Gravestock.  Approximate Bayesian model selection with the deviance statistic.  _Statistical Science_ , pages 242–257, 2015.  Publisher: JSTOR.
  * Bedrick et al. [1996] Edward J. Bedrick, Ronald Christensen, and Wesley Johnson.  A New Perspective on Priors for Generalized Linear Models.  _Journal of the American Statistical Association_ , 91(436):1450–1460, 1996.  ISSN 0162-1459.  doi:10.2307/2291571.  URL http://www.jstor.org/stable/2291571.  Publisher: [American Statistical Association, Taylor & Francis, Ltd.].
  * Hoffman and Gelman [2014] Matthew D Hoffman and Andrew Gelman.  The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo.  _Journal of Machine Learning Research_ , 15:1593–1623, April 2014.
  * Maechler et al. [2021] Martin Maechler, Peter Rousseeuw, Christophe Croux, Valentin Todorov, Andreas Ruckstuhl, Matias Salibian-Barrera, Tobias Verbeke, Manuel Koller, Eduardo L. T. Conceicao, and Maria Anna di Palma.  robustbase: Basic Robust Statistics, 2021.  URL http://robustbase.r-forge.r-project.org/.
  * Stan Development Team [2020] Stan Development Team.  RStan: the R interface to Stan, 2020.  URL http://mc-stan.org/.
  * Merigan et al. [1991] Thomas C. Merigan, David A. Amato, James Balsley, Maureen Power, William A. Price, Sharon Benoit, Angela Perez-Michael, Alan Brownstein, Amy Simon Kramer, Doreen Brettler, Louis Aledort, Margaret V. Ragni, W. Abe Andes, Joan C. Gill, Jonathan Goldsmith, Sally Stabler, Natalie Sanders, George Gjerset, and Jeanne Lusher.  Placebo-Controlled Trial to Evaluate Zidovudine in Treatment of Human Immunodeficiency Virus Infection in Asymptomatic Patients With Hemophilia.  _Blood_ , 78(4):900–906, August 1991.  ISSN 0006-4971.  doi:10.1182/blood.V78.4.900.900.  URL https://www.sciencedirect.com/science/article/pii/S0006497120841688.