# Bayesian inference for generalized linear models via quasi-posteriors

Davide Agnoletto  Department of Statistical Sciences, University of Padova,
Via C. Battisti 241, Padova, Italy  Tommaso Rigon  Department of Economics,
Management and Statistics, University of Milano–Bicocca, 20126 Milano, Italy
David B. Dunson  Department of Statistical Science, Duke University, Durham,
NC, USA

###### Abstract

Generalized linear models (glms) are routinely used for modeling relationships
between a response variable and a set of covariates. The simple form of a glm
comes with easy interpretability, but also leads to concerns about model
misspecification impacting inferential conclusions. A popular semi-parametric
solution adopted in the frequentist literature is quasi-likelihood, which
improves robustness by only requiring correct specification of the first two
moments. We develop a robust approach to Bayesian inference in glms through
quasi-posterior distributions. We show that quasi-posteriors provide a
coherent generalized Bayes inference method, while also approximating so-
called coarsened posteriors. In so doing, we obtain new insights into the
choice of coarsening parameter. Asymptotically, the quasi-posterior converges
in total variation to a normal distribution and has important connections with
the loss-likelihood bootstrap posterior. We demonstrate that it is also well-
calibrated in terms of frequentist coverage. Moreover, the loss-scale
parameter has a clear interpretation as a dispersion, and this leads to a
consolidated method of moments estimator.

Keywords: C-Bayes; Generalized Bayes; Model misspecification; Quasi-
likelihood; Robustness

##  1 Introduction

Generalized linear models (glms; Nelder and Wedderburn, 1972; McCullagh and
Nelder, 1989) are a canonical statistical framework for modeling a broad
variety of data. Despite their wide use in practice, they often incur
misspecification that could compromise inferential conclusions (Müller, 2013;
Grünwald and Van Ommen, 2017). For instance, assuming that the distribution of
the response variable belongs to a specific exponential family may be
untenable in many applications. A practical example regards the case of
heteroscedastic continuous data, as illustrated in Figure 1: in this situation
the homoscedasticity assumption of the standard Gaussian linear model is
wrong, resulting in a posterior which fails to concentrate around the “true”
parameter value. Another common case of misspecification occurs when binary or
count observations show larger variability than the one assumed by the model.
This is the so-called overdispersion issue, for which several model-based
solutions exist, such as the negative binomial model for count data or the
beta-binomial for binary data (Gelman et al., 2013). However, increased
complexity may lead to computational bottlenecks. Moreover, a model-based
solution can result again in misspecification problems. To bypass these
difficulties, one may adopt a nonparametric approach that does not make
assumptions about the underlying distribution of the data. On the one hand
this allows for greater flexibility, but it can also lead to an increased
computational cost and a loss of statistical efficiency and interpretability.

Figure 1: Comparison of the posterior distributions obtained via Gaussian
linear model and quasi-likelihood with $V(\mu_{i})=\exp\\{\mu_{i}\\}$ on a
synthetic dataset of $n=500$ heteroscedastic observations generated from a
Gaussian distribution with variance much larger than the mean. Left plot: data
points. Right plot: posterior distribution for the two-dimensional regression
parameter.

Here we rely on a semi-parametric approach, which represents a midpoint
between parametric and nonparametric methods, making weaker assumptions on the
distribution of the response variable while preserving computational
tractability. The quasi-likelihood method of Wedderburn (1974) makes second-
order assumptions, meaning that one only needs to specify the mean and
variance of the response variable. Let
$Y_{i}\in\mathcal{Y}\subseteq\mathds{R}$ denote a response variable and
$\bm{x}_{i}\in\mathcal{X}\subseteq\mathds{R}^{p}$ be a vector of covariates
for $i=1,\ldots,n$. Letting $\bm{\beta}\in\mathds{R}^{p}$ be the parameter of
interest, we assume that

$\begin{split}\mathds{E}(Y_{i})&=\mu_{i}=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}),\\\
\text{var}(Y_{i})&=\psi V(\mu_{i}),\end{split}$ (1)

where $\mu_{i}\in\mathcal{M}\subseteq\mathds{R}$, $g(\cdot)$ is a strictly
monotone and differentiable link function, $V(\cdot)>0$ is a continuous
variance function, and $\psi\in(0,\infty)$ is a dispersion parameter. The
flexibility brought by the second-order assumptions allows the quasi-
likelihood to handle the consequences of possible misspecification better than
model-based competitors. Indeed, as shown in the example of Figure 1, our
quasi-posterior adjusts for the substantial bias in the parametric glm
posterior, and includes the true parameter value in the center of the
distribution. Within the frequentist literature quasi-likelihood is well-
established. However, in the Bayesian setting such an approach is much more
complex to implement for conceptual and practical reasons. Indeed, within the
Bayesian paradigm, a statistician needs to specify all the probabilistic
aspects of the model, including those that one may want to leave unspecified.
This calls for generalized Bayesian paradigms, which may take into account
these practical considerations while maintaining a valid and formal
inferential procedure.

In the last decade, interest in generalized Bayesian methods has grown
considerably. The paper of Bissiri et al. (2016) represents a cornerstone in
this sense, showing that the update of prior beliefs via an exponentiated
negative loss in place of the likelihood function is valid and coherent, and
hence can be considered a genuine posterior distribution. The proportion of
information coming from the loss function and the prior is regulated by a
loss-scale parameter. Its calibration is delicate and, despite several
solutions (see Holmes and Walker, 2017; Lyddon et al., 2019), its choice
remains a difficult problem. Recently, Fong and Holmes (2020) used a
generalized Bayes approach to highlight connections between marginal
likelihood and cross-validation, while Rigon et al. (2023) developed methods
for clustering. Jewson and Rossell (2022) provide a general Bayesian selection
criterion for the loss function choice based on Hyvärinen score. For other
examples of generalized posteriors using loss functions, see Chernozhukov and
Hong (2003); Royall and Tsou (2003); Zhang (2006); Jiang and Tanner (2008);
Yin (2009).

Over the years a broad variety of generalized posteriors have been proposed
based on different notions of generalized likelihood functions. These include
partial likelihoods (Raftery et al., 1996; Sinha et al., 2003; Kim and Kim,
2009; Ventura and Racugno, 2016), composite likelihoods (Smith and Stephenson,
2009; Pauli et al., 2011; Ribatet et al., 2012), restricted likelihoods (Hoff,
2007), and empirical likelihoods (Lazar, 2003). These proposals can be placed
within the framework of Bissiri et al. (2016) by viewing the log generalized
likelihood as a loss. While the generalized Bayes paradigm focuses on Bayesian
updating without the requirement of a true model, the coarsened posterior
(c-posterior) framework by Miller and Dunson (2018) formalizes the problem of
having a small perturbation between the assumed model and the distribution of
the observed data. Finally, Grünwald and Van Ommen (2017) propose a tempered
likelihood solution for the misspecification problem in linear models.

In this paper, we focus on Bayesian modeling of the first two moments of the
unknown data generator without making any specific distributional assumption.
To do so, we rely on the second-order conditions (1) and we update the prior
beliefs on the parameter of interest through a generalized Bayes approach.
More precisely, inference will be based on the quasi-posterior distribution

$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)\propto
p(\bm{\beta})\exp\left\\{\sum_{i=1}^{n}\int_{a}^{\mu_{i}(\bm{\beta})}\frac{y_{i}-t}{\psi
V(t)}\text{d}t\right\\},$ (2)

where $\bm{y}=(y_{1},\ldots,y_{n})^{\top}$ is the vector of observed
responses, $\bm{X}=(\bm{x}_{1}^{\top},\ldots,\bm{x}_{n}^{\top})$ the matrix of
covariates, $a\in\mathcal{M}$ is an arbitrary constant that does not depend on
$\bm{\beta}$, and $p(\bm{\beta})$ the prior distribution for the parameter of
interest. There is a literature on related formulations to the quasi-posterior
(Lin, 2006; Greco et al., 2008; Ventura et al., 2010; Ventura and Racugno,
2016). However, our work provides novel contributions in formalizing a precise
generalized Bayes theory, showing connections with Gibbs posteriors (Bissiri
et al., 2016) and c-posteriors (Miller and Dunson, 2018), and proving several
new asymptotic properties. Moreover, although the starting assumptions are
similar, the quasi-posterior (2) is conceptually and practically different
from the Bayesian generalized methods of moments of Yin (2009).

Our quasi-posterior takes the form of the c-posterior of Miller and Dunson
(2018) when $V(\mu_{i})$ in (1) corresponds to the variance function of an
exponential family model. In this case, our quasi-posterior based on second-
order assumptions is equivalent to a fully specified model affected by small
perturbations. This result provides a new perspective on the coarsening
parameter of Miller and Dunson (2018), which we show to be directly connected
with the dispersion parameter $\psi$. This insight has relevant practical
implications: c-posteriors have become a popular approach for robustifying
Bayesian inferences to model misspecification but there are only heuristic
algorithms available for choosing the amount of coarsening. In contrast, we
propose a simple method of moments estimator. Even when the above variance
function condition is violated, the quasi-posterior is still interpretable as
a Gibbs posterior and we can use the same method of moments estimator in
choosing the $\psi$ parameter. Moreover, the quasi-posterior has desirable
asymptotic properties, since it converges to a normal distribution centered
around the true parameter value and shows correct frequentist coverage
asymptotically.

##  2 Bayesian inference for GLMs under second-order assumptions

###  2.1 Background on quasi-likelihood

Let $\bm{Y}=(Y_{1},\ldots,Y_{n})^{\top}\in\mathcal{Y}^{n}$ be a vector of
random variables whose realization is $\bm{y}$. Let
$\bm{x}_{i}\in\mathcal{X}$, $i=1,\ldots,n$, be vectors of linear independent
covariates such that $\bm{X}$ has full rank for any $n\geq p$. Standard glms
assume that observations $y_{i}$ are independent realizations of
$Y_{i}\mid\bm{x}_{i}$, whose distribution belongs to the exponential
dispersion family. Here, we rely on the second-order assumptions (1), which
specify only the first two moments instead of the entire distribution for
$Y_{i}\mid\bm{x}_{i}$. We remark on the fact that under the second-order
conditions we are assuming the existence of a true value for
$(\bm{\beta},\psi)$, say $(\bm{\beta}_{0},\psi_{0})$, and consequently the
existence of a true probability measure
${F_{0}(\mathrm{d}y\mid\bm{x})=F(\mathrm{d}y\mid\bm{x},\bm{\beta}_{0},\psi_{0})}$
on $\mathcal{Y}$ and a joint measure $F_{0}(\mathrm{d}\bm{y}\mid\bm{X})$ on
$\mathcal{Y}^{n}$. However, in this case, we only make assumptions about the
form of mean and variance of the distribution.

A frequentist semiparametric solution for inference under second-order
assumptions was introduced by Wedderburn (1974). Let

$\frac{1}{\psi}\sum_{i=1}^{n}\frac{y_{i}-\mu_{i}}{V(\mu_{i})}\frac{\partial\mu_{i}}{\partial\beta_{r}},$
(3)

be a $p$-dimensional combinant, for $r=1,\ldots,p$. It coincides with the
general form of the score function for an exponential dispersion family and it
is an unbiased estimating equation. Moreover, when the derivative of (3) with
respect to $\bm{\beta}$ is a symmetric matrix, it is possible to specify a
function reminiscent of a log-likelihood function (known as log-quasi-
likelihood) which has (3) as gradient but does not necessarily correspond to
any distribution. In particular, Wedderburn (1974) proposed

$\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)=\sum_{i=1}^{n}\ell_{\textsc{q}}(\bm{\beta};y_{i},\bm{x}_{i},\psi)=\frac{1}{\psi}\sum_{i=1}^{n}\int_{a}^{\mu_{i}(\bm{\beta})}\frac{y_{i}-t}{V(t)}\mathrm{d}t.$
(4)

Variance function $V(\cdot)$ must be chosen so that each integral in (4) is
finite. Log-quasi-likelihood is equally suitable for both continuous and
discrete data. The integral can be written in closed form for several choices
of variance function, including those corresponding to popular distributions
such as Gaussian, Poisson, and binomial, as shown in Table 1. A detailed
discussion of the properties of quasi-likelihood function is contained in
McCullagh (1983).

$V(\mu)$ | $\ell_{\textsc{q}}(\mu;y,\psi)$ | Restrictions
---|---|---
$1$ | $\frac{1}{\psi}(y\mu-\frac{1}{2}\mu^{2})$ | $-$
$\mu$ | $\frac{1}{\psi}(y\log\mu-\mu)$ | $\mu>0,\;y\geq 0$
$\mu^{2}$ | $\frac{1}{\psi}(-y/\mu-\log\mu)$ | $\mu>0,\;y\geq 0$
$\mu^{p}$ | $\frac{1}{\psi}\big{[}(1-p)^{-1}y\mu^{1-p}-(2-p)^{-1}\mu^{2-p}\big{]}$ | $p>2,\;\mu>0,\;y\geq 0$
$e^{\mu}$ | $\frac{1}{\psi}(\mu-y+1)e^{-\mu}$ | $-$
$\mu(1-\mu)$ | $\frac{1}{\psi}\big{[}y\log\mu+(1-y)\log(1-\mu)\big{]}$ | $\mu\in(0,1),\;y\in[0,1]$
$\mu^{2}(1-\mu)^{2}$ | $\frac{1}{\psi}\big{[}(2y-1)\log\frac{\mu}{1-\mu}+\frac{2y-1}{1-\mu}-\frac{y}{\mu(1-\mu)}\big{]}$ | $\mu\in(0,1),\;y\in[0,1]$
$\mu+\mu^{2}/k$ | $\frac{1}{\psi}\big{(}y\log\frac{\mu}{k+\mu}+k\log\frac{1}{k+\mu}\big{)}$ | $\mu>0,\;y\geq 0$
Table 1: Closed form of the log-quasi-likelihood function for several choices
of variance functions.

###  2.2 Quasi-posterior

In classical Bayesian statistics, inference about $\bm{\beta}$ is carried out
using the posterior distribution, which updates the prior distribution with
the likelihood function. Under only the second-order assumptions, without
additional conditions, a genuine likelihood is not available. However, a
posterior distribution for $\bm{\beta}$ can be obtained using the quasi-
likelihood.

######  Definition 1.

Let $\exp\\{\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)\\}$ be the quasi-
likelihood function and $p(\bm{\beta})$ be the prior distribution for
$\bm{\beta}$. We refer to

$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)\propto
p(\bm{\beta})\exp\left\\{\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)\right\\}=p(\bm{\beta})\exp\left\\{\frac{1}{\psi}\sum_{i=1}^{n}\int_{a}^{\mu_{i}(\bm{\beta})}\frac{y_{i}-t}{V(t)}\text{d}t\right\\}$

as the quasi-posterior distribution for $\bm{\beta}$.

The quantity $p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)$ is a well-
defined probability density function provided that the normalizing constant is
such that $0<\int
p(\bm{\beta})\exp\left\\{\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)\right\\}\mathrm{d}\bm{\beta}<\infty$.
Although this condition is not always true in general, it can be easily
verified that it holds for a wide range of variance functions, including those
most commonly used in practice.

The quasi-posterior is a rational update of a belief distribution under a
generalized Bayesian model (Bissiri et al., 2016). In a regression context,
the target of inference for a generalized Bayes procedure is
$\bm{\beta}_{\textsc{opt}}$, defined as the minimizer of an integrated loss
function,

$\bm{\beta}_{\textsc{opt}}=\arg\min_{\bm{\beta}}\int\ell(\bm{\beta};\bm{y},\bm{X})F_{0}(\mathrm{d}\bm{y}\mid\bm{X}),$
(5)

where $\ell(\bm{\beta};\bm{y},\bm{X})$ is a loss function and
$\bm{\beta}_{\textsc{opt}}$ is assumed to be unique. Bissiri et al. (2016)
used a decision-theoretic approach to determine a coherent update of the prior
beliefs about $\bm{\beta}_{\textsc{opt}}$ without knowledge of the likelihood
function. Letting $\mathcal{P}$ be the space of all conditional distributions
of $\bm{\beta}$ given the data, the following loss function is used to measure
how well $\nu\in\mathcal{P}$ represents posterior beliefs:

$\mathscr{L}\\{\nu(\bm{\beta})\\}=\frac{1}{\psi}\mathds{E}\left\\{\ell(\bm{\beta};\bm{y},\bm{X})\right\\}+\textsc{kl}\\{\nu(\bm{\beta})\mid\mid
p(\bm{\beta})\\},$

where the expectation is taken with respect to $\nu(\bm{\beta})$ and kl
represents the Kullback-Leibler divergence. The first term quantifies the
discrepancy between $\nu$ and the information about the parameter coming from
the observed data, while the second quantifies the closeness to the prior. The
loss-scale parameter $\psi>0$ regulates the balance between these two terms.
The only coherent update of the prior beliefs $p(\bm{\beta})$ about
$\bm{\beta}_{\textsc{opt}}$, after observing the data, to the posterior
$p_{\psi}(\bm{\beta}\mid\bm{y},\bm{X})$ is

$p_{\psi}(\bm{\beta}\mid\bm{y},\bm{X})=\arg\min_{\nu}\mathscr{L}\\{\nu(\bm{\beta})\\}\propto
p(\bm{\beta})\exp\left\\{-\frac{1}{\psi}\ell(\bm{\beta};\bm{y},\bm{X})\right\\},$
(6)

where
$\ell(\bm{\beta};\bm{y},\bm{X})=\sum_{i=1}^{n}\ell(\bm{\beta};y_{i},\bm{x}_{i})$
is an additive loss function. It is straightforward to see that the quasi-
posterior of Definition 1 corresponds to the minimizer (6) with

$\ell(\bm{\beta};\bm{y},\bm{X})=-\sum_{i=1}^{n}\int_{a}^{\mu_{i}(\bm{\beta})}\frac{y_{i}-t}{V(t)}\mathrm{d}t.$
(7)

Taking a generalized Bayesian approach, the target of quasi-posterior
inference is $\bm{\beta}_{\textsc{opt}}$ defined in (5). The second-order
conditions assume the existence of a true value $\bm{\beta}_{0}$ for the
regression parameter. The relation between this two quantities is clarified by
the following Theorem.

######  Theorem 1.

For quasi-posteriors, suppose the target of inference
$\bm{\beta}_{\textsc{opt}}$ defined in equation (5) is unique. Then, it must
coincide with the true value $\bm{\beta}_{0}$ assumed under the second-order
conditions

$\bm{\beta}_{\textsc{opt}}=\arg\max_{\bm{\beta}}\int\ell_{Q}(\bm{\beta};\bm{y},\bm{X},\psi)F_{0}(\mathrm{d}\bm{y}\mid\bm{X})=\bm{\beta}_{0}.$

The crucial assumption of the Bissiri et al. (2016) framework is that
$\bm{\beta}_{\textsc{opt}}$ should represent a quantity of practical interest.
Theorem 1 shows that this is indeed the case for our quasi-posterior - the
quantity targeted in the generalized Bayesian paradigm coincides with the true
parameter value. Hence, the quasi-posterior represents a worthwhile
representation of posterior beliefs about $\bm{\beta}$.

##  3 Connection with c-posteriors

C-posteriors provide an alternative approach for robustifying Bayesian
inferences to model misspecification. In this section, we formally connect our
quasi-posterior with c-posteriors, which are briefly reviewed here. Start by
choosing a parametric model for the data generator such that
$Y_{i}\mid\bm{x}_{i}$, $i=1,\ldots,n$, are drawn from the density function
$f(y\mid\bm{\beta},\bm{x})$. The observed responses $\bm{y}$ are assumed to be
perturbed realizations of $\bm{Y}$, in the sense that
${d(\bm{Y},\bm{y}\mid\bm{X})<R}$ for some measure of discrepancy
$d(\cdot,\cdot\mid\cdot)$ and some $R\in(0,\infty)$. Suppose that the observed
data $y_{i}\mid\bm{x}_{i}$, $i=1,\ldots,n$, behave like independent samples
from some unknown perturbed density $\tilde{f}(y\mid\bm{x})$. Inference on
$\bm{\beta}$ using the standard posterior is equivalent to conditioning on the
event $\bm{Y}=\bm{y}$. However, this is incorrect when there is a
perturbation. With this motivation, Miller and Dunson (2018) propose to infer
$\bm{\beta}$ conditioning on the event $d(\bm{Y},\bm{y}\mid\bm{X})<R$. The
resulting posterior distribution is the so-called c-posterior, defined as
${p\left(\bm{\beta}\mid d(\bm{Y},\bm{y}\mid\bm{X})<R\right)}$. It is important
to clarify we are not assuming that $\tilde{f}(\mathrm{d}y\mid\bm{x})$
corresponds to the density $f_{0}(y\mid\bm{x})$ of
$F_{0}(\mathrm{d}y\mid\bm{x})$: the first is the perturbed law that generates
the observations according to Miller and Dunson (2018) approach, while the
second is the unknown true mechanism of which we are assuming the first two
moments.

A reasonable choice of $d(\cdot,\cdot\mid\cdot)$ is the estimator of the
conditional relative entropy between $\tilde{f}(y\mid\bm{x})$ and
$f(y\mid\bm{\beta},\bm{x})$, that is
$d(\bm{Y},\bm{y}\mid\bm{X})=n^{-1}\sum_{i=1}^{n}\big{[}\log\tilde{f}(y_{i}\mid\bm{x}_{i})-\log
f(y_{i}\mid\bm{\beta},\bm{x}_{i})\big{]}$. Then, if $R\sim\text{Exp}(\alpha)$,
the regression c-posterior can be approximated with a power posterior of the
form

$p_{\alpha}(\bm{\beta}\mid\bm{y},\bm{X})\propto
p(\bm{\beta})\prod_{i=1}^{n}f(y_{i}\mid\bm{\beta},\bm{x}_{i})^{\frac{\alpha}{\alpha+n}},$

where $\alpha/(\alpha+n)$ represents the amount of coarsening. The choice of
$\alpha$ is crucial since it regulates the amount of coarsening and
consequently the concentration of the posterior: intuitively, to raise a
likelihood to the power of $\alpha/(\alpha+n)$ corresponds to adjusting the
sample size from $n$ to $n\alpha/(\alpha+n)$, so the c-posterior is as
concentrated as a regular posterior having sample size $n\alpha/(\alpha+n)$.
When $\alpha\gg n$, the power posterior will be similar to the standard
posterior, while the amount of coarsening tends to $\alpha/n$ if $\alpha\ll
n$. Miller and Dunson (2018) propose heuristics for choosing $\alpha$,
including simply letting $\alpha=n$, and a graphical criterion based on a
calibration curve which compares, for a fixed $\alpha$, the posterior expected
log-likelihood, chosen as a measure of goodness of fit, and the posterior
expected model complexity. However, current approaches for choosing $\alpha$
lack theoretical support.

A solution can be found in the connection between quasi- and c-posteriors.
Indeed, the power posterior approximation of the c-posterior is recovered as a
particular case of the quasi-posterior when the assumed model is an
exponential family. The crucial point is the relation between the quasi-
likelihood and the genuine likelihood for the exponential family, which is
recalled in the following Theorem.

######  Theorem 2 (Wedderburn, 1974).

For a single observation $Y$, the log-likelihood function has score function
$(y-\mu)V(\mu)^{-1}$, where $\mu=\mathds{E}(Y)$ and $V(\mu)=\mathrm{var}(Y)$,
if and only if the density of $Y$ with respect to some measure can be written
in the exponential form $\exp\\{y\theta-b(\theta)\\}$, where
$\theta=\theta(\mu)$ is the so-called natural parametrization.

In other words, the genuine log-likelihood of an exponential family coincides
with the corresponding log-quasi-likelihood, up to a multiplicative factor.
Table 1 contains the examples for Gaussian, Poisson, and binomial
distributions corresponding to $V(\mu)=1$, $V(\mu)=\mu$, and
$V(\mu)=\mu(1-\mu)$, respectively. The connection between quasi- and
c-posterior arises naturally and is summarized in the following Proposition.

######  Proposition 1.

Assume the conditional distribution of $Y_{i}\mid\bm{x}_{i}$, $i=1,\ldots,n$,
to be an exponential family of order 1, that is

$f(y_{i}\mid\bm{\beta},\bm{x}_{i})=h(y_{i})\exp\left\\{\theta(\mu_{i})y_{i}-b(\theta(\mu_{i}))\right\\},$
(8)

where $\theta(\cdot)$ is the natural parametrization,
$\mu_{i}=\mu_{i}(\bm{\beta})=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta})$, and
$h(y_{i})$ the normalizing constant. If $V(\cdot)$ in (1) corresponds to the
variance function of (8), then the corresponding quasi-posterior is equivalent
to the power posterior approximation of the c-posterior based on (8), with
$\psi=(\alpha+n)/\alpha$.

It follows directly that the quasi-posterior is approximately proportional to
the regression c-posterior distribution based on conditional relative entropy.
This relation links posteriors based on two different notions of model
misspecification: using log-quasi-likelihood as a loss function within a Gibbs
posterior (Bissiri et al., 2016) only requires specification of the first two
moments, while the c-posterior requires a complete specification of the
likelihood of the data but allows this likelihood to be slightly misspecified.
The above Proposition unifies these two approaches for robust Bayesian
inference for glms. This unification highlights the relation between amount of
coarsening and dispersion in the data, allowing estimation of $\alpha$ based
on $\psi$ as discussed in Section 4.2. The following Example clarifies the
link between $\alpha$ and $\psi$ for count data.

######  Example 1 (Overdispersed Poisson regression).

Consider the relative entropy regression c-posterior for count response with
the assumed model a Poisson distribution with canonical link function. The
corresponding power posterior approximation is

$p_{\alpha}(\bm{\beta}\mid\bm{y},\bm{X})\propto
p(\bm{\beta})\exp\left\\{\frac{\alpha}{\alpha+n}\sum_{i=1}^{n}\left(y_{i}\bm{x}_{i}^{\top}\bm{\beta}-\exp\\{\bm{x}_{i}^{\top}\bm{\beta}\\}\right)\right\\}.$

Assuming $\mathds{E}(Y_{i})=\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$ and
$\mathrm{var}(Y_{i})=\psi\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$, the
quasi-likelihood (4) is
$\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)\propto\psi^{-1}\sum_{i=1}^{n}\big{(}y_{i}\bm{x}_{i}^{\top}\bm{\beta}-\exp\\{\bm{x}_{i}^{\top}\bm{\beta}\\}\big{)}$,
which leads to the quasi-posterior

$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)\propto
p(\bm{\beta})\exp\left\\{\frac{1}{\psi}\sum_{i=1}^{n}\left(y_{i}\bm{x}_{i}^{\top}\bm{\beta}-\exp\\{\bm{x}_{i}^{\top}\bm{\beta}\\}\right)\right\\}.$
(9)

The two posteriors are equivalent when

$\alpha=\frac{n}{\psi-1}.$

This relation offers some insights into the phenomenon that c-posteriors and
quasi-posteriors can be robust against. In this example, $\psi$ controls the
difference between the dispersion of the observations and the variance assumed
by the parametric model: for $\psi=1$ the Poisson model is well specified,
while if $\psi>1$ or $\psi\in(0,1)$ it is affected by overdispersion or
underdispersion, respectively. Both the quasi-posterior and c-posterior are
robust against overdispersion, since $\psi>1$ implies $\alpha\in(0,\infty)$,
and they recover the standard posterior when $\psi\rightarrow 1$ and so
$\alpha\rightarrow\infty$. However, the quasi-posterior is also robust against
underdispersion while the c-posterior is not, since $\alpha>0$ by
construction.

##  4 Asymptotic properties

In this section we study large sample properties of quasi-posteriors. The
literature contains some related work on generalized posteriors; for example,
studying convergence in distribution (see Greco et al., 2008; Ventura and
Racugno, 2016) or in total variation (Chernozhukov and Hong, 2003) to a normal
distribution. Here we provide results on limiting distributions and
frequentist coverage of quasi posteriors. In the following theorem, based on
Theorems 4 and 5 of Miller (2021), we present the asymptotic behavior of the
quasi-posterior distribution.

######  Theorem 3.

Let $p(\bm{\beta})$ be continuous at $\bm{\beta}_{0}$ and such that
$p(\bm{\beta}_{0})>0$. Let $Y_{i}\mid\bm{x}_{i}$, $i=1,\ldots,n$ be
independent random variables generated according to
$F_{0}(\mathrm{d}y\mid\bm{x})$ and assume second-order conditions of equation
(1) are satisfied. Let
$\hat{\bm{\beta}}=\arg\max_{\bm{\beta}}\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$
be the empirical risk minimizer and
$i(\bm{\beta})=\lim_{n\rightarrow\infty}(\psi
n)^{-1}\bm{X}^{\top}\bm{D}\bm{X}$ where
$\bm{D}=\mathrm{diag}(d_{1},\ldots,d_{n})$,
$d_{i}=[V(\mu_{i})g^{\prime}(\mu_{i})^{2}]^{-1}$, for $i=1,\ldots,n$. Then,
under regularity conditions, the density $p_{n}(\bm{\beta})$ of the centered
and scaled quasi-posterior $\sqrt{n}(\bm{\beta}-\hat{\bm{\beta}})$ is such
that

$\left|p_{n}(\bm{\beta})-\mathcal{N}_{p}\left(\bm{\beta}\mid
0,i(\bm{\beta}_{0})^{-1}\right)\right|_{\text{TV}}\overset{F_{0}}{\longrightarrow}0,\qquad
n\rightarrow\infty,$

almost surely. Moreover, the estimator $\hat{\bm{\beta}}$ is consistent, that
is $\hat{\bm{\beta}}\rightarrow\bm{\beta}_{0}$ almost surely.

The regularity conditions are discussed in Appendix A.3. Theorem 3 shows that
the quasi-posterior converges in total variation to a Gaussian distribution
centered around the maximum quasi-likelihood estimator $\hat{\bm{\beta}}$,
which is a consistent estimator of $\bm{\beta}_{0}$. Hence, the quasi-
posterior for $\bm{\beta}$ asymptotically concentrates around the true
parameter value. The quantity $\hat{\bm{\beta}}$ is the solution of the
unbiased estimating equation obtained with (3), and is an M-estimator.
Furthermore, when (3) coincides with the score function of an exponential
family, it is equivalent to the maximum likelihood estimator (mle). As a
consequence of Theorem 1, $\hat{\bm{\beta}}$ is also the empirical counterpart
of $\bm{\beta}_{\textsc{opt}}$, often called the empirical risk minimizer. The
covariance of the asymptotic distribution is determined by the expected
information matrix

$I(\bm{\beta})=\mathds{E}\left[-\nabla^{2}\ell_{\textsc{q}}(\bm{\beta};\bm{Y},\bm{X},\psi)\right]=\psi^{-1}\bm{X}^{\top}\bm{D}\bm{X},$
(10)

where $\nabla$ denotes the gradient operator with respect to $\bm{\beta}$,
computed for $\bm{\beta}=\bm{\beta}_{0}$. Consequently, $\psi$ governs the
concentration of quasi-posterior around $\hat{\bm{\beta}}$: as $\psi$ grows,
the asymptotic variance is inflated, forcing the quasi-posterior to be
flatter. Theorem 3 is stated for a general form of quasi-likelihood, but more
specific conclusions for several popular choices of link function and variance
function could be obtained, e.g.

  1. 1.

$\mathds{E}(Y_{i})=\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$,
$\mathrm{var}(Y_{i})=\psi_{0}\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$
(overdispersed Poisson regression);

  2. 2.

$\mathds{E}(Y_{i})=\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}(1+\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\})^{-1}=\mu_{i}$,
$\mathrm{var}(Y_{i})=\psi_{0}\mu_{i}(1-\mu_{i})$ (overdispersed logistic
regression);

  3. 3.

$\mathds{E}(Y_{i})=\bm{x}_{i}^{\top}\bm{\beta}_{0}$,
$\mathrm{var}(Y_{i})=\psi_{0}\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$
(heteroscedastic continuous data regression).

While for the overdispersed Poisson regression and overdispersed logistic
regression the result of Theorem 3 may not be surprising since the derivatives
are equivalent to the parametric case up to a constant, it remarkably holds,
under mild regularity conditions, even for the case of heteroscedastic
continuous data regression.

######  Example 2 (Heteroscedastic continuous data regression).

Assume $\mathds{E}(Y_{i})=\bm{x}_{i}^{\top}\bm{\beta}_{0}$ and
$\mathrm{var}(Y_{i})=\psi_{0}\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$,
$i=1,\ldots,n$. The corresponding quasi-posterior distribution is

$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)\propto
p(\bm{\beta})\exp\left\\{\frac{1}{\psi}\sum_{i=1}^{n}(\bm{x}_{i}^{\top}\bm{\beta}-y_{i}+1)\exp\\{-\bm{x}_{i}^{\top}\bm{\beta}\\}\right\\}.$

Thus, the quasi posterior distribution of $\bm{\beta}$ is approximately a
$p$-variate Gaussian with mean $\hat{\bm{\beta}}$ and covariance matrix
$I(\bm{\beta}_{0})^{-1}$, with
$I(\bm{\beta}_{0})=\psi^{-1}\bm{X}^{\top}\tilde{\bm{D}}\bm{X}$ and
$\tilde{\bm{D}}=\mathrm{diag}(\tilde{d}_{1},\ldots,\tilde{d}_{n})$, where
$\tilde{d}_{i}=\exp\\{-\bm{x}_{i}^{\top}\bm{\beta}_{0}\\}$, for
$i=1,\ldots,n$. The matrix $I(\bm{\beta}_{0})$ can be practically evaluated
using the plug-in estimate $I(\hat{\bm{\beta}})$.

###  4.1 Correct frequentist coverage

A notable property of the asymptotic distribution of the quasi-posterior is
correct frequentist coverage. We say that the quasi-posterior is well-
calibrated in terms of frequentist coverage when, asymptotically, the
posterior credible sets of level $\rho$ contain the true parameter value
$\bm{\beta}_{0}$ with probability $\rho$. Generalized posteriors do not
automatically have this property, so quasi-posteriors are special in this
regard. Miller (2021) provides sufficient conditions for which a generalized
posterior is asymptotically calibrated in terms of coverage. Roughly speaking,
a generalized posterior for $\bm{\beta}$ has correct coverage if converges in
total variation distance to the same normal distribution to which the
empirical risk minimizer $\hat{\bm{\beta}}$ converges in distribution.
Although this condition is not satisfied for many common choices of
generalized posteriors, it holds for quasi-posteriors. This allows the quasi-
posteriors to have asymptotically correct frequentist coverage, thanks to the
general results of Miller (2021). The proof is discussed in Appendix A.4,
together with a description of the regularity conditions.

######  Theorem 4.

Let $S_{1},S_{2},\ldots\subseteq\mathds{R}^{p}$ be a sequence of convex
credible sets for the quasi-posterior of $\bm{\beta}$ of asymptotic level
$\rho\in(0,1)$. Then, under mild regularity conditions and assuming
$\psi=\psi_{0}$

$\mathds{P}(\bm{\beta}_{0}\in S_{n}\mid\bm{y},\bm{X},\psi_{0})\rightarrow\rho$

as $n\rightarrow\infty.$

###  4.2 Connection with loss-likelihood bootstrap

The quasi-posterior has connections with the loss-likelihood bootstrap
introduced by Lyddon et al. (2019) and their proposal for setting the value of
the loss-scale parameter. The loss-likelihood bootstrap is a version of the
Bayesian bootstrap (Rubin, 1981) which, essentially, enables generation of
samples from an approximate posterior distribution built from a loss function.
Asymptotically, the loss-likelihood bootstrap posterior converges to a normal
distribution with covariance matrix that does not depend on $\psi$.
Conversely, the asymptotic normal distribution of a generalized posterior
having the same loss function depends on the loss-scale parameter; an example
is the result of Theorem 3 for the quasi-posterior. For this reason, Lyddon et
al. (2019) proposed to choose $\psi$ by matching the amount of information of
the two asymptotic distributions quantified in terms of the Fisher information
number, that, for a normal distribution, is asymptotically equal to the trace
of the Fisher information matrix. The matching value is

$\psi_{\textsc{llb}}=\frac{\text{tr}\\{j(\bm{\beta}_{0})\\}}{\text{tr}\\{j(\bm{\beta}_{0})h(\bm{\beta}_{0})^{-1}j(\bm{\beta}_{0})\\}},$
(11)

where

$j(\bm{\beta})=\lim_{n\rightarrow\infty}n^{-1}\mathds{E}\left[\nabla^{2}\ell(\bm{\beta};\bm{Y},\bm{X})\right],\qquad
h(\bm{\beta})=\lim_{n\rightarrow\infty}n^{-1}\mathds{E}\left[\nabla\ell(\bm{\beta};\bm{Y},\bm{X})\nabla\ell(\bm{\beta};\bm{Y},\bm{X})^{\top}\right].$

This result refers to a generic loss function $\ell(\cdot)$. However, for
quasi-posteriors we obtain a substantial simplification.

######  Proposition 2.

If
$\mathds{E}(Y_{i})=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})=\mu_{i}(\bm{\beta}_{0})$
and $\mathrm{var}(Y_{i})=\psi_{0}V(\mu_{i}(\bm{\beta}_{0}))$, with both
operators taken with respect to $F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})$, and
loss function in equation (7), then $\psi_{\textsc{llb}}=\psi_{0}$.

This means that under the second order conditions, the value of $\psi$ that
equates the Fisher information number of the two asymptotic distributions is
the true value $\psi_{0}$. In practice, of course, $\psi_{0}$ is unknown and
the common usage in the generalized Bayes literature is to assign a value to
$\psi$ based on the data before approximate sampling from the generalized
posterior (see Bissiri et al., 2016; Holmes and Walker, 2017; Rigon et al.,
2023). For example, Lyddon et al. (2019) propose the estimator
${\hat{\psi}_{\textsc{llb}}=\text{tr}\\{j_{n}(\hat{\bm{\beta}})\\}\text{tr}\\{j_{n}(\hat{\bm{\beta}})h_{n}(\hat{\bm{\beta}})^{-1}j_{n}(\hat{\bm{\beta}})\\}^{-1}}$,
where
${j_{n}(\bm{\beta})=n^{-1}\sum_{i=1}^{n}\nabla^{2}\ell(\bm{\beta};y_{i},\bm{x}_{i})}$
and
$h_{n}(\bm{\beta})=n^{-1}\sum_{i=1}^{n}\nabla\ell(\bm{\beta};y_{i},\bm{x}_{i})\nabla\ell(\bm{\beta};y_{i},\bm{x}_{i})^{\top}$
are the empirical versions of $j(\bm{\beta})$ and $h(\bm{\beta})$,
respectively. However, it is not guaranteed that $\hat{\psi}_{\textsc{llb}}$
is a consistent estimator of $\psi_{0}$. For this reason we suggest to use the
estimator based on the method of moments

$\hat{\psi}=\frac{1}{n-p}\sum_{i=1}^{n}\frac{\big{(}y_{i}-\hat{\mu}_{i}\big{)}^{2}}{V\big{(}\hat{\mu}_{i}\big{)}},$
(12)

where $\hat{\mu}_{i}=\mu_{i}(\hat{\bm{\beta}})$. Equation (12) produces a
consistent estimator of the true dispersion parameter and it represents the
standard choice in the frequentist quasi-likelihood literature (Wedderburn,
1974; McCullagh, 1983; McCullagh and Nelder, 1989). An exemplification for the
comparison of $\hat{\psi}$ and $\hat{\psi}_{\textsc{llb}}$ is the Gaussian
loss with $\beta\in\mathds{R}$ illustrated in Example 3.

######  Example 3.

Let $\mathds{E}(Y_{i})=x_{i}\beta$ and $\mathrm{var}(Y_{i})=\psi$, for
$i=1,\ldots,n$. Then the method of moments estimator for $\psi$ is the well
known variance estimator
${\hat{\psi}=(n-1)^{-1}\sum_{i=1}^{n}(y_{i}-x_{i}\hat{\beta})^{2}}$, while the
one based on loss-likelihood bootstrap is
${\hat{\psi}_{\textsc{llb}}=(\sum_{i=1}^{n}x_{i}^{2})^{-1}\sum_{i=1}^{n}x_{i}^{2}(y_{i}-x_{i}\hat{\beta})^{2}}$.
However, while $\hat{\psi}$ is unbiased, $\hat{\psi}_{\textsc{llb}}$ is not,
since
$\mathds{E}(\hat{\psi}_{\textsc{llb}})=\psi\big{[}1-\sum_{i=1}^{n}x_{i}^{4}/(\sum_{i=1}^{n}x_{i}^{2})^{2}\big{]}$.

A further justification for (12) is provided by the coverage result of Theorem
4. Indeed, a condition for having correct coverage is that the quasi-posterior
and empirical risk minimizer have the same asymptotic covariance matrix, which
depends on $\psi$. Since the standard frequentist choice for its estimation is
(12), the method of moments estimator represents the solution that keeps this
equivalence valid in practice.

Importantly, $\hat{\bm{\beta}}$ does not depend on $\psi$. Thus, by estimating
$\psi$ using the empirical risk minimizer, we avoid need for an iterative
algorithm. The estimated value can be interpreted as the observed amount of
dispersion, based on the value of $\bm{\beta}$ that minimizes the quasi-
likelihood function in the sample.

##  5 Illustrations

Quasi-posteriors are a useful tool in many practical situations. In this
section we illustrate some simulated and real data examples of quasi-
posteriors in different contexts where traditional model-based approaches
struggle.

###  5.1 Simulation example: heteroscedastic continuous data

To show consequences of model misspecification on estimates of regression
coefficients, we apply a Gaussian linear model to heteroscedastic continuous
data and compare performance to a quasi-posterior with mean and variance
correctly specified. We produce 100 simulated datasets, each of size $n=300$.
In each simulation, the true data generator for $y_{i}\mid\bm{x}_{i}$ is
$\mathcal{N}(\bm{x}_{i}^{\top}\bm{\beta}_{0},\psi_{0}\exp\\{\bm{x}_{i}^{\top}\bm{\beta}_{0}\\})$,
$i=1,\ldots,n$, with true parameter value $\bm{\beta}_{0}=(-3,2,1.5,1)$ and
$\psi_{0}=2.5$. The covariates $\bm{x}_{i}$ are sampled from independent
standard normal distributions.

Table 2: Heteroscedastic continuous data example. Frequentist coverage of credible intervals with nominal level $\rho\in\\{0.90,0.95,0.99\\}$ for Gaussian linear model (lm) and quasi-posterior (qp). |  | $\rho=0.90$ |  | $\rho=0.95$ |  | $\rho=0.99$
---|---|---|---|---|---|---
|  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$
lm |  | 0.92 | 0.65 | 0.71 | 0.76 |  | 0.95 | 0.71 | 0.82 | 0.84 |  | 0.99 | 0.86 | 0.93 | 0.93
qp |  | 0.94 | 0.91 | 0.90 | 0.93 |  | 0.96 | 0.95 | 0.96 | 0.95 |  | 0.99 | 1.00 | 0.99 | 1.00
Figure 2: Heteroscedastic continuous data example. Distribution of posterior
mean of $\bm{\beta}$ obtained with Gaussian linear model and quasi-posterior.

The log-quasi-likelihood function contribution for the $i$th observation is

${\ell_{\textsc{q}}(\bm{\beta};y_{i},\bm{x}_{i},\psi)=\psi^{-1}(\bm{x}_{i}^{\top}\bm{\beta}-y_{i}+1)\exp\\{-\bm{x}_{i}^{\top}\bm{\beta}\\}},$

which leads to the quasi-posterior presented in Example 2. In both the
approaches we put an uninformative prior distribution on $\bm{\beta}$, and for
the normal linear model an uninformative prior is put on the variance. For the
quasi-posterior, coherently with Section 4.2, we estimate $\psi$ through the
method of moments estimator (12), where $\hat{\bm{\beta}}$ is the maximum
quasi-likelihood estimate. To approximately sample from both the posterior
distributions we use a Hamiltonian Monte Carlo (HMC) algorithm provided by the
Stan software (Carpenter et al., 2017). For each simulation, we run 3 chains
of length 3,000 iterations discarding the first 1,000, for a total of 6,000
samples considered for inferential evaluations.

We assess performance by comparing accuracy of the estimated posterior means
of $\bm{\beta}$ and frequentist coverage of credible intervals. The
distributions of the posterior means across the 100 simulations are reported
in Figure 2. Table 2 shows empirical coverage of credible intervals. The
quasi-posterior means outperformed the posterior means under the Gaussian
linear model. Misspecification of the variance function in the normal linear
model analysis leads to a loss of efficiency of the posterior means and poor
frequentist coverage of credible intervals. In contrast, even with a modest
sample size, the quasi-posterior is able to produce well calibrated credible
intervals.

###  5.2 Simulation example: overdispersed count data

It is very common for count data to exhibit over-dispersion relative to what
simple models assume. To show effectiveness of quasi-posteriors in this
context, we provide a simulation study comparing our approach with Poisson and
negative binomial glms, with the latter representing the usual model-based
solution to overdispersion. In particular, we consider the setting of Example
1.

We generate $S=100$ synthetic datasets of sample size $n=1,000$. In each
simulation, the vector of observed responses is generated in two steps: i) for
$i=1,\ldots,n$, we draw $\tilde{y}_{i}$ from
$\mathrm{Ga}(\psi_{0}^{-1}\mu_{i,0},\psi_{0}^{-1})$, where $\psi_{0}=3.5$,
$\mu_{i,0}=\exp\\{x_{i}^{\top}\bm{\beta}_{0}\\}$,
$\bm{\beta}_{0}=(3.5,1.5,-1,0.5)$ and $\bm{x}_{i}$ are sampled from
independent standard normal distributions; ii) the corresponding observed
response $y_{i}$ is obtained rounding $\tilde{y}_{i}$ to the closest integer
value. As shown by Tricker (1984), the effect of rounding on the mean of the
gamma distribution is negligible and on the variance is very mild. Thus, the
$i$th observed response is drawn from a data generator that has mean
$\mu_{i,0}$ and variance approximately equal to $\psi_{0}\mu_{i,0}$. For the
quasi-posterior, the dispersion parameter $\psi$ is estimated using (12) where
in this case $\hat{\bm{\beta}}$ coincides with the maximum likelihood
estimate. Inference on posteriors and quasi-posteriors is carried out with an
HMC algorithm using the Stan software. For each simulation, we run 3 chains of
length 1,000 iterations discarding the first 500, for a total of 1,500 samples
considered for inferential evaluations. Figure 3 compares the distributions of
the posterior mean for $\bm{\beta}$ for the three models across the $S$
simulations, while Table 3 reports coverage of the corresponding $0.95$
credible intervals.

Poisson model and quasi-posterior have similar results in terms of efficiency.
However, coverage differs substantially, with the Poisson model having
frequentist coverage much lower than the nominal level. As discussed in
Example 1, quasi-posterior (9) is equivalent to the power posterior obtained
raising the Poisson likelihood to $\psi^{-1}$. Poisson model and quasi-
posterior express the relation between $\bm{\beta}$ and $\bm{y}$ with the same
structure; consequently it is reasonable to expect that the two posteriors
have similar means. However, the posterior concentration varies substantially.
The misspecified Poisson model produces smaller credible intervals, which tend
to exclude the true parameter value if the posterior mean is not close to
$\bm{\beta}_{0}$. Conversely, the quasi-posterior is less concentrated with
wider credible intervals. As expected, the method of moments estimator is
effective, leading to credible intervals that are well calibrated in their
frequentist coverage. This agrees with the asymptotic results of Section 4.1.

The comparison between negative binomial and quasi-posterior settings offers
additional interesting insights. In terms of coverage calibration, both the
models show good performance since they take into account the dispersion of
the observations, although in two different manners: the first treats the
dispersion as a model parameter and puts a prior on it, while the second uses
its preliminary estimate to regulate the scale of the information coming from
the data within a generalized Bayes approach. However, the negative binomial
glm is misspecified, leading to the loss of efficiency shown in Figure 3. The
quasi-posteriors are substantially more concentrated than the negative
binomial model posteriors, with narrower credible intervals that maintain well
calibrated coverage.

Table 3: Overdispersed count data example. Frequentist coverage of credible intervals with nominal level $\rho\in\\{0.90,0.95,0.99\\}$ for negative binomial model (nb), Poisson model (poi) and quasi-posterior (qp). |  | $\rho=0.90$ |  | $\rho=0.95$ |  | $\rho=0.99$
---|---|---|---|---|---|---
|  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$ |  | $\beta_{1}$ | $\beta_{2}$ | $\beta_{3}$ | $\beta_{4}$
nb |  | 0.79 | 0.84 | 0.92 | 0.94 |  | 0.88 | 0.95 | 0.99 | 0.97 |  | 0.91 | 1.00 | 0.99 | 0.99
poi |  | 0.59 | 0.62 | 0.63 | 0.64 |  | 0.72 | 0.70 | 0.72 | 0.75 |  | 0.89 | 0.80 | 0.86 | 0.85
qp |  | 0.94 | 0.94 | 0.92 | 0.88 |  | 0.96 | 0.97 | 0.97 | 0.95 |  | 1.00 | 1.00 | 0.99 | 0.99
Figure 3: Overdispersed count data example. Distribution of posterior mean of
$\bm{\beta}$ obtained with negative binomial model (nb), Poisson model (poi)
and quasi-posterior (qp).

###  5.3 Application: Rhynchosporium secalis data

Rhynchosporium secalis is a fungus that can infect barley leaves causing leaf
blotches. Wedderburn (1974) reported a study concerning the incidence of
Rhynchosporium secalis on $10$ different varieties of barley grown in $9$
sites, for a total of $n=90$ observations. The response variable is a
continuous proportion in $[0,1]$, representing the percentage of leaf area
infected by the fungus.

Figure 4: Preliminary check of model appropriateness. The first row shows the
comparison between the observed values and the means of the expected value
obtained with mle for each site (top-left) and variety (top-right). The second
row shows the comparison of the observed relation between mean and variance
with the theoretical ones for each site (bottom-left) and for each variety
(bottom-right). For the quasi-posteriors, $\psi$ is estimated through the
method of moments, while for beta regression, $\phi$ is estimated using the
mle. Sites and varieties are sorted in ascending order with respect to the
corresponding mean of the response.

We first considered a full likelihood specification using beta regression,
assuming a beta distribution with logit link function for the proportion of
leaf area infected. Based on fitting such a model, we find that the variance
function is sufficient for small values of the mean but as the mean increases
there is underestimation of variance. In addition, small values of the mean
seem to be overestimated.

This motivated us to consider quasi-likelihood. A first attempt used a quasi-
binomial model with logit link function and
$\mathrm{var}(Y_{i})=\psi\mu_{i}(1-\mu_{i})$, $i=1,\ldots,n$. Results of
preliminary model checking are shown in Figure 4. The quasi-binomial model
exhibits good performance in fitting the mean. However, the choice of variance
function seems poor, as there is clear underestimation of the relationship
between the mean and variance. In attempting to choose a better variance
function, we considered the proposal of Wedderburn (1974) to set
$\mathrm{var}(Y_{i})=\psi\mu_{i}^{2}(1-\mu_{i})^{2}$. Unfortunately, the
variance estimates are similar to the beta model, while the mean estimates are
similar to the quasi-binomial.

An appealing aspect of quasi-likelihood is the great flexibility in choosing
the variance function to better fit data having complex relationships. For
these data, we find that a good choice for the variance is
$\mathrm{var}(Y_{i})=\psi\mu_{i}^{9/4}(1-\mu_{i})^{9/4}$, as shown in Figure
4. In this case both the estimated mean and variance provide a good fit to the
observed data.

Figure 5: Comparison of the posterior expectation for $\mu$ in three
illustrative combinations of site and variety with small, central and large
observed values, respectively.

We conduct detailed comparisons of Bayesian implementations of beta regression
and quasi-posteriors assuming a logistic regression for the mean and variance
functions that are such that either
$\mathrm{var}(Y_{i})=\psi\mu_{i}^{2}(1-\mu_{i})^{2}$ or
$\mathrm{var}(Y_{i})=\psi\mu_{i}^{9/4}(1-\mu_{i})^{9/4}.$ We put a weakly
informative prior on $\bm{\beta}$ and use HMC and random walk Metropolis
Hastings algorithms for sampling from the posterior and the quasi-posteriors,
respectively. The estimates of posterior mean for $\bm{\beta}$ differ across
the approaches. To assess goodness-of-fit, we check whether the 0.95 credible
intervals for each $\mu_{i}$ contain the corresponding observed responses
$y_{i}$. For the quasi-posterior with Wedderburn variance function, the
proportion of intervals containing the observed response is only $0.078.$ This
proportion is $0.456$ for beta regression and $0.656$ for the proposed quasi-
posterior, suggesting significant improvements in fit for the proposed
approach. Comparisons for three combinations of site and variety are
illustrated in Figure 5. Wedderburn quasi-posterior tends to overestimate the
mean, particularly for observations not close to one. The proposed quasi-
posterior works better than beta regression in each variety and in sites with
either small or large mean.

###  5.4 Application: willow warbler abundance in Finland

As an additional application, we analyze data from an ecological monitoring
study of birds focusing on occurrence rates of the willow warbler
(Phylloscopus trochilus) in Finland (Lindström et al., 2015) from 2006 to
2008. The dataset consists of $n=244$ observations collected in $J=171$
different sites. For each observation, the data include a set of environmental
covariates that provide infomation about the habitat type (classified as
either broadleaved forest, coniferous forest, open habitat, urban habitat, or
wetlands) and the mean temperature registered in April-May.

The willow warbler is a long distance migrant passerine and is common in
Finland certain times of the year. After spending the winters in sub-Saharan
Africa, this species arrives in the Fennoscandia (Baltic Shield) region in
late April or beginning of May for the breeding season. For each year,
occurrences were surveyed across all sites between the end of May and
beginning of July, which is after migration ends. Willow warblers arriving in
Finland have substantially different migratory routes than those breeding in
Norway and southern Sweden (Hedenström and Pettersson, 1987). This species
thrives in deciduous forests, but has adapted to living in heavily developed
areas (Cramp and Brooks, 1992). Willow warbler abundance varies substantially
across Finland, and our focus is on assessing to what extent this variation is
explained by temperature and habitat features.

Standard practice in modeling of bird species abundance data focuses on using
Poisson glms that include site-specific random effects in addition to
environmental covariates (Ovaskainen and Abrego, 2020). However, bird
occurrence counts tend to be overdispersed and zero-inflated relative to the
Poisson. For the willow warbler data, the estimated dispersion is
${\hat{\psi}=6.57}$, suggesting substantial overdispersion relative to the
Poisson. Stoklosa et al. (2022) advocates for the general use of negative
binomial models for species abundance in ecology. Hence, we focus on a
negative binomial hierarchical regression model as a competitor to our quasi-
posterior approach. As zero counts are rare in these data ($0.4\%$ of the
observations), zero-inflated and hurdle models do not seem appropriate. We
assume the following mean and variance functions in constructing our quasi-
posterior:

$\mathds{E}(Y_{ij})=\mu_{ij},\quad\mathrm{var}(Y_{ij})=\psi\mu_{ij},\quad\mu_{ij}=\exp\\{\bm{x}_{ij}^{\top}\bm{\beta}+\delta_{j}\\},\quad\delta_{j}\sim\mathcal{N}(0,\sigma),$

for $i=1,\ldots,n_{j}$, $j=1,\ldots,J$, $\sum_{j=1}^{J}n_{j}=n$. We focus on
the log link, as it leads to easy interpretation of regression coefficients
and is standard in ecological abundance modeling.

The covariate vector $\bm{x}_{ij}$ includes temperature, habitat, and year.
Specifically, the temperature variable corresponds to the average temperature
(in degrees Celsius) recorded during April and May in each site for the
corresponding year. The average temperature during these months may impact
bird preferences for breeding sites and movement patterns. The different
habitats mentioned above are included as dummy indicators, with broadleaved
forest considered as the baseline as the primary native habitat of this
species. Average spring temperature and habitat for each location are reported
in Appendix B. Finally, the survey year is included to account for possible
annual effects. This is encoded in the covariates into two dummy variables,
taking 2006 as baseline.

Table 4: Posterior mean and 0.95 credible intervals for $\bm{\beta}$ components and $\sigma$. Parameters in italics refer to habitat categories (broadleaved forest is taken as baseline), while 2006 is the baseline for year components. |  | poi |  | nb |  | qp
---|---|---|---|---|---|---
Intercept |  | $4.10$ |  [$3.80$, $4.38$] |  | $3.71$ |  [$3.43$, $3.98$] |  | $3.74$ |  [$3.47$, $3.99$]
Temperature |  | $-0.07$ |  [$-0.10$, $-0.03$] |  | $0.03$ |  [$-0.01$, $0.07$] |  | $0.02$ |  [$-0.01$, $0.05$]
Year |  |  |  |  |  |  |  |  |
2007 |  | $0.28$ |  [$0.21$, $0.35$] |  | $0.30$ |  [$0.11$, $0.49$] |  | $0.27$ |  [$0.11$, $0.41$]
2008 |  | $-0.03$ |  [$-0.10$, $0.05$] |  | $0.03$ |  [$-0.14$, $0.20$] |  | $0.00$ |  [$-0.15$, $0.14$]
Habitat |  |  |  |  |  |  |  |  |
Coniferous f. |  | $-0.26$ |  [$-0.50$, $-0.01$] |  | $-0.25$ |  [$-0.47$, $-0.03$] |  | $-0.24$ |  [$-0.45$, $-0.06$]
Open habitat |  | $-0.74$ |  [$-1.07$, $-0.43$] |  | $-0.44$ |  [$-0.72$, $-0.13$] |  | $-0.42$ |  [$-0.69$, $-0.17$]
Urban habitat |  | $-0.54$ |  [$-0.88$, $-0.19$] |  | $-0.77$ |  [$-1.09$, $-0.46$] |  | $-0.73$ |  [$-1.03$, $-0.42$]
Wetlands |  | $-0.37$ |  [$-0.79$, $0.03$] |  | $-0.16$ |  [$-0.54$, $0.22$] |  | $-0.18$ |  [$-0.52$, $0.17$]
Site |  |  |  |  |  |  |  |  |
$\sigma$ |  | $0.34$ |  [$0.25,0.45]$ |  | $0.07$ | $(0.00,0.14]$ |  | $0.34$ | $[0.27,0.42]$

The inclusion of random intercepts
$\bm{\delta}=(\delta_{1},\ldots,\delta_{J})$ accounts for the hierarchical
structure of the data and heterogeneity among the locations in willow warbler
occurrence rates; for example, this heterogeneity may be due to differences in
habitat that are not captured by measured covariates. The common Gaussian
prior centered at zero shrinks the site-specific random intercepts toward the
population mean.

We compare our quasi-posterior with Poisson and negative binomial regression
models having an identical mean structure to ours, including site-specific
random effects. We put weakly informative priors on $\bm{\beta}$ and $\sigma$.
For each model, we sample from the posterior distribution using an HMC
algorithm in Stan running 4 chains of length 2,000 iterations discarding the
first 500, for a total of 6,000 samples considered for inferential
evaluations. To assess goodness-of-fit, we compare the in-sample standardized
mean squared error (sMSE) of the three models computed using Pearson
residuals, which remove the effect of assumed heteroscedasticity (Agresti,
2015). As expected, the Poisson model has the highest sMSE ($1.72$) since it
underestimates the dispersion, while the quasi-posterior and negative binomial
model have substantially better performances since both account for
overdispersion, with a better fit for the quasi-posterior ($0.59$ sMSE vs
$0.65$ for the negative binomial). For this reason, we focus on the
inferential conclusions of these two.

Figure 6: Panel (a) represents the quasi-posterior means of site-specific
random effects and their geographical location in Finland. Site 3 is
highlighted in blue and its posteriors distributions are represented in panel
(b). Site 59 is highlighted in red and its posteriors distributions are
represented in panel (c).

Posterior means and 0.95 credible intervals for $\bm{\beta}$ are reported in
Table 4. As expected, the negative binomial model and quasi-posterior provide
similar estimates for the regression coefficients, and their posterior
distributions have higher variance than those under the Poisson model. The
results suggest that 2007 was favorable for willow warbler abundance in
Finland compared to 2006 and 2008. Causes of the differences in abundance
across years may be due to unmeasured weather covariates. Willow warbler
abundance is lower in coniferous forests, open and urban habitats compared to
broadleaved forests. These conclusions are in line with literature noting that
willow warblers prefer to feed in deciduous forests (Virkkala, 1991) and that
in urban areas there is a reduction of the number of Phylloscopus trochilus
males due to traffic (Reijnen and Foppen, 1991, 1994). Moreover, Phylloscopus
trochilus diffusion does not significantly vary according to the average
spring temperature across the locations. This indicates that the temperature
the birds encounter in Finland after their spring migration does not affect
their spreading over the whole territory, which is influenced mainly by the
habitat.

The posterior distribution of the random intercepts helps identify differences
in willow warbler abundance across locations that are not explained by our
observed habitat and temperature covariates. In interpreting these results, we
focus on our quasi-posterior for brevity. As shown in Figure 6a, sites having
unusually low willow warbler abundance are mainly located in southern Finland
and upper Lapland, while sites with unusually high abundance concentrate along
the Bothnian Bay coast and in lower Lapland. It is likely that these spatial
patterns reflect additional more nuanced habitat features not captured in the
coarse habitat classification provided in the data for each site. In future
work it will be interesting to look more into these patterns to infer these
latent habitat factors to learn more about what drives bird population
dynamics, and expand habitat classification in future studies.

Probing the posterior on the random effects allows identification of
interesting locations. For example, Figure 6b and 6c show two sites with
unexpected high abundances - one in urban Helsinki and another in a central
Finland coniferous forest. The former may represent a success story in terms
of providing a good habitat within a city environment. This success story may
inform future policies for protecting green areas in cities. The later would
be interesting to follow-up to assess whether there is something special about
that coniferous forest causing high abundance of a deciduous forest species.

##  6 Discussion

Our proposed quasi-posteriors provide a generalized Bayes solution to overcome
misspecification problems in Bayesian glms. Quasi-posteriors additionally have
interesting links to c-posteriors, giving new insights into coarsening. We
have demonstrated that our quasi-posteriors are highly versatile, allowing one
to focus on mean and variance functions in model specification, simplifying
the search for an appropriate model for data with challenging
heteroscedasticity. An interesting application is to count response data,
providing notable improvements over not only Poisson regression but also
negative binomial-based glms. We also see substantial improvements for
continuous response data having variance changing with covariates. Addition of
random effects is straightforward, making our method convenient for many
different applications. The code to replicate our analyses is available online
in the repository https://github.com/davideagno/Quasi-posteriors.

Quasi-posteriors are designed for inference on regression coefficients. For
simplicity and ease in general implementation, we use a method of moments
estimator for the dispersion parameter $\psi$, following standard practice in
the frequentist literature. An interesting area of future research is to
develop approaches for uncertainty quantification for $\psi$. Potentially a
prior can be chosen for $\psi$ and a Bayesian approach can be developed for
inference, building on methods in Bissiri et al. (2016). However, there are
some associated subtleties that would need to be addressed.

An additional interesting area for future research is the development of
formal model selection procedures. We illustrated one heuristic strategy in
our analysis in Section 5.3; in particular, in this case we considered a
variety of different variance functions and selected the one with the best
goodness-of-fit to the data. As alternative model performance criteria for use
in model comparisons, one could rely on leave-one-out cross-validation and
marginal quasi-likelihood as in Fong and Holmes (2020), or adaptations of
deviance information criteria and related ideas (Gelman et al., 2013). This
can be applied to selection of both the mean and variance functions.
Interesting next steps include theoretical support for such model selection
strategies, including frequentist asymptotic theory and a Bayesian decision
theoretic justification.

## Funding

This research was partially supported by the Lifeplan project funded by the
European Research Council under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No 856506). Additional funding came from
the United States National Institutes of Health R01-ES028804 and R01-ES035625.

##  Appendix A Proofs

###  A.1 Proof of Theorem 1

The inferential target for Gibbs posterior is

$\bm{\beta}_{\textsc{opt}}=\arg\max_{\bm{\beta}}\int\ell_{\textsc{q}}(\bm{\beta},\bm{y},\bm{X},\psi)F_{0}(\mathrm{d}\bm{y}\mid\bm{X})=\arg\max_{\bm{\beta}}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi),$

where $F_{0}(\mathrm{d}\bm{y}\mid\bm{X})$ denotes the probability measure
corresponding to the density
$f_{0}(\bm{y}\mid\bm{X})=f(\bm{y}\mid\bm{X},\bm{\beta}_{0},\psi_{0})=\prod_{i=1}^{n}f(y_{i}\mid\bm{x}_{i},\bm{\beta}_{0},\psi_{0})$,
since $Y_{i}\mid\bm{x}_{i}$ are independent by assumption, and having defined

$\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)=\sum_{i=1}^{n}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{x}_{i},\psi)=\sum_{i=1}^{n}\int_{\mathcal{Y}}\int_{a}^{\mu_{i}}\frac{y_{i}-t}{\psi
V(t)}\mathrm{d}tF_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i}),$

where
$\mu_{i}=\mu_{i}(\bm{\beta})=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta})\in\mathcal{M}\subseteq\mathds{R}$.
By Fubini-Tonelli theorem, the integral

$\int_{\mathcal{Y}\times[a,\mu_{i}]}\left|\frac{y_{i}-t}{\psi
V(t)}\right|\mathrm{d}tF_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})=\int_{\mathcal{Y}}\int_{a}^{\mu_{i}}\left|\frac{y_{i}-t}{\psi
V(t)}\right|\mathrm{d}tF_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})=\int_{a}^{\mu_{i}}\int_{\mathcal{Y}}\left|\frac{y_{i}-t}{\psi
V(t)}\right|F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})\mathrm{d}t$

is well defined for any $i=1,\ldots,n$, albeit potentially infinite. By
assumption, the variance function $\psi V(t)$ is strictly positive and
$\int|y_{i}|F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})\leq\left(\int
y_{i}^{2}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})\right)^{1/2}=c_{i}<\infty$
because the second moment is finite for all $i=1,\ldots,n$. Then, as an
application of the triangular inequality, we obtain the following chain of
inequalities

$\int_{a}^{\mu_{i}}\int_{\mathcal{Y}}\left|\frac{y_{i}-t}{\psi
V(t)}\right|F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})\mathrm{d}t\leq\int_{a}^{\mu_{i}}\frac{1}{\psi
V(t)}\left(\int_{\mathcal{Y}}|y_{i}|F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})+|t|\right)\mathrm{d}t\leq\int_{a}^{\mu_{i}}\frac{c_{i}-|t|}{\psi
V(t)}\mathrm{d}t<\infty,$

where the last inequality follows because the function
$(c_{i}-|t|)\psi^{-1}V(t)^{-1}$ is continuous and bounded for any
$t\in[a,\mu_{i}]$ and $i=1,\ldots,n$. Thus, the function
$(y_{i}-t)\psi^{-1}V(t)^{-1}$ is integrable in the product space
$\mathcal{Y}\times[a,\mu_{i}]$ with respect to the product measure and the
order of integration can be interchanged, giving

$\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)=\sum_{i=1}^{n}\int_{a}^{\mu_{i}}\int_{\mathcal{Y}}\frac{y_{i}-t}{\psi
V(t)}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})\mathrm{d}t=\sum_{i=1}^{n}\int_{a}^{\mu_{i}}\frac{g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})-t}{\psi
V(t)}\mathrm{d}t,$

where
$g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})=\int_{\mathcal{Y}}y_{i}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})$.
The gradient vector

$\nabla\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)=\frac{\partial}{\partial\bm{\beta}}\sum_{i=1}^{n}\int_{a}^{\mu_{i}}\frac{g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})-t}{\psi
V(t)}\mathrm{d}t=\sum_{i=1}^{n}\frac{g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})-g^{-1}(\bm{x}_{i}^{\top}\bm{\beta})}{\psi
V(\mu_{i})g^{\prime}(\mu_{i})}\bm{x}_{i}$

is zero if $\bm{\beta}=\bm{\beta}_{0}$. The vector $\bm{\beta}_{0}$ is a
maximum because the second derivative
$\nabla^{2}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)$, evaluated in
$\bm{\beta}_{0}$, is negative definite. In fact, the expected information
defined in (10) evaluated in $\bm{\beta}_{0}$ equals to

$I(\bm{\beta}_{0})=-\nabla^{2}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)|_{\bm{\beta}=\bm{\beta}_{0}}=\psi^{-1}\bm{X}^{\top}\tilde{\bm{D}}\bm{X},$

where $\tilde{\bm{D}}=\text{diag}(\tilde{d}_{1},\ldots,\tilde{d}_{n})$,
$\tilde{d}_{i}=\big{[}V(\mu_{i}(\bm{\beta}_{0}))g^{\prime}(\mu_{i}(\bm{\beta}_{0}))^{2}\big{]}^{-1}$
for $i=1,\ldots,n$, and is positive definite because the design matrix
$\bm{X}$ has full rank for any $n\geq p$ and the entries $\tilde{d}_{i}>0$ are
strictly positive by assumption.

The value $\bm{\beta}_{0}$ is the only solution in many concrete scenarios.
For instance, suppose the identity $g^{\prime}(\mu)=1/V(\mu)$ holds,
corresponding to a quasi-likelihood with canonical link such as overdispersed
Poisson regression with log link and logistic regression; then the integrated
loss $\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)$ is concave, because
$I(\bm{\beta})=-\nabla^{2}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)$
for all values of $\bm{\beta}$. In this case, $\bm{\beta}_{0}$ is the unique
solution.

###  A.2 Proof of Proposition 1

The power posterior approximation of the corresponding c-posterior based on
(8) is

$p_{\alpha}(\bm{\beta}\mid\bm{y},\bm{X})\propto
p(\bm{\beta})\prod_{i=1}^{n}\exp\left\\{\frac{\alpha}{\alpha+n}\Big{[}\theta(\mu_{i})y_{i}-b(\theta(\mu_{i}))\Big{]}\right\\}.$
(13)

By assumption, the second order conditions (1) are such that
$\mathrm{var}(Y_{i})=\psi V(\mu_{i})$, $i=1,\ldots,n$, with

$V(\mu_{i})=\frac{\partial^{2}}{\partial\theta^{2}}b(\theta)\bigg{|}_{\theta=\theta(\mu_{i})},$

and, by Theorem 2, the corresponding quasi-posterior is

$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)\propto
p(\bm{\beta})\prod_{i=1}^{n}\exp\left\\{\frac{1}{\psi}\Big{[}\theta(\mu_{i})y_{i}-b(\theta(\mu_{i}))\Big{]}\right\\},$

that is equivalent to (13) for $\psi=(\alpha+n)/\alpha$.

###  A.3 Proof of Theorem 3 and regularity conditions

The result follows from the general theory developed in Theorems 4 and 5 of
Miller (2021), in which are presented abstract sufficient conditions for which
a generalized posterior converges in total variation to a Gaussian
distribution. We review these conditions and show how they can be specialized
in the context of quasi posteriors. Let $B\subseteq\mathds{R}^{p}$ be an open
and convex set such that $\bm{\beta}_{0}\in B$. In the first place, we need to
ensure that the function $\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$ is
well-behaved, namely that $\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$
has continuous third derivatives on $B$ and that the third derivative is
uniformly bounded on $B$; these are very mild conditions that are satisfied
for any reasonable choice of the link $g(\cdot)$ and variance $V(\cdot)$
functions. In the second place, it is required in Miller (2021, Theorem 5)
that the following quantity

$\frac{1}{n}\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)=\frac{1}{n}\sum_{i=1}^{n}\ell_{\textsc{q}}(\bm{\beta};y_{i},\bm{x}_{i},\psi)=\frac{1}{n}\sum_{i=1}^{n}\frac{1}{\psi}\int_{a}^{\mu_{i}(\bm{\beta})}\frac{y_{i}-t}{V(t)}\mathrm{d}t,$

converges pointwise and almost surely to a well defined limit, say
$\ell_{\infty}(\bm{\beta})$, for any $\bm{\beta}\in B$, as
$n\rightarrow\infty$. In quasi-posteriors, the contributions
$\ell_{\textsc{q}}(\bm{\beta};y_{1},\bm{x}_{1},\psi),\dots,\ell_{\textsc{q}}(\bm{\beta};y_{n},\bm{x}_{n},\psi)$
are independent realizations, therefore extensions of the law of large numbers
can be invoked. Note that each
$\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{x}_{i},\psi)=\mathds{E}(\ell_{\textsc{q}}(\bm{\beta};y_{i},\bm{x}_{i},\psi))$
is finite, as already discussed in Section A.1. Hence, sufficient conditions
for pointwise convergence to $\ell_{\infty}(\bm{\beta})$ are

$\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{x}_{i},\psi)=\ell_{\infty}(\bm{\beta})<\infty,\qquad\sum_{n=1}^{\infty}\frac{\text{var}(\ell_{\textsc{q}}(\bm{\beta};Y_{n},\bm{x}_{n},\psi))}{n^{2}}<\infty.$

as a consequence of Proposition 6.13, Cinlar (2011, Chapter 3). Both these
conditions are very mild and essentially requires that the quasi-log-
likelihood is well-defined at the limit. Thirdly, regularity conditions on the
limiting function $\ell_{\infty}(\bm{\beta})$ are needed. In particular, it is
required that $i(\bm{\beta}_{0})=-\nabla^{2}\ell_{\infty}(\bm{\beta}_{0})$ is
a positive definite matrix. As shown in the proof of Appendix A.1,
$I(\bm{\beta}_{0})=-\nabla^{2}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)|_{\bm{\beta}=\bm{\beta}_{0}}=\psi^{-1}\bm{X}^{\top}\tilde{\bm{D}}\bm{X}$
is positive definite. Hence, recalling that
$\ell_{\infty}(\bm{\beta})=\lim_{n\rightarrow\infty}n^{-1}\sum_{i=1}^{n}\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{x}_{i},\psi)$,
we assume that
$i(\bm{\beta}_{0})=\lim_{n\rightarrow\infty}n^{-1}I(\bm{\beta}_{0})=\lim_{n\rightarrow\infty}(n\psi)^{-1}\bm{X}^{\top}\tilde{\bm{D}}\bm{X}$
is positive definite. Finally, either of the following two conditions must be
satisfied:

  1. i)

$\ell_{\infty}(\bm{\beta})<\ell_{\infty}(\bm{\beta}_{0})$ for all
$\bm{\beta}\in B_{*}\setminus\\{\bm{\beta}_{0}\\}$ and
$\lim\sup_{n}\sup_{\bm{\beta}\in\mathds{R}^{p}\setminus
B_{*}}\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)<\ell_{\infty}(\bm{\beta}_{0})$
for some compact set $B_{*}\subseteq B$ with $\bm{\beta}_{0}$ in the interior
of $B_{*}$;

  2. ii)

$\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$ is concave and
$\nabla\ell_{\infty}(\bm{\beta}_{0})=0$.

The requirement (i) is a Wald-type condition and essentially it means that the
function $\ell_{\infty}(\bm{\beta})$ is maximized in $\bm{\beta}_{0}$.
Condition (ii) is stronger and it implies (i), but is satisfied in many
practical occasions. Indeed, paralleling the reasoning of the proof of
Appendix A.1, if the identity $g^{\prime}(\mu)=1/V(\mu)$ holds, then the loss
$\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$ is concave, because we get
$I(\bm{\beta})=-\nabla^{2}\ell_{\textsc{q}}(\bm{\beta};\bm{y},\bm{X},\psi)$,
which is positive definite for all values of $\bm{\beta}$. Moreover, note that
$\nabla\ell_{\infty}(\bm{\beta}_{0})=\lim_{n\rightarrow\infty}n^{-1}\nabla\mathcal{L}(\bm{\beta};\bm{\beta}_{0},\bm{X},\psi)|_{\bm{\beta}=\bm{\beta}_{0}}=0$,
as shown in Appendix A.1.

###  A.4 Proof of Theorem 4

The result follows from Theorem 8 in Miller (2021), that provides the
sufficient conditions for which a generalized posterior exhibits correct
frequentist coverage. In the first place, we require the convergence in
distribution of $-\sqrt{n}(\hat{\bm{\beta}}-\bm{\beta}_{0})$ to a random
variable $Z$ as $n\rightarrow\infty$. In our case this is true because, under
the regularity assumptions of Theorem 3, we have
$\sqrt{n}(\hat{\bm{\beta}}-\bm{\beta}_{0})\xrightarrow{d}Z$, with $Z\sim
N(0,i(\bm{\beta}_{0})^{-1})$; see e.g. McCullagh (1983). The second condition
requires that the distribution $\pi_{n}$ of
$\sqrt{n}(\bm{\beta}-\hat{\bm{\beta}})$ converges in total variation to the
law of $Z$, which is guaranteed by Theorem 3 under the regularity conditions
presented in Appendix A.3. The third condition in Miller (2021) requires the
pointwise convergence of the sequence of centered and rescaled credible sets
$C_{n}=\\{\sqrt{n}(\bm{\beta}-\hat{\bm{\beta}})\in\mathds{R}^{p}:\bm{\beta}\in
S_{n}\\}$ to a fixed set $C\subseteq\mathds{R}^{p}$ with finite nonzero
Lebesgue measure, that is, for all $\bm{u}\in\mathds{R}^{p}$,
$\mathds{1}(\bm{u}\in C_{n})\rightarrow\mathds{1}(\bm{u}\in C)$ almost surely
as $n\rightarrow\infty$. The validity of this requirement depends on the
specific choice of the credible sets. As an example, fix
$C=\\{\bm{u}\in\mathds{R}^{p}:\mathcal{N}(\bm{u}\mid
0,i(\bm{\beta}_{0})^{-1})\geq k_{\rho}\\}$, where $\mathcal{N}(\cdot\mid
0,i(\bm{\beta}_{0})^{-1})$ denotes the density function of $Z$, and $k_{\rho}$
is such that $\int_{C}\mathcal{N}(\bm{u}\mid
0,i(\bm{\beta}_{0})^{-1})\mathrm{d}\bm{u}=\rho$. Suppose
$p_{\textsc{q}}(\bm{\beta}\mid\bm{y},\bm{X},\psi)$ is unimodal and $S_{n}$ are
highest posterior density credible sets. Then
$C_{n}=\\{\bm{u}\in\mathds{R}^{p}:p_{n}(\bm{u})\geq k_{\rho}\\}$, where
$p_{n}(\cdot)$ is the density function of
$\sqrt{n}(\bm{\beta}-\hat{\bm{\beta}})$. Then,

$\displaystyle\mathds{1}(\bm{u}\in C_{n})=1\quad$ $\displaystyle\text{if and
only if}\quad p_{n}(\bm{u})\geq k_{\rho},$ $\displaystyle\mathds{1}(\bm{u}\in
C)=1\quad$ $\displaystyle\text{if and only if}\quad\mathcal{N}(\bm{u}\mid
0,i(\bm{\beta}_{0})^{-1})\geq k_{\rho}.$

Therefore, for all $\bm{u}\in\mathds{R}^{p}$, $\mathds{1}(\bm{u}\in
C_{n})\xrightarrow{a.s.}\mathds{1}(\bm{u}\in C)$ as $n\rightarrow\infty$,
because

$\int_{\mathds{R}^{p}}\left|p_{n}(\bm{u})-\mathcal{N}(\bm{u}\mid
0,i(\bm{\beta}_{0})^{-1})\right|\mathrm{d}\bm{u}\rightarrow 0$

almost surely under $f_{0}$ by our Theorem 3. The final condition in Miller
(2021) is automatically satisfied as $Z$ is absolutely continuous, being a
multivariate Gaussian.

###  A.5 Proof of Proposition 2

Recalling (7), we can write
$\ell_{\textsc{q}}(\bm{\beta},\bm{y},\bm{X},\psi)=\sum_{i=1}^{n}\ell_{\textsc{q}}(\bm{\beta},y_{i},\bm{x}_{i},\psi)=-\psi^{-1}\sum_{i=1}^{n}\ell(\bm{\beta};y_{i},\bm{x}_{i})$
with

$\ell(\bm{\beta};y_{i},\bm{x}_{i})=-\int_{a}^{\mu_{i}}\frac{y_{i}-t}{V(t)}\mathrm{d}t,$

where
$\mu_{i}=\mu_{i}(\bm{\beta})=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta})\in\mathcal{M}\subseteq\mathds{R}$.
Then

$\displaystyle J(\bm{\beta})$
$\displaystyle=\int\nabla^{2}\ell(\bm{\beta};\bm{y},\bm{X})F_{0}(\mathrm{d}\bm{y}\mid\bm{X})=\sum_{i=1}^{n}\int\nabla^{2}\ell(\bm{\beta};y_{i},\bm{x}_{i})F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i}),$
$\displaystyle H(\bm{\beta})$
$\displaystyle=\int\nabla\ell(\bm{\beta};\bm{y},\bm{X})\nabla\ell(\bm{\beta};\bm{y},\bm{X})^{\top}F_{0}(\mathrm{d}\bm{y}\mid\bm{X})=\sum_{i=1}^{n}\int\nabla\ell(\bm{\beta};y_{i},\bm{x}_{i})\nabla\ell(\bm{\beta};y_{i},\bm{x}_{i})^{\top}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i}).$

Recall that $\int
y_{i}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})=g^{-1}(\bm{x}_{i}^{\top}\bm{\beta}_{0})=\mu_{i}(\bm{\beta}_{0})$
and
$\int[y_{i}-\mu_{i}(\bm{\beta}_{0})]^{2}F_{0}(\mathrm{d}y_{i}\mid\bm{x}_{i})=\psi_{0}V(\mu_{i}(\bm{\beta}_{0}))$.
Hence, after some algebra (Wedderburn, 1974), we obtain

$J(\bm{\beta}_{0})=\bm{X}^{\top}\tilde{\bm{D}}\bm{X},\qquad
H(\bm{\beta}_{0})=\psi_{0}\bm{X}^{\top}\tilde{\bm{D}}\bm{X}=\psi_{0}J(\bm{\beta}_{0}).$

Consequently,

$\displaystyle j(\bm{\beta}_{0})$
$\displaystyle=\lim_{n\rightarrow\infty}n^{-1}\mathds{E}\left[\nabla^{2}\ell(\bm{\beta}_{0};\bm{Y},\bm{X})\right]=\lim_{n\rightarrow\infty}n^{-1}J(\bm{\beta}_{0}),$
$\displaystyle h(\bm{\beta}_{0})$
$\displaystyle=\lim_{n\rightarrow\infty}n^{-1}\mathds{E}\left[\nabla\ell(\bm{\beta}_{0};\bm{Y},\bm{X})\nabla\ell(\bm{\beta}_{0};\bm{Y},\bm{X})^{\top}\right]=\lim_{n\rightarrow\infty}n^{-1}H(\bm{\beta}_{0})=\psi_{0}j(\bm{\beta}_{0}).$

Therefore, (11) can be simplified as

$\psi_{\textsc{llb}}=\frac{\text{tr}\\{j(\bm{\beta}_{0})\\}}{\text{tr}\\{j(\bm{\beta}_{0})h(\bm{\beta}_{0})^{-1}j(\bm{\beta}_{0})\\}}=\psi_{0}.$

##  Appendix B Summary of environmental covariates of willow warblers
abundance analysis

Figure 7 contains the summary of the environmental covariates used in Section
5.4

Figure 7: Habitat type (left) and average temperature in April-May across
2006-2008 (right) for each site location.

## References

  * Agresti (2015) Agresti, A. (2015).  Foundations of linear and generalized linear models.  John Wiley & Sons.
  * Bissiri et al. (2016) Bissiri, P. G., C. C. Holmes, and S. G. Walker (2016).  A general framework for updating belief distributions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 78(5), 1103–1130.
  * Carpenter et al. (2017) Carpenter, B., A. Gelman, M. D. Hoffman, D. Lee, B. Goodrich, M. Betancourt, M. A. Brubaker, J. Guo, P. Li, and A. Riddell (2017).  Stan: A probabilistic programming language.  Journal of statistical software 76(1).
  * Chernozhukov and Hong (2003) Chernozhukov, V. and H. Hong (2003).  An MCMC approach to classical estimation.  Journal of econometrics 115(2), 293–346.
  * Cinlar (2011) Cinlar, E. (2011).  Probability and Stochastics.  Springer.
  * Cramp and Brooks (1992) Cramp, S. and D. Brooks (1992).  Handbook of the birds of Europe, the Middle East and North Africa. The birds of the western Palearctic, vol. VI. Warblers, Volume 906.  Oxford: Oxford university press.
  * Fong and Holmes (2020) Fong, E. and C. C. Holmes (2020).  On the marginal likelihood and cross-validation.  Biometrika 107(2), 489–496.
  * Gelman et al. (2013) Gelman, A., J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin (2013).  Bayesian Data Analysis (3rd ed.).  Chapman and Hall/CRC.
  * Greco et al. (2008) Greco, L., W. Racugno, and L. Ventura (2008).  Robust likelihood functions in Bayesian inference.  Journal of Statistical Planning and Inference 138(5), 1258–1270.
  * Grünwald and Van Ommen (2017) Grünwald, P. and T. Van Ommen (2017).  Inconsistency of Bayesian inference for misspecified linear models, and a proposal for repairing it.  Bayesian Analysis 12(4), 1069–1103.
  * Hedenström and Pettersson (1987) Hedenström, A. and J. Pettersson (1987).  Migration routes and wintering areas of willow warblers phylloscopus trochilus (l.) ringed in fennoscandia.  Ornis Fennica 64(4), 137–143.
  * Hoff (2007) Hoff, P. D. (2007).  Extending the rank likelihood for semiparametric copula estimation.  Annals of Applied Statistics 1(1), 265–283.
  * Holmes and Walker (2017) Holmes, C. C. and S. G. Walker (2017).  Assigning a value to a power likelihood in a general Bayesian model.  Biometrika 104(2), 497–503.
  * Jewson and Rossell (2022) Jewson, J. and D. Rossell (2022).  General bayesian loss function selection and the use of improper models.  Journal of the Royal Statistical Society Series B: Statistical Methodology 84(5), 1640–1665.
  * Jiang and Tanner (2008) Jiang, W. and M. A. Tanner (2008).  Gibbs posterior for variable selection in high-dimensional classification and data mining.  Annals of Statistics 36(5), 2207–2231.
  * Kim and Kim (2009) Kim, Y. and D. Kim (2009).  Bayesian partial likelihood approach for tied observations.  Journal of statistical planning and inference 139(2), 469–477.
  * Lazar (2003) Lazar, N. A. (2003).  Bayesian empirical likelihood.  Biometrika 90(2), 319–326.
  * Lin (2006) Lin, L. (2006).  Quasi Bayesian likelihood.  Statistical Methodology 3(4), 444–455.
  * Lindström et al. (2015) Lindström, Å., M. Green, M. Husby, J. A. Kålås, and A. Lehikoinen (2015).  Large-scale monitoring of waders on their boreal and arctic breeding grounds in northern europe.  Ardea 103(1), 3–15.
  * Lyddon et al. (2019) Lyddon, S. P., C. C. Holmes, and S. G. Walker (2019).  General Bayesian updating and the loss-likelihood bootstrap.  Biometrika 106(2), 465–478.
  * McCullagh (1983) McCullagh, P. (1983).  Quasi-likelihood functions.  Annals of Statistics 11(1), 59–67.
  * McCullagh and Nelder (1989) McCullagh, P. and J. A. Nelder (1989).  Generalized Linear Models.  Chapman & Hall.
  * Miller (2021) Miller, J. W. (2021).  Asymptotic normality, concentration, and coverage of generalized posteriors.  The Journal of Machine Learning Research 22(1), 7598–7650.
  * Miller and Dunson (2018) Miller, J. W. and D. B. Dunson (2018).  Robust Bayesian inference via coarsening.  Journal of the American Statistical Association 114(527), 1113–1125.
  * Müller (2013) Müller, U. K. (2013).  Risk of Bayesian inference in misspecified models, and the sandwich covariance matrix.  Econometrica 81(5), 1805–1849.
  * Nelder and Wedderburn (1972) Nelder, J. A. and R. W. Wedderburn (1972).  Generalized linear models.  Journal of the Royal Statistical Society: Series A (General) 135(3), 370–384.
  * Ovaskainen and Abrego (2020) Ovaskainen, O. and N. Abrego (2020).  Joint species distribution modelling: With applications in R.  Cambridge University Press.
  * Pauli et al. (2011) Pauli, F., W. Racugno, and L. Ventura (2011).  Bayesian composite marginal likelihoods.  Statistica Sinica, 149–164.
  * Raftery et al. (1996) Raftery, A. E., D. Madigan, and C. T. Volinsky (1996).  Accounting for model uncertainty in survival analysis improves predictive performance.  Bayesian Statistics 5, 323–349.
  * Reijnen and Foppen (1991) Reijnen, R. and R. Foppen (1991).  Effect of road traffic on the breeding site-tenacity of male willow warblers (phylloscopus trochilus).  Journal für Ornithologie 132, 291–295.
  * Reijnen and Foppen (1994) Reijnen, R. and R. Foppen (1994).  The effects of car traffic on breeding bird populations in woodland. i. evidence of reduced habitat quality for willow warblers (phylloscopus trochilus) breeding close to a highway.  Journal of Applied ecology, 85–94.
  * Ribatet et al. (2012) Ribatet, M., D. Cooley, and A. C. Davison (2012).  Bayesian inference from composite likelihoods, with an application to spatial extremes.  Statistica Sinica, 813–845.
  * Rigon et al. (2023) Rigon, T., A. H. Herring, and D. B. Dunson (2023).  A generalized Bayes framework for probabilistic clustering.  Biometrika 10(3), 559–578.
  * Royall and Tsou (2003) Royall, R. and T.-S. Tsou (2003).  Interpreting statistical evidence by using imperfect models: robust adjusted likelihood functions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 65(2), 391–404.
  * Rubin (1981) Rubin, D. B. (1981).  The Bayesian bootstrap.  Annals of Statistics, 130–134.
  * Sinha et al. (2003) Sinha, D., J. G. Ibrahim, and M.-H. Chen (2003).  A Bayesian justification of Cox’s partial likelihood.  Biometrika 90(3), 629–641.
  * Smith and Stephenson (2009) Smith, E. L. and A. G. Stephenson (2009).  An extended Gaussian max-stable process model for spatial extremes.  Journal of Statistical Planning and Inference 139(4), 1266–1275.
  * Stoklosa et al. (2022) Stoklosa, J., R. V. Blakey, and F. K. Hui (2022).  An overview of modern applications of negative binomial modelling in ecology and biodiversity.  Diversity 14(5), 320.
  * Tricker (1984) Tricker, A. (1984).  Effects of rounding on the moments of a probability distribution.  Journal of the Royal Statistical Society Series D: The Statistician 33(4), 381–390.
  * Ventura et al. (2010) Ventura, L., S. Cabras, and W. Racugno (2010).  Default prior distributions from quasi- and quasi-profile likelihoods.  Journal of Statistical Planning and Inference 140(11), 2937–2942.
  * Ventura and Racugno (2016) Ventura, L. and W. Racugno (2016).  Pseudo-likelihoods for Bayesian inference.  In Topics on methodological and applied statistical inference, pp. 205–220. Springer.
  * Virkkala (1991) Virkkala, R. (1991).  Spatial and temporal variation in bird communities and populations in north-boreal coniferous forests: a multiscale approach.  Oikos, 59–66.
  * Wedderburn (1974) Wedderburn, R. W. (1974).  Quasi-likelihood functions, generalized linear models, and the Gauss-Newton method.  Biometrika 61(3), 439–447.
  * Yin (2009) Yin, G. (2009).  Bayesian generalized methood of moments.  Bayesian Analysis 4(2), 119–208.
  * Zhang (2006) Zhang, T. (2006).  Information-theoretic upper and lower bounds for statistical estimation.  IEEE Transactions on Information Theory 52(4), 1307–1321.