# Bayesian Nonparametric Quasi Likelihood

Antonio R. Linero Department of Statistics and Data Sciences, University of
Texas at Austin<EMAIL_ADDRESS>
###### Abstract

A recent trend in Bayesian research has been revisiting generalizations of the
likelihood that enable Bayesian inference without requiring the specification
of a model for the data generating mechanism. This paper focuses on a Bayesian
nonparametric extension of Wedderburn’s quasi-likelihood, using Bayesian
additive regression trees to model the mean function. Here, the analyst posits
only a structural relationship between the mean and variance of the outcome.
We show that this approach provides a unified, computationally efficient,
framework for extending Bayesian decision tree ensembles to many new settings,
including simplex-valued and heavily heteroskedastic data. We also introduce
Bayesian strategies for inferring the dispersion parameter of the quasi-
likelihood, a task which is complicated by the fact that the quasi-likelihood
itself does not contain information about this parameter; despite these
challenges, we are able to inject updates for the dispersion parameter into a
Markov chain Monte Carlo inference scheme in a way that, in the parametric
setting, leads to a Bernstein-von Mises result for the stationary distribution
of the resulting Markov chain. We illustrate the utility of our approach on a
variety of both synthetic and non-synthetic datasets.

Key words: Bayesian additive regression trees; decision tree ensembles;
nonparametric regression; robust statistics; Markov chain Monte Carlo.

##  1 Introduction

We consider the problem of performing Bayesian inference when the analyst
would like to avoid specifying a model for the data generating mechanism. To
perform fully-Bayesian inference, one typically begins by specifying (i) a
prior distribution $\pi(\theta)$ for the parameters $\theta$ and (ii) a
likelihood $L(\theta\mid\mathcal{D})$ for the data $\mathcal{D}$. While a
common complaint is that inferences might be sensitive to the choice of
$\pi(\theta)$, a more consequential problem is that it is rarely the case that
one has a strong justification for $L(\theta\mid\mathcal{D})$ either. To make
inferences robust to misspecification of $L(\theta\mid\mathcal{D})$, many
recent works have aimed to bypass the specification of the likelihood
entirely, and instead model the information contained in the data using, for
example, moment equations (Siddhartha Chib and Simoni,, 2018; Yin,, 2009;
Schennach,, 2005) or loss functions (Lyddon et al.,, 2019; Bissiri et al.,,
2016).

In the context of regression, one possible solution is to specify a _quasi-
likelihood_ $Q(\theta)$ (Nelder and Pregibon,, 1987; McCullaugh and Nelder,,
1989), which only requires specifying a mean function
$\mathbb{E}_{\theta}(Y_{i}\mid X_{i}=x)=\mu_{\theta}(x)$ and a mean-variance
relation that relates $\operatorname{Var}_{\theta}(Y_{i}\mid X_{i}=x)$ to
$\mu_{\theta}(x)$. A _quasi-posterior_ for $\theta$ can then be derived by
replacing $L(\theta\mid\mathcal{D})$ with $Q(\theta)$ in Bayes rule, i.e.,
$\widetilde{\pi}(\theta)\propto\pi(\theta)\,Q(\theta)$. Despite this, quasi-
likelihood methods have seen limited attention in the Bayesian literature;
notable exceptions include Agnoletto et al., (2023), Lin, (2006) and several
works of Ventura (Greco et al.,, 2008; Ventura et al.,, 2010; Ventura and
Racugno,, 2016). This may be because most common problems solved by quasi-
likelihood — specifically, accounting for overdispersion for count and
binomial type data — have alternate Bayesian solutions that instead specify a
full data generating mechanism.

This manuscript revisits the quasi-likelihood, primarily for nonparametric
Bayesian inference, and argues that it offers a convenient, robust, and
computationally efficient alternative to the use of full probabilistic models
for routine data analysis problems. We make the following contributions:

  1. 1.

We develop nonparametric models for the mean function $\mu_{\theta}(x)$ using
Bayesian additive regression tree (BART) models, and show that these
approaches are computationally efficient and perform very well on both real
and synthetic data. This effectively reduces our modeling assumptions to a
single assumption about the mean-variance relation. Computations proceed via
simple modifications to the Bayesian backfitting algorithms of Chipman et al.,
(2010) and Murray, (2021). As a byproduct, we obtain extensions of BART to
both simplex-valued outcome data and structured variance models that were
previously unavailable.

  2. 2.

A limitation of using $Q(\theta)$ for inference is that it cannot be used
directly for inferences about the _dispersion parameter_ $\phi$ of the model
(Davidian and Carroll,, 1988; Nelder and Pregibon,, 1987), a parameter that is
analogous to the error variance in linear regression. We address this problem
for Bayesian methods by introducing three different schemes for updating
$\phi$ during Markov chain Monte Carlo (MCMC) based inference, and further
extend our approaches to allow for estimation of the mean-variance relation as
well.

  3. 3.

A theoretical issue is that the best-performing methods for updating $\phi$
during MCMC correspond to updates from “full conditionals” that are not
compatible with the quasi-likelihood (see Gelman and Speed,, 1993, for a
discussion of compatibility). We therefore give theoretical justification to
our approach in the parametric setting by establishing a Bernstein-von Mises
result for our quasi-posterior, with the quasi-posterior of the regression
parameters being asympototically normal with a correctly-calibrated variance.

We illustrate the use of our Bayesian nonparametric quasi-likelihood methods
on synthetic data and data from the Medical Expenditure Panel Survey (MEPS),
which is an ongoing panel survey of the use of the healthcare system in the
United States. On synthetic data, we find that BART-based quasi-likelihood
methods perform extremely well and are no more computationally intensive to
use than the original methods of Chipman et al., (2010). On the MEPS dataset,
we confirm existing results in the literature that the variance in medical
expenditures scales roughly as $V(\mu)=\mu^{\kappa}$ for some
$\kappa\in(1,2)$, and further use our techniques to understand which variables
are most material to forecasting an individual’s medical expenditures.

###  1.1 Notation and Literature Review

We consider data $\mathcal{D}=(\bm{Y},\bm{X})$ where
$\bm{Y}=(Y_{1},\ldots,Y_{N})^{\top}$ denotes the outcome of interest and
$\bm{X}\in\mathbb{R}^{N\times P}=(X_{1},\ldots,X_{N})^{\top}$ is a design
matrix of covariates, with each $X_{i}\in\mathbb{R}^{P}$. Associated to each
outcome $Y_{i}$ is a positive weight $\omega_{i}$; this allows for us to
accommodate binomial-type data with differing counts per observation,
heteroskedastic data with known weights, and survey sampling weights. We
define $\mathbb{E}_{\theta}(Y_{i}\mid X_{i}=x)=\mu_{\theta}(x)$ and assume
that the mean and variance are related through the _mean-variance relation_
$\operatorname{Var}_{\theta}(Y_{i}\mid
X_{i},\omega)=\phi\,V_{\eta}\\{\mu_{\theta}(X_{i})\\}/\omega_{i}$; for
example, for Poisson data we have $V_{\eta}(\mu)=\mu$ while for gamma
distributed data we have $V_{\eta}(\mu)=\mu^{2}$. Here, $\theta$ denotes a
collection of all of the model parameters; it consists of a regression
function $r(\cdot)$ that determines $\mu_{\theta}(\cdot)$, a dispersion
parameter $\phi>0$ and, possibly, a parameter $\eta$ for the _variance
function_ $V_{\eta}$.

Quasi-likelihood (QL) methods were introduced by Wedderburn, (1974) shortly
after the introduction of generalized linear models (GLMs) by Nelder and
Wedderburn, (1972). Rather than defining an exponential dispersion family as
done with GLMs, we use the specified mean and variance functions to define a
_quasi-likelihood function_

$\displaystyle q_{\theta}(y\mid x,\omega)=\exp\left\\{-\frac{\omega
D\\{y,\mu_{\theta}(x)\\}}{2\phi}\right\\}\qquad\text{and}\qquad
D(y,\mu)=2\int_{\mu}^{y}\frac{y-t}{V_{\eta}(t)}\ dt,$

where $D(y,\mu)$ is referred to as the _quasi-deviance_. In the special case
where $[Y_{i}\mid X_{i},\omega_{i},\theta]$ is modeled with an exponential
dispersion family, it can be shown that the quasi-likelihood coincides with
the associated density up-to a term that does not depend on $\mu_{\theta}(x)$.
This definition is motivated by the fact that, when the mean and variance
functions are correctly specified, the quasi-likelihood mimics the usual score
and information equations

$\displaystyle\begin{split}\mathbb{E}_{\theta}\left\\{\frac{\partial}{\partial\mu}\log
q_{\theta}(Y_{i}\mid X_{i},\omega_{i})\mid
X_{i}\right\\}&=0\qquad\text{and}\qquad\\\
\operatorname{Var}_{\theta}\left\\{\frac{\partial}{\partial\mu}\log
q_{\theta}(Y_{i}\mid X_{i},\omega_{i})\mid
X_{i}\right\\}&=-\mathbb{E}_{\theta}\left\\{\frac{\partial^{2}}{\partial\mu^{2}}\log
q_{\theta}(Y_{i}\mid X_{i},\omega_{i})\mid X_{i}\right\\}.\end{split}$ (1)

Because these equations are effectively the same as those used to justify
likelihood-based inference for GLMs, it is unsurprising that the quasi-
likelihood $Q(\theta)=\prod_{i}q_{\theta}(Y_{i}\mid X_{i},\omega)$ behaves in
many ways like a genuine likelihood function even when $q_{\theta}(y\mid
x,\omega)$ is not a valid density; for example, the quasi-likelihood can be
used to perform likelihood ratio tests and construct confidence intervals as
though it corresponded to a genuine likelihood (McCullaugh and Nelder,, 1989,
Chapter 9).

Despite the appeal of quasi-likelihood methods, Bayesian variants are
surprisingly underdeveloped. One possible strategy for Bayesian inference is
to define a _quasi-posterior_ based on data as
$\pi(\theta\mid\mathcal{D})=Q(\theta)\,\pi(\theta)/\int
Q(\theta)\,\pi(\theta)\ d\theta$. In addition to having an interpretation in
terms of quasi-likelihood, following Bissiri et al., (2016) the quasi-
posterior can be justified as giving a coherent update of $\theta$ based on
the prior and the loss function $D(y,\mu)$. When $\phi$ and $\eta$ are fixed
by the prior, this approach works quite well; for example, Agnoletto et al.,
(2023) prove a Bernstein-von Mises result when the mean and variance are
correctly specified and $\phi$ is known.

###  1.2 Dispersion and Variance Function Estimation

A complication of the use of the quasi-posterior is that the quasi-likelihood
cannot be used to infer $\phi$ (McCullaugh and Nelder,, 1989, Section 9.6).
The reason for this is that quasi-likelihoods and genuine likelihoods only
agree up to a constant $c(y,\phi)$, and while this constant is immaterial to
estimating $\mu_{\theta}(x)$ it is very important for estimating $\phi$. While
there have been some attempts at specifying the form of $c(y,\phi)$ for quasi-
likelihood methods to allow estimation of $\phi$ (Nelder and Pregibon,, 1987),
these solutions may not lead to consistent estimation of $\phi$ (Davidian and
Carroll,, 1988).

The most common solution, and the default in most software, is to replace
$\phi$ with a moment-based estimator

$\displaystyle\widehat{\phi}=\frac{1}{N-P}\sum_{i=1}^{N}\frac{\omega_{i}(Y_{i}-\widehat{\mu}_{i})^{2}}{V(\widehat{\mu}_{i})},$
(2)

where the $\widehat{\mu}_{i}$’s are the maximum likelihood estimators of the
$\mu_{\theta}(X_{i})$’s and $P$ is the number of parameters in the
$\mu_{\theta}(\cdot)$ model; importantly, the $\widehat{\mu}_{i}$’s do not
depend on $\phi$. This strategy is effective in parametric non-Bayesian
settings, but it is not clear how to apply it in Bayesian settings where it is
unclear what to replace $N-P$ with and estimation of the
$\mu_{\theta}(X_{i})$’s is no longer decoupled with estimation of $\phi$.

To aaddress this, we introduce three methods — extended quasi-posteriors
(EQPs), pseudo-likelihood posteriors (PLPs), and Bayesian bootstrapped quasi-
likelihood (BBQ) — for estimating $(\phi,\eta)$ during MCMC-based inference.
As there are subtle issues with the compatibility of the PLP and BBQ methods
with the quasi-posterior, we also give some theoretical justification for
these methods.

###  1.3 Outline

In Section 2, we describe further our Bayesian nonparametric quasi-likelihood
approach, discuss its extension to multivariate quasi-likelihood (MQL), and
review the use of Bayesian additive regression trees (BART) to flexibly model
$\mu_{\theta}(x)$. In Section 3 we introduce our approaches to estimating the
dispersion parameter and variance function, and in Section 4 we give
Frequentist justification to these approaches. In Section 5 we provide
illustrations of the use of Bayesian nonparametric quasi-likelihood on both
the MEPS dataset and synthetic data. We close in Section 6 with a discussion.

##  2 Bayesian Nonparametric Quasi-Likelihood

Let $\theta=(r,\phi,\eta)$ denote the parameter vector under the moment
assumptions $g\\{\mu_{\theta}(x)\\}=r(x)$ and
$\operatorname{Var}_{\theta}(Y_{i}\mid
X_{i},\omega_{i})=\phi\,V_{\eta}\\{\mu_{\theta}(X_{i})\\}$. The monotonically
increasing function $g(\cdot)$ is referred to as the _link function_
(McCullaugh and Nelder,, 1989), and it is assumed to be known. Temporarily, we
will also assume that $\phi$ and $\eta$ are known.

Our point of departure is to use the quasi-likelihood to construct a _quasi-
posterior_ distribution for $r$:

$\displaystyle\widetilde{\Pi}(dr\mid\phi,\eta,\mathcal{D})=\frac{Q(\theta)\,\Pi(dr)}{\int
Q(\theta)\,\Pi(r)\ dr}\qquad\text{where}\qquad
Q(\theta)=\exp\left\\{-\frac{1}{2\phi}\sum_{i}\omega_{i}D\\{Y_{i},\mu_{\theta}(X_{i})\\}\right\\}.$

To remove restrictive assumptions on $\mu_{\theta}(x)$, we specify a
nonparametric prior on $r$; we focus here on Bayesian additive regression
trees (BART), but in principle there is no obstruction to modeling $r(x)$
with, for example, a Gaussian process (Rasmussen and Williams,, 2005).

Model | $\omega$ | $\mu$ | $V(\mu)$
---|---|---|---
Quasi-binomial | $n_{i}$ | $(1+e^{r})^{-1}$ | $\mu(1-\mu)$
Quasi-Poisson | 1 | $e^{r}$ | $\mu$
Quasi-gamma | 1 | $e^{-r}\text{ or }e^{r}$ | $\mu^{2}$
Quasi-power | 1 | $e^{r}$ | $\mu^{k}$
Quasi-multinomial | $n_{i}$ | $\frac{e^{r_{j}}}{\sum_{k}e^{r_{k}}}$ | $D_{\mu}-\mu\mu^{\top}$
Table 1: Quasi-likelihood and models for which we can provide simple MCMC
algorithms when $r(x)$ is modeled with BART. Model gives the quasi-likelihood
model name, $\omega$ denotes the form of the weight, $\mu$ denotes the formula
for $\mu_{\theta}(X_{i})$ in terms of $r(X_{i})$, and $V(\mu)$ denotes the
variance function. For the quasi-multinomial,
$D_{\mu}=\operatorname{diag}(\mu_{1},\ldots,\mu_{K})$ where $K$ is the
dimension of $Y_{i}$.

In Table 1 we provide a list of models that we will consider, with the
different choice of weights $\omega_{i}$, choice of $g^{-1}$ to relate
$\mu_{\theta}(x)$ to $r(x)$, and variance function. Before proceeding, we
briefly describe the common use cases for each of these models. The quasi-
binomial and quasi-Poisson models are generally used to address overdispersion
(see, for example McCullaugh and Nelder,, 1989, Chapter 4.5 and 5.5) for
binomial-type and count-type data respectively, as the binomial and Poisson
distributions directly imply that $\phi=1$; additionally, the beta
distribution satisfies the quasi-binomial variance relation, and so the quasi-
binomial model can also be used to perform beta regression. The quasi-gamma
model is in fact formally equivalent to the usual gamma regression model, but
we continue with the term “quasi-gamma” to make clear first that we are not
assuming that the underlying data is gamma distributed and second that the
dispersion parameter $\phi$ is not estimated using the gamma likelihood;
beyond the gamma distribution, the quasi-gamma model naturally addresses scale
regression models of the form $Y_{i}=e^{r(X_{i})}\,\epsilon_{i}$ with
$\mathbb{E}(\epsilon_{i})=1$ and $\operatorname{Var}(\epsilon_{i})=\phi$, and
so is valid for many settings of interest where the data are not gamma
distributed. The quasi-power model generalizes both the quasi-Poisson and
quasi-gamma models to allow for an arbitrary power in the variance, and is
appropriate for the medical expenditure data we analyze here; these are
sometimes referred to as _Tweedie models_ (Smyth and Verbyla,, 1999). Finally,
the quasi-multinomial model is appropriate both for multinomial-type data to
handle overdispersed outcomes, but is also appropriate for simplex-valued
outcome data such as Dirichlet-outcome data, where outcomes will tend to be
_underdispersed_ relative to the multinomial distribution.

###  2.1 Review Bayesian Additive Regression Tree

Conveniently, when the _Bayesian additive regression tress_ (BART) framework
is used to model $r(x)$, all of the quasi-likelihood models given in Table 1
are amenable to simple Gibbs sampling algorithms using “Bayesian backfitting”
(Hastie and Tibshirani,, 2000; Chipman et al.,, 2010). The BART we framework
we use sets

$\displaystyle\mu_{\theta}(x)=g^{-1}\\{r(x)\\}\qquad\text{where}\qquad
r(x)=\sum_{t=1}^{T}\operatorname{Tree}(x;\mathcal{T}_{t},\mathcal{M}_{t}).$

The function $\operatorname{Tree}(x;\mathcal{T},\mathcal{M})$ denotes a
_regression tree_ with decision tree $\mathcal{T}$ and leaf node parameters
$\mathcal{M}$, which returns the value $\lambda_{\ell}$ if $x$ is associated
to leaf $\ell$ of tree $\mathcal{T}$. A schematic of such a regression tree is
given in Figure 1, with $x=(x_{1},x_{2})^{\top}\in[0,1]^{2}$. It is apparent
from this schematic that a regression tree on $[0,1]^{P}$ is associated with a
piecewise constant function, with $r(x)$ corresponding to a basis function
expansion of many such decision shallow trees.

$x_{1}\leq 0.6$ $x_{2}\leq 0.4$ $\lambda_{1}$ $\lambda_{2}$ $\lambda_{3}$
TrueFalseTrueFalse$\lambda_{1}$$\lambda_{2}$$\lambda_{3}$$x_{1}$$x_{2}$ Figure
1: Schematic of a regression tree decision tree (left) and the induced step-
function on $[0,1]^{2}$ induced by the tree.

BART is similar to other decision tree ensembling approaches, such as gradient
boosting (Friedman,, 2001), with the main differences being (i) that we place
prior distributions on the $(\mathcal{T}_{t},\mathcal{M}_{t})$’s and (ii) that
inference proceeds via Markov chain Monte Carlo rather than in a greedy
fashion. We favor the use of BART as a default prior for $r(\cdot)$, as it is
highly flexible, performs very well in practice, and is admits “default
priors” that allow the analyst to automatically choose the model
hyperparameters; see, for example, Chipman et al., (2010) or Linero and Yang,
(2018) for the performance of BART on benchmark datasets relative to random
forests, gradient boosting, the lasso. BART also possesses a number of
desirable theoretical properties, such as obtaining near-minimax rates of
posterior concentration within the space of functions that consist primarily
of main effects and low-order interactions (Ročková and van der Pas,, 2020;
Linero and Yang,, 2018). See Hill et al., (2020) and Linero, (2017) for
comprehensive reviews on BART and Bayesian decision tree-based methods.

###  2.2 Prior Specification and Bayesian Backfitting Algorithms

Assuming that $V_{\eta}(\cdot)$ and $\phi$ are fixed a-priori, our approach
only requires specifying a prior for the
$(\mathcal{T}_{t},\mathcal{M}_{t})$’s. We use the default branching process
prior on the $\mathcal{T}_{t}$’s described by Chipman et al., (2010), which
encourages most decision trees to be of depth one or two; we defer the
interested reader to that work for details. These priors for $\mathcal{T}_{t}$
are empirically robust choices, and are taken by default by most BART software
packages (Kapelner and Bleich,, 2016; Dorie,, 2023; Sparapani et al.,, 2021).
In our illustrations and software we also use the Dirichlet splitting rule
prior introduced by Linero, (2018) to allow the model to automatically remove
irrelevant predictors.

The only remaining detail for prior specification is the choice of prior for
the $\mathcal{M}_{t}$’s. We assume that the $\lambda_{t\ell}$’s are iid from
some distribution $\pi_{\lambda}(\lambda)$. In order to obtain a conditionally
conjugate model, we take $\pi_{\lambda}(\lambda)$ to be a
$\log\operatorname{Gam}(a,b)$ distribution for the quasi-binomial, quasi-
Poisson, quasi-gamma, and quasi-multinomial models, with $a$ and $b$ chosen so
that $\mathbb{E}(\lambda)=0$ and
$\operatorname{Var}(\lambda)=\sigma^{2}_{\lambda}$ for some user-specified
value of $\sigma_{\lambda}$. We then follow Chipman et al., (2010) and set
$\sigma_{\lambda}=3/(k\sqrt{T})$ with $k=2$ by default; see Chipman et al.,
(2010) for a justification of this choice, and in particular for the
normalization in the prior by $\sqrt{T}$. Because it is unclear that the
performance of the defaults extend to the new models we propose, we also
recommend using a half-Cauchy distribution with scale $3/(k\sqrt{T})$ as a
prior for $\sigma_{\lambda}$; this allows the model to estimate
$\sigma_{\lambda}$, while encouraging it to be near the default of Chipman et
al., (2010).

The models in Table 1 can be fit via Markov chain Monte Carlo using the
generalized Bayesian backfitting algorithm given by Hill et al., (2020), with
the exception of the quasi-power model. These algorithms work by first
approximately sampling $\mathcal{T}_{t}$ from its conditional distribution
$[\mathcal{T}_{t}\mid\mathcal{T}_{-t},\mathcal{M}_{-t},\phi,\eta,\mathcal{D}]$
using a Metropolis-Hastings algorithm and second sampling the leaf node
parameters $\mathcal{M}_{t}$ from the full conditional distribution
$[\mathcal{M}_{t}\mid\mathcal{T}_{t},\mathcal{T}_{-t},\mathcal{M}_{-t},\phi,\eta,\mathcal{D}]$;
here, $\mathcal{T}_{-t}$ and $\mathcal{M}_{-t}$ denote the decision trees and
leaf node parameters of all of the regression trees except tree $t$. Gibbs
samplers for these models are given in the appendix.

By contrast, the quasi-power model cannot be fit using any simple extensions
of existing BART sampling algorithms. We instead use the Laplace
approximation, which can be computed analytically, to approximate the relevant
marginal likelihoods and full conditionals. After applying the Laplace
approximation, implementing the quasi-power model is no more computationally
intensive than the original BART model. Moreover, as shown by Linero, (2022),
it is easy to use the Laplace approximation to construct a valid proposal in a
Metropolis-Hastings algorithm that holds the target posterior exactly
invariant. Details of this approximation are given in the appendix.

###  2.3 Extension to Multivariate Quasi-Likelihood Models

Defining the quasi-multinomial model in Table 1 requires extending the quasi-
likelihood to the multivariate setting (McCullaugh and Nelder,, 1989, Section
9.3). Consider the pair of moment equations $\mathbb{E}_{\theta}(Y_{i}\mid
X_{i}=x,\omega_{i})=\mu_{\theta}(x)$ and
$\operatorname{Var}_{\theta}(Y_{i}\mid
X_{i}=x,\omega_{i})=\frac{\phi}{\omega_{i}}\,V_{\eta}\\{\mu_{\theta}(x)\\}$
where $Y_{i}$ and $\mu_{\theta}(x)$ are now a $K$-dimensional vectors and
$V_{\eta}(\mu)$ is a $K\times K$ covariance matrix. The _multivariate quasi-
likelihood_ (MQL) is defined as

$\displaystyle q_{\theta}(y\mid
x,\omega)=\exp\left\\{-\frac{\omega\,D\\{y,\mu_{\theta}(x)\\}}{2\phi}\right\\}\qquad\text{where}\qquad
D(y,\mu)=-2\int_{\mu}^{y}(y-t)^{\top}V_{\eta}(t)^{-1}\ dt,$

where both $\mu$ and $y$ are $K$-dimensional vectors, $V(t)$ is a $K\times K$
covariate matrix, and the integral is a path-independent line integral.

The MQL model we consider takes $g\\{\mu_{\theta}(X_{i})\\}=r(X_{i})$ where
$r(x)=\sum_{t=1}^{T}\operatorname{Tree}(x;\mathcal{T}_{t},\mathcal{M}_{t})$ is
now a sum of _multivariate decision trees_ ; see, for example, Linero et al.,
(2020) for how multivariate decision trees can be combined using BART to model
multivariate outcomes.

We focus on the quasi-multinomial model, which takes
$V(\mu)=D_{\mu}-\mu\mu^{\top}$ where
$D_{\mu}=\operatorname{diag}(\mu_{1},\ldots,\mu_{K})$ is a diagonal matrix.
The quasi-likelihood for this model is

$\displaystyle q_{\theta}(y\mid
x,n)=\frac{\exp\left\\{z^{\top}r(x)/\phi\right\\}}{[\sum_{k}\exp\\{r_{k}(x)\\}]^{n/\phi}}$

where $\omega=n$ corresponds to the sample size of a multinomial count vector
$z$ and $y=z/n$ is a vector of proportions. This model reduces to the usual
multinomial likelihood in the special case when $\phi=1$, and can be used to
also model Dirichlet response data when $n=1$ and $\phi<1$. Note that, when
applied to multinomial data, the differences $r_{k}(x)-r_{j}(x)$ correspond to
the log-odds ratio associated to a comparison between category level $k$ and
category level $j$.

For a general MQL model, assuming the $\mu_{i}$’s are known a moment estimator
of $\phi$ is given by

$\displaystyle\widehat{\phi}=\frac{1}{NK}\sum_{i}\omega_{i}(Y_{i}-\mu_{i})^{\top}V(\mu_{i})^{-1}(Y_{i}-\mu_{i}).$
(3)

In the special case of the quasi-multinomial model, this moment estimator
becomes

$\displaystyle\widehat{\phi}=\frac{1}{N(K-1)}\sum_{i}\sum_{k}\frac{n_{i}(Y_{ik}-\mu_{ik})^{2}}{\mu_{ik}}$

with the modification of $K$ to $K-1$ due to the loss of a degree of freedom
from the constraint $\sum_{k}Y_{ik}=1$; it is straight-forward to verify that
$\widehat{\phi}$ is unbiased for $\phi$.

To obtain a Bayesian nonparametric quasi-multinomial model, we set $r=r(x)$ in
the above model where $r(x)$ is a sum of multivariate decision trees
$r(x)=\sum_{t}\operatorname{Tree}(x;\mathcal{T}_{t},\mathcal{M}_{t})$ where
$\operatorname{Tree}(x;\mathcal{T},\mathcal{M})=(\lambda_{\ell
1},\ldots,\lambda_{\ell K})^{\top}$ if $x$ is associated to leaf node $\ell$
of $\mathcal{T}$. For identifiability purposes in the quasi-multinomial model,
it is customary to set $r_{1}(x)\equiv 0$. However, as the individual
functions $r_{k}(x)$ are not of direct interest, and to ensure that the model
is symmetric across all categories, we place no restrictions on the
$\lambda_{\ell k}$’s. The $\lambda_{\ell k}$’s are given the same log-gamma
prior described in Section 2.2, with the default
$\sigma_{\lambda}\sim\text{half-Cauchy}\\{0,3/(k\sqrt{2T})\\}$ with $k=2$ by
default. The choice $3/k\sqrt{2T}$ is chosen to represent a prior belief that,
for each $x$ the log-odds ratio $r_{k}(x)-r_{j}(x)$ is approximately normally
distributed with mean $0$ and standard deviation $3/k$.

##  3 Inference for the Dispersion Parameter

We now describe some strategies for estimating $\phi$ and discuss their
merits. For the sake of exposition, the updates are roughly ordered from least
to most recommended.

#### Extended Quasi-Posteriors (EQP)

Introduced by Wedderburn, (1974), the _extended quasi-likelihood_ (EQL)
augments $Q(\theta)$ with an additional term $\sqrt{2\pi\phi V_{\eta}(Y_{i})}$
to obtain

$\displaystyle Q^{+}(\theta)=\prod_{i}\frac{1}{\sqrt{2\pi\phi
V_{\eta}(Y_{i})}}\exp\left(-\frac{\omega_{i}D\\{Y_{i},\mu_{\theta}(X_{i})\\}}{2\phi}\right).$

The extended quasi-posterior is then given by
$\Pi(dr,d\phi,d\eta\mid\mathcal{D})\propto
Q^{+}(\theta)\,\Pi(dr,d\phi,d\eta)$. Given a gamma prior
$\phi^{-1}\sim\operatorname{Gam}(a,b)$, the EQP is convenient in that the full
conditional for $\phi$ is given by
$\phi^{-1}\sim\operatorname{Gam}\\{a+N/2,b+\sum_{i}\omega_{i}D\\{Y_{i},\mu_{\theta}(X_{i})\\}/2\\}$.
The EQP is not recommended because it generally does not produce a consistent
estimator of $\phi$ (Davidian and Carroll,, 1988). The primary advantage of
EQL is that it has none of the compatibility issues of the PLP and BBQ
approaches discussed below.

#### Pseudo-Empirical Bayes

A simple strategy is to use the moment estimator (2) obtained by fitting the
model with (say) $\phi=1$, and then refit the model with this new value of
$\phi$. This might then be iterated: we refit the model with the new value of
$\phi$, apply (2), and repeat this some modest number of times. We do not
recommend this approach for Bayesian nonparametric models because it is not
computationally efficient and we lack of obvious replacement for $P$ in (2).
Instead, inspired by the pseudo-likelihood (PL) approach of Davidian and
Carroll, (1988), we propose a pseudo-empirical Bayes strategy, which
iteratively updates $(\phi,\eta)$ by setting

$\displaystyle\begin{split}(\phi_{t+1},\eta_{t+1})&\leftarrow\arg\max_{\phi,\eta}\int\text{PL}(\phi,\eta)\
\
\widetilde{\Pi}(r\mid\phi_{t},\eta_{t},\mathcal{D})\qquad\text{where}\qquad\\\
\text{PL}(\phi,\eta)&=\prod_{i}\frac{1}{\sqrt{2\pi\phi\,V_{\eta}\\{\mu(X_{i})\\}/\omega_{i}}}\exp\left[-\frac{\omega_{i}\\{Y_{i}-\mu(X_{i})\\}^{2}}{2\,\phi\,V_{\eta}\\{\mu(X_{i})\\}}\right].\end{split}$
(4)

To approximately evaluate the integral, one can use empirical Bayes Gibbs
sampling (Casella,, 2001). This can be iterated a modest number of times to
obtain plausible values of $(\phi,\eta)$. This is computationally expensive
and does not provide uncertainty quantification for $(\phi,\eta)$, and so
still not ideal.

#### Pseudo-Likelihood Posteriors (PLP)

The pseudo-likelihood posterior (PLP) approach extends the pseudo-empirical
Bayes approach to obtain a posterior distribution for $(\phi,\eta)$, and is
both computationally simpler and obtains (uncalibrated) uncertainty
quantification. The main idea is to run Bayesian backfitting on $r(x)$, but
include updates for $(\phi,\eta)$ based on the partial likelihood; that is, we
use a two-stage Gibbs sampler:

  1. 1.

Sample $r$ from $\widetilde{\Pi}(dr\mid\phi,\eta,\mathcal{D})\propto
Q(\theta)\,\Pi(dr)$.

  2. 2.

Sample $(\phi,\eta)$ from $\widetilde{\Pi}(d\phi,d\eta\mid
r,\mathcal{D})\propto\text{PL}(\phi,\eta)\,\Pi(d\phi,d\eta)$.

As with the EQL approach, under the gamma prior
$\phi^{-1}\sim\operatorname{Gam}(a,b)$, the PL approach results in a simple
full conditional
$\phi^{-1}\sim\operatorname{Gam}(a+N/2,b+\sum_{i}Z_{i}^{2}/2)$, where
$Z_{i}=\frac{\omega_{i}^{1/2}\\{Y_{i}-\mu_{\theta}(X_{i})\\}}{V\\{\mu_{\theta}(X_{i})\\}^{1/2}}$
is the _standardized residual_. The PL approach generally produces good point
estimates of $(\phi,\eta)$ and uncertainty quantification for $r(\cdot)$, but
may not produce good uncertainty quantification for $(\phi,\eta)$ and so
should be used with caution when uncertainty in the variance function is of
intrinsic interest. The PL approach also samples from incompatible
conditionals for $r$ and $(\phi,\eta)$ in the sense that there is no joint
distribution for which the update distributions are the conditionals.

#### Bayesian Bootstrap Quasi-Likelihood (BBQ)

A simple, robust, approach to sampling $\phi$ is to use an approach inspired
by the loss-likelihood bootstrap (Lyddon et al.,, 2019). Again let
$Z_{i}=\frac{\omega_{i}^{1/2}\\{Y_{i}-\mu_{\theta}(X_{i})\\}}{V\\{\mu_{\theta}(X_{i})\\}^{1/2}}$
denote the standardized residual and note that
$\mathbb{E}_{\theta}(Z_{i}^{2})=\phi$. We model the marginal distribution of
the $Z_{i}$’s using a Bayesian bootstrap
$F_{Z}(dz)=\sum_{i=1}^{N}p_{i}\,\delta_{Z_{i}}(dz)$ with the full conditional
$p\sim\operatorname{Dirichlet}(1,\ldots,1)$ and set

$\displaystyle\phi=\sum_{i=1}^{N}p_{i}\,Z_{i}^{2}.$ (5)

MCMC proceeds as with the PL approach, except that $\phi$ is updated by
sampling $p\sim\operatorname{Dirichlet}(1,\ldots,1)$ and using (5). This
approach can further accommodate estimation of $\eta$ by replacing (5) with
the update
$(\phi,\eta)=\arg\max_{\phi,\eta}\sum_{i}p_{i}\log\operatorname{Normal}\big{(}Y_{i}\mid\mu_{\theta}(X_{i}),\phi\,V_{\eta}\\{\mu_{\theta}(X_{i})\\}/\omega_{i}\big{)}$.
The advantage of the BBQ approach over the PLP approach is that we intuitively
expect the use of the bootstrap to provide better calibrated uncertainty
quantifification for $(\phi,\eta)$. The downsides of the BBQ approach is that
the conditionals of $r$ and $(\phi,\eta)$ are again incompatible and that it
may be difficult to jointly optimize $(\phi,\eta)$ in a robust fashion.

##  4 Theoretical Properties

We now examine the theoretical properties of the PL approach (see Agnoletto et
al.,, 2023, for results for the quasi-posterior with known $\phi$). Proofs of
these results can be found in the appendix. For simplicity, we consider the
parametric setting with $r(x)=x^{\top}\beta$ and assume $\eta$ to be known.
Let $g(\phi\mid\beta)$ denote conditional density that $\phi$ is sampled from
in our Gibbs sampler and let $\widetilde{\pi}(\beta\mid\phi)\propto
Q(\beta,\phi)\,\pi(\beta)$. We restrict attention to the ideal _two-step Gibbs
sampler_ :

  1. i.

Sample $\beta\sim\widetilde{\pi}(\beta\mid\phi)$.

  2. ii.

Sample $\phi\sim g(\phi\mid\beta)$ where $g$ is the PL update density for
$\phi$.

We show that the sequence $\beta^{(1)},\beta^{(2)},\cdots$ is a Markov chain
with a stationary distribution that is expressible as a mixture
$\widetilde{\pi}(\beta)=\int\widetilde{\pi}(\beta\mid\phi)\,m(\phi)\ d\phi$.
The mixing distribution $m(\phi)$ is also the marginal distribution of $\phi$
in the stationary distribution, and is well-behaved in the sense that it
concentrates near $\phi_{0}$, which effectively causes
$\widetilde{\pi}(\beta)$ to be well-approximated by a
$\widetilde{\pi}(\beta\mid\phi_{0})$ density, leading to a Bernstein-von Mises
result for $\widetilde{\pi}(\beta)$ due to existing results of Agnoletto et
al., (2023) and Miller, (2021).

A minor inconvenience with the BBQ and PL approaches is that there is
typically no joint distribution $\widetilde{\pi}(\beta,\phi)$ whose
conditional distributions are $\widetilde{\pi}(\beta\mid\phi)$ and
$g(\phi\mid\beta)$; a necessary condition for this to occur is that
factorization
$\widetilde{\pi}(\beta\mid\phi)/g(\phi\mid\beta)=u(\beta)\,v(\phi)$ holds for
some functions $u(\cdot),v(\cdot)$ (Arnold and Press,, 1989, Theorem 4.1),
which is not the case for any of the models we have discussed. The following
result states that, despite these challenges, the chain has a well-defined
stationary distribution.

######  Theorem 1.

Suppose that Condition A and Condition P in Appendix E hold. Then the two-step
Gibbs sampler initialized at any point $\phi>0$ has a unique stationary
distribution with density
$\widetilde{\pi}(\beta,\phi)=\widetilde{\pi}(\beta)\,g(\phi\mid\beta)$ where
$\widetilde{\pi}(\beta)=\int\widetilde{\pi}(\beta\mid\phi)\,m(\phi)\ d\phi$
for some mixing distribution $m(\phi)$. Additionally,
$m(\phi)=\int\widetilde{\pi}(\beta)\,g(\phi\mid\beta)\,d\beta$.

With the issue of the existence of the stationary distribution dispensed with,
we now present our Bernstein-von Mises result.

######  Theorem 2.

Under Condition P, Condition A, Condition F, Conditions D1 – D4, and Condition
G listed in the Appendix hold, we have

$\displaystyle\int|\widetilde{\pi}(\beta)-N\\{\beta\mid\widehat{\beta},\phi_{0}N^{-1}H(\beta_{0})^{-1}\\}|\
d\beta\to 0\qquad\text{almost surely},$

where $H(\beta_{0})$ is given in Condition D and
$\phi_{0}N^{-1}H(\beta_{0})^{-1}$is the asymptotic covariance matrix of the
maximum quasi-likelihood estimator
$\widehat{\beta}=\arg\min_{\beta}\sum_{i}\omega_{i}D\\{Y_{i};g^{-1}(X_{i}^{\top}\beta)\\}$.
Additionally, $m(\phi)$ converges in distribution to a point-mass at
$\phi_{0}$ almost surely.

When the model is correctly specified, Theorem 2 implies correct asymptotic
coverage of posterior credible intervals for $\beta$. While we do not prove
this, under misspecification we generally expect that a version of Theorem 2
will hold, but $\phi_{0}\,H(\beta_{0})^{-1}/N$ will no longer be the
asymptotic covariance matrix of $\widehat{\beta}$, and so inferences will not
be correctly calibrated.

######  Remark 1.

Theorem 2 is proved for the PLP rather than the BBQ posterior simply because
it is more analytically tractable, however we expect Theorem 2 to also hold
for the BBQ posterior. As argued by Davidian and Carroll, (1988), the extended
quasi-likelihood does not lead to consistent estimates of $(\phi,\eta)$; in
our own experiments, it was instead observed that the EQP concentrates around
$\beta_{0}$ but does not possess good Frequentist properties, and that the
posteriors for $\phi$ and $\eta$ do not concentrate around their true values.
This is unfortunate because the EQP avoids the problems listed above:
inference for all components of $\theta$ is based on the single likelihood
$Q^{+}(\theta)$, and so the stationary distribution of the Markov chain will
exist and be equal to $Q^{+}(\theta)\,\Pi(d\theta)/\int
Q^{+}(\theta)\,\Pi(d\theta)$ so long as the denominator is finite. Because of
these issues, we do not study the theoretical properties of the EQP posterior.

######  Remark 2.

The regularity conditions used in the appendix are optimized for ease of proof
rather than generality, and our results should be regarded primarily as a
proof-of-concept. The most substantial simplification we make is that
$\phi_{0}$ is assumed to lie in an interval $[a,b]$, and $\beta_{0}$ in some
bounded set, with both sets known a-priori.

##  5 Illustrations

###  5.1 Experiments With Parametric Models

We first provide some numerical support for the Bernstein-von Mises results
given in Section 4. That a Bernstein-von Mises result will hold does not
follow from our theory, even for parametric models, because we usually cannot
use the ideal two-step Gibbs sampler. We consider data generated according to
the quasi-Poisson model with

$\displaystyle Y_{i}=\phi\,Z_{i}\qquad\text{where}\qquad
Z_{i}\sim\operatorname{Poisson}\\{\mu_{\theta}(X_{i})/\phi\\}$

and $\mu_{\theta}(x)=\exp(\beta_{0}+\beta_{1}\,x)$. Data generated according
to this model satisfies the quasi-Poisson relation
$\operatorname{Var}(Y_{i}\mid X_{i}=x)=\phi\,\mu_{\theta}(x)$. We considered a
flat prior $\pi(\beta)\propto 1$, $X_{i}\sim\operatorname{Normal}(0,1)$,
$N=2000$, and use the BBQ approach to quantify uncertainty in $\phi$.
Hamiltonian Monte Carlo (HMC, Neal,, 2011) was used to sample from the quasi-
posterior.

Figure 2: Joint and marginal posterior densities of the parameters from the
quasi-Poisson simulation in Section 5.1. The quasi-posterior under the BBQ
prior is given by “BBQ” while the asymptotic distribution with mean equal to
the MLE and variance equal to the inverse Fisher information corresponds to
“glm.”

The posterior distribution is displayed in Figure 2 and compared with the
asymptotic normal distribution obtained from standard quasi-likelihood theory
using the pilot estimator
$\widehat{\phi}=(N-2)^{-1}\sum_{i}(Y_{i}-\widehat{\mu}_{i})^{2}/\widehat{\mu}_{i}$.
We see strong agreement between the quasi-posterior and the asymptotic normal
distribution, giving some further evidence that the good theoretical
properties of the two-step Gibbs sampler likely extend to methods that replace
sampling $\widetilde{\pi}(\beta\mid\phi)$ with sampling a Markov transition
density that merely leaves $\widetilde{\pi}(\beta\mid\phi)$ invariant.

Next, we consider a quasi-power model with $V_{\eta}(\mu)=\mu^{\kappa}$. This
model contains the quasi-Poisson $(\eta=1)$ and quasi-gamma $(\eta=2)$ as
submodels. We set $\mu_{\theta}(x)=\exp(x^{\top}\beta)$ with
$\beta_{j}=1.2/\sqrt{5}$ for the $\kappa=1$ setting and $\beta_{j}=2/\sqrt{5}$
for the $\kappa=2$ setting; these values were chosen to highlight differences
between the quasi-gamma and quasi-Poisson models without being unrealistic. We
evaluate the following models:

  * •

Quasi-Poisson: a quasi-Poisson GLM fit using the glm function in R via maximum
quasi-likelihood.

  * •

Quasi-gamma: a quasi-gamma GLM fit using the glm function in R via maximum
quasi-likelihood.

  * •

Bayesian Bootstrap Poisson: samples of $\beta$ are collected by first sampling
$p\sim\operatorname{Dirichlet}(1,\ldots,1)$ and then solving the equation
$\sum_{i}\frac{p_{i}(Y_{i}-\mu_{\beta}(X_{i}))}{\mu_{\beta}(X_{i})}\,X_{i}=0$.

Figure 3: Width of 95% intervals (top), root mean-squared error (middle), and
bias (bottom) results for the power-variance experiment. BBP denotes the
Bayesian Bootstrap Poisson model, BBQ denotes Bayesian Bootstrap Quasi-
likelihood, QG denotes quasi-gamma, and QP denotes quasi-Poisson. Coverage
rates are given above the widths in the top row.

Figure 3 displays the root mean squared error, estimated bias, and interval
width/coverage associated to nominal 95% interval estimates. We see that the
BBQ approach is the only approach that performs uniformly well across all
settings, and that the performance of all of the methods is not greatly
impacted by $\phi$. When $\kappa=1$ (quasi-Poisson holds) we see that the BBQ,
Quasi-Poisson, and BB Poisson methods are all roughly equivalent, while the
quasi-gamma estimates have higher mean-squared errors and have substantially
wider intervals to obtain coverage rates just below the nominal. On the other
hand, when $\kappa=2$ (quasi-gamma holds), the methods based on the Poisson
likelihood perform poorly in terms of interval width, coverage, RMSE, and
bias, while the BBQ approach essentially matches the performance of the quasi-
gamma model.

###  5.2 Bayesian Nonparametric Quasi-Gamma Model

We next consider a quasi-gamma model where $V(\mu)=\mu^{2}$ under a data
generating mechanism with $\mu_{\theta}(x)=e^{r(x)}$ where
$r(x)=\sin(\pi\,x_{1}\,x_{2})+2(x_{3}-0.5)^{2}+x_{4}+x_{5}/2$; the form of
$r(x)$ is taken from Friedman, (1991), where it was used to illustrate the
benefits of multivariate adaptive regression splines (MARS). For the outcome,
we set $[Y_{i}^{-1}\mid
X_{i}=x]\sim\operatorname{Gam}\\{\alpha,(\alpha-1)\,\mu_{\theta}(x)\\}$ which
corresponds to a quasi-gamma model with $\phi=(\alpha-2)^{-1}$.

We compare the quasi-gamma model to the conventional (misspecified) gamma
model that takes $[Y_{i}\mid
X_{x}=x]\sim\operatorname{Gam}\\{\alpha,\alpha/\mu_{\theta}(x)\\}$. We note
that the gamma and quasi-gamma models differ _only_ in how the parameter
$\phi$ is estimated, with the quasi-gamma model using BBQ and the gamma model
using the posterior derived from the gamma likelihood.

Figure 4: Boxplots of mean squared error of the gamma model relative to the
quasi-gamma model for the experiment described in Section 5.2.

The models were compared in terms of relative mean squared error, defined as
$\frac{\sum_{i}\\{r(X_{i})-\widehat{r}_{\text{G}}(X_{i})\\}^{2}}{\sum_{i}\\{r(X_{i})-\widehat{r}_{\text{QG}}\\}^{2}}$,
where $\widehat{r}_{\text{G}}(x)$ and $\widehat{r}_{\text{QG}}(x)$ are the
Bayes estimates of $r(x)$ under the gamma and quasi-gamma models,
respectively. The results of 200 replications of this experiment under the
settings $N=250$, $P=10$, and $\phi\in\\{0.5,1,2\\}$ are given in Figure 4. We
see that the quasi-gamma model typically outperforms the gamma model on this
problem, with a larger gap existing for larger values of $\phi$. In terms of
practical differences, we found that the gamma model tends to substantially
underestimate $\phi$, while BBQ tends to slightly overestimate $\phi$, with
this accounting for the mild differences in predictive performance.

###  5.3 Bayesian Nonparametric Quasi-Power Model

To assess the reliability of our implementation of the quasi-power model, we
use the quasi-power ground truth $[Y_{i}\mid
X_{i}=x]\sim\operatorname{Gam}\\{\alpha(x),\beta(x)\\}$ where
$\alpha(x)=e^{r(x)/2}/\phi$ and $\beta(x)=e^{-r(x)/2}/\phi$; this corresponds
to a quasi-power model with $\mu_{\theta}(x)=e^{r(x)}$ and $\kappa=1.5$. We
take $r(x)$ to be the same as in Section 5.2.

We compare the quasi-power model to the quasi-Poisson and quasi-gamma models,
which correspond to special cases of the quasi-power model with $\kappa=1$ and
$2$ respectively. We set $N=250$, $P=10$, $\phi\in\\{0.5,1,2\\}$, and
replicated the experiment $200$ times. Models were compared in terms of mean
squared error
$N^{-1}\sum_{i}\\{\mu_{\theta}(X_{i})-\widehat{\mu}(X_{i})\\}^{2}$ where
$\widehat{\mu}(x)$ is the Bayes estimator of $\mu_{\theta}(x)$, as well as the
coverage per-experiment of the pointwise nominal 95% credible intervals for
the $\mu_{\theta}(X_{i})$’s and the average length of these intervals.

Figure 5: Results from simulation study of Section 5.3. Boxplots displays the
squared error, coverage, and interval length associated to $\mu(X_{i})$ for
$i=1,\ldots,N$, averaged over $N$, for the 200 replications.

Results are given in Figure 5. As anticipated, the quasi-power modestly
improves upon the quasi-gamma and quasi-Poisson models in terms of predictive
performance. Interestingly, all models appear to have conservative credible
intervals, with the quasi-Poisson model coming the closest to the nominal
coverage levels.

Figure 6: Mixing of $\phi$ and $\kappa$ for a single replication of the
simulation experiment in Section 5.3 with $N=2500$.

Figure 6 displays the mixing of $\phi$ and $\kappa$ on a single representative
simulation iteration with $N=2500$. We see that $\phi$ and $\kappa$ have
posteriors centered around their true values, and mix reasonable well over
2,000 post-burn-in samples.

###  5.4 Quasi-Multinomial With Dirichlet Outcome Data

Our last synthetic data experiment applies the quasi-multinomial model to
Dirichlet outcome data. We consider the data generating process
$Y_{i}\sim\operatorname{Dirichlet}\\{\rho\,\mu_{1}(X_{i}),\rho\,\mu_{2}(X_{i}),\rho\,\mu_{3}(X_{i})\\}$
where $\mu_{k}(x)\propto\exp\\{r_{k}(x)\\}$ and

$\displaystyle r_{1}(x)$ $\displaystyle=2x_{1}+x_{2},$ $\displaystyle
r_{2}(x)$ $\displaystyle=x_{1}+4x_{2}x_{3}$ $\displaystyle r_{3}$
$\displaystyle=x_{2}+2x_{3}.$

We consider $N=1000$ observations,
$X_{i}\sim\operatorname{Uniform}([0,1]^{5})$, and set $\rho=0.5$, which
corresponds to a quasi-multinomial model with $\phi=2/3$.

Figure 7: Top: comparison of the Bayes estimate ($x$-axis) of $\mu_{1}(X_{i})$
to the true value of $\mu_{1}(X_{i})$ ($y$-axis) for the BART normal
regression model (left) and quasi-Multinomial model (right); the line $y=x$ is
provided for comparison. Bottom: traceplot of $\phi$ for the quasi-Multinomial
model.

Figure 7 gives results for a single simulated dataset, comparing the quasi-
multinomial model to the standard BART regression model implemented in the
dbarts package that directly estimates $\mu_{1}(x)$; results for this single
simulated dataset are typical of replications. We see that the quasi-
multinomial model generally performs better than BART, with no obvious bias
and estimates falling quite close to the lines $x=y$.

Method | RMSE | Width | Coverage
---|---|---|---
BART | 0.08 | 0.28 | 0.92
QMN | 0.05 | 0.20 | 0.95
Table 2: Results of the simulation experiment in Section 5.4

This experiment was replicated 100 times, and results are presented in Table
2. Here, RMSE refers to the RMSE in estimating $\mu_{1}(X_{i})$ given by
$\sqrt{\frac{1}{N}\sum_{i}\\{\mu_{1}(X_{i})-\widehat{\mu}_{1}(X_{i})\\}^{2}}$,
Width refers to the average interval width of the 95% highest posterior
density (HPD) credible interval, and Coverage refers to the average coverage
of the HPD interval. By properly accounting for the presence of
heteroskedasticity we see that the quasi-multinomial model outperforms BART on
all three metrics, with a much lower RMSE being attained and nominal coverage
being achieved with shorter intervals.

###  5.5 The Medical Expenditure Panel Survey

The medical expenditure panel survey (MEPS) is an ongoing survey on the usage
of the healthcare system in the United States by individuals, insurance
companies, and hospitals. We consider as an outcome $Y_{i}$ the total medical
expenditure incurred during 2011 by a subset of 9,426 adult women who incurred
at least some medical expenditure. Parametric quasi-likelihood methods have
been used in the past for this data and it has been demonstrated that the data
is highly heteroskedastic with a mean-variance relation well-described by a
quasi-power relation $V(\mu)=\mu^{\kappa}$ (Natarajan et al.,, 2008, who found
$\kappa\approx 1.5$ on the subset of the data they examined). Our interest is
in determining what factors lead to higher forecasts of medical expenditures,
with variables of interest including demographic variables, measurements of an
individual’s health, and direct measures of the usage of the healthcare
system. Predictors include age, insurance status, race, whether an individual
is a regular smoker, income, self-assessed physical health status, mental
health status, body mass index, number of prescription refills during the
year, education level, number of dentist trips during the year, and indicators
of major health diagnoses (diabetes, heart attack, stroke, cancer, and
arthritis).

Figure 8: Traceplots for the parameters $\phi$ and $\kappa$ on the MEPS
dataset for five independent Markov chains.

We fit the quasi-power model to this data, with $\log\mu_{\theta}(x)=r(x)$ and
$V_{\eta}(\mu)=\mu^{\kappa}$, using the Laplace approximation to update the
tree topologies as described in Section 2.2 and the pseudo-likelihood approach
to update $(\phi,\kappa)$. For this analysis we ran five Markov chains in
parallel for a total of 2000 iterations per chain, with the first 1000
iterations of each chain discarded to burn-in. Traceplots for the parameters
$\phi$ and $\kappa$ are displayed in Figure 8, and we see that the chains
converged rapidly and there were no notable issues with the mixing. Bearing in
mind that our analysis is based on a different subpopulation than Natarajan et
al., (2008), our results are quite similar to theirs in terms of the plausible
values of $\kappa$.

Figure 9: Posterior summarization of variable importances (left) and posterior
variable inclusion probabilities (right). Variable importances correspond to
the average number of splits in the ensemble that make use of the listed
variable.

Having verified that the chain is mixing well and that the results for
$\kappa$ are broadly in agreement with the existing literature, we move onto
summarizing the posterior distribution. We begin by examining the _posterior
inclusion probabilities_ for each variable (Linero,, 2018) as well as their
variable importances as defined by Chipman et al., (2010) (defined to be the
posterior mean of the number of times a variable is used in a splitting rule).
Results are displayed in Figure 9, which shows roughly that the most important
indicators of future medical expenditures are first direct usage of the
healthcare system (number of prescription refills and number of dentist
trips), followed by indicators of poor health (cognitive limitations, cancer,
self-perceived health status etc). Variables that had lower posterior
inclusion probabilities include race, diabetes status, and depression status.

Next, we use the posterior projection approach of Woody et al., (2021) to
create summaries of the function $r(x)$, which can be useful both for
interpreting the model and for identifying interesting subpopulations. For
each sampled $r(x)$ we compute an optimal summary $\widetilde{r}(x)$ in a
given class of summary functions $\mathcal{H}$, defined by
$\widetilde{r}=\arg\min_{g\in\mathcal{H}}\sum_{i}\\{r(X_{i})-\widetilde{r}(X_{i})\\}^{2}$.
We use the _summary $R^{2}$_

$\displaystyle
R^{2}=1-\frac{\sum_{i}\\{r(X_{i})-\widetilde{r}(X_{i})\\}^{2}}{\sum_{i}\\{r(X_{i})-\bar{r}\\}^{2}}\qquad\text{where}\qquad\bar{r}=N^{-1}\sum_{i}r(X_{i})$

to measure the quality of the summary, with $R^{2}$ close to $1$ implying that
$\widetilde{r}(x)$ is a good approximation of $r(x)$ at the observed values of
the covariates.

In preliminary analysis, we found that the class of generalized additive
models (GAMs) that take $\widetilde{r}(x)=\sum_{j=1}^{P}h_{j}(x)$ provided
very strong summaries of $r(x)$, with the _partial effect_ $h_{j}(x)$ given by
a natural cubic spline (see Wood,, 2006, for a review) for the continuous
predictors. The summary $R^{2}$ for the GAM summaries were found to
concentrate between 95% and 98%.

Figure 10: Summary of the quasi-power BART model fit to the MEPS dataset.
Continuous variables display the posterior mean and confidence bands for the
associated term of the projection, while categorical variables display the
posterior median, 50% credible interval, and 95% credible interval using a
boxplot.

The posterior distributions of the $h_{j}$’s from the GAM summary of the BART
fit are given in Figure 10; log transformations were used prior to fitting the
GAM summary for the number of perscription refills (RX), income, and poverty
level. For the continuous variables, we see clear increasing trends in the
number of dentist visits, perceived health, and number of prescription
refills. For the categorical variables, the strongest effects are associated
to cancer diagnosis, cognitive limitations, and insurance status.

##  6 Discussion

In this paper, we propose to combine quasi-likelihood methods with Bayesian
nonparametric priors on the regression function $\mu_{\theta}(x)$ to obtain
methods that are both highly practical and robust to model misspecification.
The methods we propose depend only on correctly specifying a mean-variance
relation of the form $\operatorname{Var}_{\theta}(Y_{i}\mid
X_{i}=x)=\frac{\phi}{\omega}V\\{\mu_{\theta}(x)\\}$. We develop Bayesian
nonparametric versions of the quasi-binomial, quasi-Poisson, and quasi-gamma
models, and show that all of these approaches tare amenable to simple Bayesian
backfitting algorithms (Chipman et al.,, 2010) when $\mu_{\theta}(x)$ is
modeled using a suitable BART prior.

An obvious limitation of our approach is that we are still required to
correctly specify the variance function in order to obtain valid inferences.
The variance model can fail either because (i) the form of $V(\mu)$ is not
correctly specified or (ii) the variance depends on $X_{i}$ through features
other than $\mu_{\theta}(X_{i})$. We partially alleviate this by allowing for
the variance function itself to be estimated jointly with the dispersion
parameter $\phi$ using the Bayesian bootstrap quasi-likelihood (BBQ) update.
Both of these issues can be partially fixed by allowing the dispersion
parameter $\phi$ to depend on covariates $\phi(x)$ using, for example, the
_multiplicative decision tree_ model of Pratola et al., (2020). We defer such
extensions to future work.

This work is part of a larger program of attempting to perform Bayesian
inference while avoiding the specification of a data generating mechanism.
Quasi-likelihood uses only a specification of the mean function (which we
model nonparametrically) and a mean-variance function relationship. One might
hope further to either (i) eliminate the specification of the mean-variance
relationship or (ii) develop a method that is _robust_ to misspecification of
the mean-variance relation. A simple approach would be to fully commit to the
Bayesian bootstrap and _define_ $\theta$ as the minimum of a population-level
objective function; this works well for parametric problems, as we have
already seen here in our simulation experiments, but it is not clear how to
apply such methods in a non-parametric setting. Other approaches, such as
Bayesian exponentially tilted likelihood (BETEL) or Bayesian method of moments
(BMOM) are limited in scope due to computational challenges (Schennach,, 2005;
Siddhartha Chib and Simoni,, 2018; Yin,, 2009). A large part of the appeal of
quasi-likelihood methods for us is that they allow for the use of simple
Bayesian backfitting algorithms to fit BART-based models, which can no longer
be used with these other techniques. Future work will pursue computationally
efficient alternatives to BETEL to further reduce the modeling assumptions.

## References

  * Agnoletto et al., (2023) Agnoletto, D., Rigon, T., and Dunson, D. B. (2023).  Bayesian inference for generalized linear models via quasi-posteriors.
  * Arnold and Press, (1989) Arnold, B. C. and Press, S. J. (1989).  Compatible conditional distributions.  Journal of the American Statistical Association, 84(405):152–156.
  * Bissiri et al., (2016) Bissiri, P. G., Holmes, C. C., and Walker, S. G. (2016).  A General Framework for Updating Belief Distributions.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 78(5):1103–1130.
  * Casella, (2001) Casella, G. (2001).  Empirical Bayes Gibbs sampling.  Biostatistics, 2(4):485–500.
  * Chipman et al., (2010) Chipman, H. A., George, E. I., and McCulloch, R. E. (2010).  BART: Bayesian additive regression trees.  The Annals of Applied Statistics, 4(1):266 – 298.
  * Davidian and Carroll, (1988) Davidian, M. and Carroll, R. J. (1988).  A Note on Extended Quasi-Likelihood.  Journal of the Royal Statistical Society: Series B (Methodological), 50(1):74–82.
  * Dorie, (2023) Dorie, V. (2023).  dbarts: Discrete Bayesian additive regression trees sampler.  R package version 0.9-23.
  * Friedman, (1991) Friedman, J. H. (1991).  Multivariate adaptive regression splines.  The Annals of Statistics, 19(1):1–67.
  * Friedman, (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  The Annals of Statistics, 29(5):1189 – 1232.
  * Gelman and Speed, (1993) Gelman, A. and Speed, T. P. (1993).  Characterizing a Joint Probability Distribution by Conditionals.  Journal of the Royal Statistical Society: Series B (Methodological), 55(1):185–188.
  * Greco et al., (2008) Greco, L., Racugno, W., and Ventura, L. (2008).  Robust likelihood functions in Bayesian inference.  Journal of Statistical Planning and Inference, 138(5):1258–1270.
  * Hastie and Tibshirani, (2000) Hastie, T. and Tibshirani, R. (2000).  Bayesian backfitting (with comments and a rejoinder by the authors.  Statistical Science, 15(3):196 – 223.
  * Hill et al., (2020) Hill, J., Linero, A., and Murray, J. (2020).  Bayesian additive regression trees: A review and look forward.  Annual Review of Statistics and Its Application, 7(1):251–278.
  * Kapelner and Bleich, (2016) Kapelner, A. and Bleich, J. (2016).  bartMachine: Machine learning with Bayesian additive regression trees.  Journal of Statistical Software, 70(4):1–40.
  * Lin, (2006) Lin, L. (2006).  Quasi Bayesian likelihood.  Statistical Methodology, 3(4):444–455.
  * Linero, (2017) Linero, A. R. (2017).  A review of tree-based Bayesian methods.  Communications for Statistical Applications and Methods, 24(6):543–559.
  * Linero, (2018) Linero, A. R. (2018).  Bayesian regression trees for high-dimensional prediction and variable selection.  Journal of the American Statistical Association, 113(522):626–636.
  * Linero, (2022) Linero, A. R. (2022).  Generalized Bayesian additive regression trees models: Beyond conditional conjugacy.
  * Linero et al., (2020) Linero, A. R., Sinha, D., and Lipsitz, S. R. (2020).  Semiparametric mixed-scale models using shared bayesian forests.  Biometrics, 76(1):131–144.
  * Linero and Yang, (2018) Linero, A. R. and Yang, Y. (2018).  Bayesian regression tree ensembles that adapt to smoothness and sparsity.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(5):1087–1110.
  * Lyddon et al., (2019) Lyddon, S. P., Holmes, C. C., and Walker, S. G. (2019).  General Bayesian updating and the loss-likelihood bootstrap.  Biometrika, 106(2):465–478.
  * McCullaugh and Nelder, (1989) McCullaugh, P. and Nelder, J. A. (1989).  Generalized Linear Models.  Chapman & Hall/CRC, 2 edition.
  * Meyn and Tweedie, (2012) Meyn, S. and Tweedie, R. (2012).  Markov Chains and Stochastic Stability.  Communications and Control Engineering. Springer London.
  * Miller, (2021) Miller, J. W. (2021).  Asymptotic normality, concentration, and coverage of generalized posteriors.  Journal of Machine Learning Research, 22(168):1–53.
  * Murray, (2021) Murray, J. S. (2021).  Log-linear Bayesian additive regression trees for multinomial logistic and count regression models.  Journal of the American Statistical Association, 116(534):756–769.
  * Natarajan et al., (2008) Natarajan, S., Lipsitz, S. R., Fitzmaurice, G., Moore, C. G., and Gonin, R. (2008).  Variance estimation in complex survey sampling for generalized linear models.  Journal of the Royal Statistical Society Series C: Applied Statistics, 57(1):75–87.
  * Neal, (2003) Neal, R. M. (2003).  Slice sampling.  The Annals of Statistics, 31(3):705 – 767.
  * Neal, (2011) Neal, R. M. (2011).  MCMC using Hamiltonian dynamics.  In Books, S., Gelman, A., Jones, G. L., and Meng, X. L., editors, Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC.
  * Nelder and Pregibon, (1987) Nelder, J. A. and Pregibon, D. (1987).  An extended quasi-likelihood function.  Biometrika, 74(2):221–232.
  * Nelder and Wedderburn, (1972) Nelder, J. A. and Wedderburn, R. W. M. (1972).  Generalized Linear Models.  Royal Statistical Society. Journal. Series A: General, 135(3):370–384.
  * Pratola et al., (2020) Pratola, M. T., Chipman, H. A., George, E. I., and McCulloch, R. E. (2020).  Heteroscedastic BART via multiplicative regression trees.  Journal of Computational and Graphical Statistics, 29(2):405–417.
  * Rasmussen and Williams, (2005) Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Ročková and van der Pas, (2020) Ročková, V. and van der Pas, S. (2020).  Posterior concentration for Bayesian regression trees and forests.  The Annals of Statistics, 48(4):2108–2131.
  * Schennach, (2005) Schennach, S. M. (2005).  Bayesian exponentially tilted empirical likelihood.  Biometrika, 92(1):31–46.
  * Siddhartha Chib and Simoni, (2018) Siddhartha Chib, M. S. and Simoni, A. (2018).  Bayesian estimation and comparison of moment condition models.  Journal of the American Statistical Association, 113(524):1656–1668.
  * Smyth and Verbyla, (1999) Smyth, G. K. and Verbyla, A. P. (1999).  Adjusted likelihood methods for modelling dispersion in generalized linear models.  Environmetrics: The official journal of the International Environmetrics Society, 10(6):695–709.
  * Sparapani et al., (2021) Sparapani, R., Spanbauer, C., and McCulloch, R. (2021).  Nonparametric machine learning and efficient computation with Bayesian additive regression trees: The BART R package.  Journal of Statistical Software, 97(1):1–66.
  * Ventura et al., (2010) Ventura, L., , Cabras, S., and Racugno, W. (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 Racugno, W. (2016).  Pseudo-likelihoods for bayesian inference.  In Di Battista, T., Moreno, E., and Racugno, W., editors, Topics on Methodological and Applied Statistical Inference, pages 205–220, Cham. Springer International Publishing.
  * Wedderburn, (1974) Wedderburn, R. W. M. (1974).  Quasi-likelihood functions, generalized linear models, and the Gauss—Newton method.  Biometrika, 61(3):439–447.
  * Wood, (2006) Wood, S. N. (2006).  Generalized Additive Models: An Introduction with R.  Chapman and Hall/CRC, second edition.
  * Woody et al., (2021) Woody, S., Carvalho, C. M., and Murray, J. S. (2021).  Model interpretation through lower-dimensional posterior summarization.  Journal of Computational and Graphical Statistics, 30(1):144–161.
  * Yin, (2009) Yin, G. (2009).  Bayesian generalized method of moments.  Bayesian Analysis, 4(2):191 – 207.

##  Appendix A Gibbs Sampler for the Quasi-Poisson Model

Let $\mathcal{L}(\mathcal{T}_{t})$ denote the set of leaves of tree
$\mathcal{T}_{t}$ and let $\lambda_{\ell}$ denote the leaf node parameter
associated to leaf $\ell$. Further, let $\mathcal{T}_{-t}$ denote the set of
trees excluding tree $t$, let $\mathcal{M}_{-t}$ denote the set of leaf node
parameters for all trees excluding tree $t$, and let $\zeta_{i}=\sum_{m\neq
t}g(X_{i};\mathcal{T}_{m},\mathcal{M}_{m})$ denote the sum of all trees except
for the $t$th tree. Following Chipman et al., (2010), we use a two-stage
sampling algorithm to update the tree parameters
$(\mathcal{T}_{1},\ldots,\mathcal{T}_{T})$. We first compute the integrated
likelihood of tree $\mathcal{T}_{t}$ given $\mathcal{T}_{-t}$ and
$\mathcal{M}_{-t}$ as

$\displaystyle
L(\mathcal{T}_{t})=\int\prod_{\ell\in\mathcal{L}(\mathcal{T}_{t})}\prod_{i:X_{i}\leadsto\ell}q_{\phi}(Y_{i}\mid\zeta_{i}+\lambda_{\ell})\
d\lambda_{\ell},$

where $X_{i}\leadsto\ell$ denotes the event that $X_{i}$ is associated to leaf
$\ell$ of tree $t$ and

$\displaystyle
q_{\phi}(y\mid\zeta)=\exp\left(-\frac{D\\{y,g^{-1}(\zeta)\\}}{2\phi}\right),$

is the quasi-likelihood of $y$ given $\zeta$.

For the quasi-Poisson model with $\mu_{\theta}(x)=\exp\\{r(x)\\}$ we have
$-D(y,\mu)/2=c(y)+y\log\mu-e^{-\mu}$ for some function $c(y)$. The likelihood
associated to $(\mathcal{M}_{t},\mathcal{T}_{t})$ is therefore proportional to

$\displaystyle\prod_{\ell\in\mathcal{L}(\mathcal{T}_{t})}\prod_{i:X_{i}\leadsto\ell}\exp\left(\frac{Y_{i}\zeta_{i}+Y_{i}\lambda_{\ell}-e^{\zeta_{i}}e^{\lambda_{\ell}}}{\phi}\right)\propto\prod_{\ell\in\mathcal{L}(\mathcal{T})}\exp\left(A_{\ell}\lambda_{\ell}-B_{\ell}e^{\lambda_{\ell}}\right),$

where $A_{\ell}=\sum_{i:X_{i}\leadsto\ell}\frac{Y_{i}}{\phi}$ and
$B_{\ell}=\sum_{i:X_{i}\leadsto\ell}\frac{e^{\zeta_{i}}}{\phi}$. Assuming that
$\lambda_{\ell}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\log\operatorname{Gam}(a,b)$,
we can compute the integrated likelihood of tree $\mathcal{T}_{t}$ from this
as

$\displaystyle
L(\mathcal{T})=\prod_{\ell\in\mathcal{L}(\mathcal{T})}\int_{-\infty}^{\infty}\frac{b^{a}}{\Gamma(a)}\exp(a\lambda-
be^{\lambda})\times\exp(A_{\ell}\lambda-B_{\ell}e^{\lambda})\
d\lambda=\prod_{\ell\in\mathcal{L}(\mathcal{T})}\frac{b^{a}}{\Gamma(a)}\times\frac{\Gamma(A_{\ell}+a)}{(B_{\ell}+b)^{A_{\ell}+a}}.$

Additionally, we can see that the full-conditional of $\lambda_{\ell}$ under
this model is given by
$\lambda_{\ell}\sim\log\operatorname{Gam}(A_{\ell}+a,b+B_{\ell})$. This leads
to the following two-step sampling algorithm for updating
$(\mathcal{T}_{t},\mathcal{M}_{t})$:

  1. 1.

Propose a tree $\mathcal{T}^{\prime}\sim
Q(\mathcal{T}_{t}\to\mathcal{T}^{\prime})$ using, for example, the Birth,
Depth, or Perturb moves described by Kapelner and Bleich, (2016) and its
appendix.

  2. 2.

Set $\mathcal{T}_{t}=\mathcal{T}^{\prime}$ if
$U\leq\frac{L(\mathcal{T}^{\prime})\,\pi(\mathcal{T}^{\prime})\,Q(\mathcal{T}^{\prime}\to\mathcal{T}_{t})}{L(\mathcal{T}_{t})\,\pi(\mathcal{T}_{t})\,Q(\mathcal{T}_{t}\to\mathcal{T}^{\prime})}$,
where $U\sim\operatorname{Uniform}(0,1)$; otherwise, leave $\mathcal{T}_{t}$
unmodified.

  3. 3.

For each $\ell\in\mathcal{L}(\mathcal{T}_{t})$, sample
$\lambda_{\ell}\sim\operatorname{Gam}(a+A_{\ell},b+B_{\ell})$.

##  Appendix B Gibbs Sampler for Quasi-Gamma

For the quasi-gamma model we have

$\displaystyle-\frac{D(y,\mu)}{2}=c(y)-\frac{y}{\mu}-\log(\mu).$

We then specify the model $\mu_{\theta}(x)=e^{-r(x)}$. Defining $\zeta_{i}$
and the other tree-based quantities as we did for the quasi-Poisson model the
integrated likelihood of $\mathcal{T}_{t}$ given
$(\mathcal{T}_{-t},\mathcal{M}_{-t})$ is therefore given by

$\displaystyle L(\mathcal{T}_{t})=$
$\displaystyle\prod_{\ell\in\mathcal{L}(\mathcal{T}_{t})}\int\prod_{i:X_{i}\leadsto\ell}\exp\left(\frac{\lambda-
Y_{i}\,e^{\zeta_{i}}e^{\lambda}}{\phi}\right)\frac{b^{a}}{\Gamma(a)}\exp(a\lambda-
be^{\lambda})\ d\lambda$
$\displaystyle\qquad\propto\prod_{\ell\in\mathcal{L}(\mathcal{T}_{t})}\int\exp\left(A_{\ell}\lambda-
B_{\ell}e^{\lambda}\right)\frac{b^{a}}{\Gamma(a)}\exp(a\lambda-be^{\lambda})\
d\lambda$
$\displaystyle\qquad=\prod_{\ell\in\mathcal{L}(\mathcal{T}_{t})}\frac{b^{a}}{\Gamma(a)}\times\frac{\Gamma(a+A_{\ell})}{(b+B_{\ell})^{a+A_{\ell}}},$

where $A_{\ell}=N_{\ell}/\phi$ where $N_{\ell}$ is the number of observations
assigned to leaf node $\ell$ of tree $\mathcal{T}_{t}$, and
$B_{\ell}=\sum_{i:X_{i}\leadsto\ell}Y_{i}\,e^{\zeta_{i}}/\phi$. Similarly, we
also have the full conditional
$\lambda_{\ell}\sim\log\operatorname{Gam}(a+A_{\ell},b+B_{\ell})$. The
algorithm is now identical to that for the quasi-Poisson model.

##  Appendix C Gibbs Sampler for the Quasi-Power Model

The quasi-power model uses the variance function $V(\mu)=\mu^{\kappa}$ where
$\kappa$ is estimated based on the data. By routine calculations, the quasi-
likelihood associated to the quasi-power model is given by

$\displaystyle
q_{\phi}(y\mid\mu,\kappa)=\exp\left\\{\frac{1}{\phi}\left(\frac{y\mu^{1-\kappa}}{1-\kappa}-\frac{\mu^{2-\kappa}}{2-\kappa}\right)\right\\},$

for $\kappa\neq 1,2$, with the special cases of $\kappa=1,2$ corresponding to
the quasi-Poisson model and the quasi-gamma model, respectively. We consider
the exponential link $\mu_{\theta}(x)=\exp\\{r(x)\\}$.

As usual, we let $\zeta_{i}=\sum_{m\neq
t}g(X_{i};\mathcal{T}_{m},\mathcal{M}_{m})$. The integrated quasi-likelihood
associated $\mathcal{T}_{t}$ given
$(\kappa,\phi,\mathcal{D},\mathcal{T}_{-t},\mathcal{M}_{-t})$ under this model
is given by

$\displaystyle\prod_{\ell\in\mathcal{T}_{t}}\int\prod_{i:X_{i}\leadsto\ell}\exp\left\\{\frac{1}{\phi}\left(\frac{Y_{i}e^{\zeta_{i}(1-\kappa)}e^{\lambda(1-\kappa)}}{1-\kappa}-\frac{e^{\zeta_{i}(2-\kappa)}e^{\lambda(2-\kappa)}}{2-\kappa}\right)\right\\}\,\operatorname{Normal}(\lambda\mid
0,\sigma^{2}_{\lambda})\ d\lambda$
$\displaystyle=\prod_{\ell\in\mathcal{T}_{t}}\int\exp\left\\{\frac{1}{\phi}\left(\frac{A_{\ell}e^{\lambda(1-\kappa)}}{1-\kappa}-\frac{B_{\ell}e^{\lambda(2-\kappa)}}{2-\kappa}\right)\right\\}\operatorname{Normal}(\lambda\mid
0,\sigma^{2}_{\lambda})\ d\lambda$

where $A_{\ell}=\sum_{i:X_{i}\leadsto\ell}Y_{i}e^{\zeta_{i}(1-\kappa)}$ and
$B_{\ell}=\sum_{i:X_{i}\leadsto\ell}e^{\zeta_{i}(2-\kappa)}$.

Next, we define
$\Lambda_{\ell}(\lambda)=\frac{1}{\phi}\left(\frac{A_{\ell}e^{\lambda(1-\kappa)}}{1-\kappa}-\frac{B_{\ell}e^{\lambda(2-\kappa)}}{2-\kappa}\right)$.
Directly optimizing $\Lambda_{\ell}(\lambda)$ shows that it is maximized at
$\widehat{\lambda}_{\ell}=\log(A_{\ell}/B_{\ell})$. Proceeding by the Laplace
approximation, we expand about $\widehat{\lambda}$ to second order to get the
approximation

$\displaystyle\begin{split}L(\mathcal{T}_{t})&\approx\prod_{\ell\in\mathcal{T}_{t}}e^{\Lambda_{\ell}(\widehat{\lambda}_{\ell})}\int\exp\left\\{-\frac{(\lambda-\widehat{\lambda}_{\ell})^{2}I(\widehat{\lambda})}{2}\right\\}\operatorname{Normal}(\lambda\mid
0,\sigma^{2}_{\lambda})\ d\lambda\\\
&=\prod_{\ell\in\mathcal{T}_{t}}e^{\Lambda_{\ell}(\widehat{\lambda}_{\ell})}\sqrt{\frac{I^{-1}(\widehat{\lambda}_{\ell})}{\sigma^{2}_{\lambda}+I^{-1}(\widehat{\lambda}_{\ell})}}\exp\left\\{-\frac{1}{2}\frac{\widehat{\lambda}^{2}}{(I^{-1}(\widehat{\lambda}_{\ell})+\sigma^{2}_{\lambda})}\right\\},\end{split}$
(6)

where
$I(\widehat{\lambda}_{\ell})=\left.\frac{\partial^{2}}{\partial\lambda^{2}}\Lambda_{\ell}(\lambda)\right|_{\lambda=\widehat{\lambda}_{\ell}}$.
We use (6) to compute an approximate acceptance probability for the Bayesian
backfitting algorithm.

After updating $\mathcal{T}_{t}$, we can sample the $\lambda_{\ell}$’s from
their full conditionals, which are proportional to

$\displaystyle\exp\left\\{\frac{1}{\phi}\left(\frac{A_{\ell}e^{\lambda_{\ell}(1-\kappa)}}{1-\kappa}-\frac{B_{\ell}e^{\lambda_{\ell}(2-\kappa)}}{2-\kappa}\right)\right\\}\operatorname{Normal}(\lambda_{\ell}\mid
0,\sigma^{2}_{\lambda}).$

Applying the Laplace approximation again, this full conditional is
approximately proportional to

$\displaystyle\exp\left\\{-\frac{(\lambda_{\ell}-\widehat{\lambda}_{\ell})^{2}I(\widehat{\lambda}_{\ell})}{2}\right\\}\operatorname{Normal}(\lambda_{\ell}\mid
0,\sigma^{2}_{\lambda}),$

which is the kernel of a normal distribution with mean and variance given by

$\displaystyle
m=\frac{I(\widehat{\lambda}_{\ell})\widehat{\lambda}_{\ell}}{I(\widehat{\lambda}_{\ell})+\sigma^{-2}_{\lambda}}\qquad\text{and}\qquad
v=\frac{1}{I(\widehat{\lambda}_{\ell})+\sigma^{-2}_{\lambda}}.$

This approximation can be further refined by applying, for example, several
rounds of slice sampling (Neal,, 2003) to the exact full conditional after
sampling $\lambda_{\ell}$.

As described above, the Gibbs sampler proposed for the quasi-power model is
approximate, and does not hold the exact quasi-posterior invariant (even
ignoring the updates for $\phi$ and $\kappa$). Linero, (2022) gives a general
recipe for converting Laplace approximations such as these into proposals that
leave the target posterior exactly invariant, and the modifications can be
applied if this is desired. In practice, we find that the Laplace
approximation is adequate on its own for identifying useful tree topologies,
and that the effect of any inaccuracy in the Laplace approximation for
updating $\lambda_{\ell}$ can be more-or-less eliminated by following up the
sampling with a single slice sampling update. Moreover, because of the large
amount of regularization applied by the prior — the prior takes
$\sigma^{2}_{\lambda}$ very small, especially when many trees are used — the
Laplace approximation is typically highly accurate without making any
corrections, even if $\Lambda_{\ell}(\lambda)$ is not close to quadratic.

##  Appendix D Gibbs Sampler for the Quasi-Multinomial

Let $Z_{i}=n_{i}\,Y_{i}$, where $Y_{i}$ is modeled with a quasi-multinomial
model. Then the quasi-likelihood is given by

$\displaystyle
Q(\theta)=\prod_{i}\frac{\exp\\{Z_{i}^{\top}r(X_{i})/\phi\\}}{\left[\sum_{k}\exp\\{r_{k}(X_{i})\\}\right]^{n_{i}/\phi}}.$

We now introduce latent variables
$\xi_{i}\sim\operatorname{Gam}(n_{i}/\phi,\sum_{k}\exp\left\\{r_{k}(X_{i})\right\\})$,
which leads to the joint likelihood

$\displaystyle Q(\theta,\xi)$
$\displaystyle\propto\prod_{i}\exp\left\\{\frac{Z_{i}^{\top}r(X_{i})}{\phi}-\xi_{i}\sum_{k}e^{r_{k}(X_{i})}\right\\}\xi_{i}^{n_{i}/\phi-1}.$

Defining $\zeta_{ik}=\sum_{m\neq
t}[g(X_{i};\mathcal{T}_{t},\mathcal{M}_{t})]_{k}$, the integrated likelihood
associated to a given tree $\mathcal{T}_{t}$ (conditional on the $\xi$’s as
well as $\mathcal{T}_{-t},\mathcal{M}_{-t}$, and $\phi$) is then proportional
to

$\displaystyle\prod_{\ell}\prod_{k}\int\prod_{i:X_{i}\leadsto\ell}\exp\left\\{\frac{Z_{ik}\,(\zeta_{ik}+\lambda)}{\phi}-\xi_{i}e^{\zeta_{ik}}e^{\lambda}\right\\}\,\pi(\lambda)\
d\lambda$
$\displaystyle\propto\prod_{\ell}\prod_{k}\int\exp\left(\frac{A_{\ell
k}\lambda}{\phi}-B_{\ell k}e^{\lambda}\right)\,\pi(\lambda)\ d\lambda,$

where $A_{\ell k}=\sum_{i:X_{i}\leadsto\ell}Z_{ik}$ and $B_{\ell
k}=\sum_{i:X_{i}\leadsto\ell}\xi_{i}e^{\zeta_{ik}}$. Assuming the
$\log\operatorname{Gam}(a,b)$ prior for the $\lambda_{\ell k}$’s, this gives

$\displaystyle
L(\mathcal{T}_{t})\propto\prod_{\ell}\prod_{k}\frac{b^{a}}{\Gamma(a)}\times\frac{\Gamma(a+A_{\ell
k}/\phi)}{(b+B_{\ell k})^{a+A_{\ell k}/\phi}},$

and similarly to previous calculations the full conditional of $\lambda_{\ell
k}$ is a $\log\operatorname{Gam}(a+A_{\ell k}/\phi,b+B_{\ell k})$
distribution.

##  Appendix E Assumptions

This section provides the regularity conditions that are required to establish
the theoretical results in Section 4. First, we make some assumptions about
the form of the quasi-likelihood $Q(\theta)$.

#### Condition P

The function $r(x)=x^{\top}\beta$ is related to $\mu_{\theta}(x)$ the relation
$g\\{\mu_{\theta}(x)\\}=r(x)$ where $g(\mu)$ is a continuous strictly-
increasing function. Additionally, the function $V_{\eta}\\{g^{-1}(\psi)\\}$
is continuous, strictly positive, and
$D\\{Y_{i},g^{-1}(\psi)\\}=2\int_{g^{-1}(\psi)}^{y}\frac{y-t}{V_{\eta}(t)}\
dt$ exists for all $\psi\in(-\infty,\infty)$.

Next, the following assumption is used to establish the existence of the
stationary distribution.

#### Condition A

The prior $\pi(\beta)$ is continuous and bounded away from $0$ on its support,
which is assumed to be a ball $\\{\beta:\|\beta\|\leq R\\}$, and for the PL
approach $g(\phi\mid\beta)$ is an inverse gamma distribution
$\text{InvGam}(N/2,\widehat{\phi}(\beta)N/2)$ truncated to an interval $[a,b]$
with $a>0$ and $b<\infty$. Additionally, the data $\mathcal{D}$ is such that
$Q(\theta)$ is well-defined for all $\theta$ with $\|\beta\|\leq R$ and
$\phi\in[a,b]$.

The following assumption states, essentially, that the model is correctly
specified.

#### Condition F

The true density of $[Y_{i}\mid X_{i},\omega_{i}]$, denoted $f_{0}(y\mid
x,\omega)$, satisfies the quasi-likelihood moment conditions for some
$\theta_{0}=(\phi_{0},\beta_{0})$ such that $\phi_{0}>0$,
$r_{0}(x)=x^{\top}\beta_{0}$ for some $\beta_{0}\in\mathbb{R}^{P}$, and link
function $g\\{\mu_{0}(x)\\}=r_{0}(x)$, with the $Y_{i}$’s independent .
Additionally, the prior $\pi(\beta)$ is continuous at $\beta_{0}$ with
$\pi(\beta_{0})>0$.

Regularity conditions on the structure of the quasi-likelihood and data
generating mechanism are also required. The following conditions are borrowed
from Appendix A.3 of Agnoletto et al., (2023) and guarantee a pointwise
Bernstein-von Mises result for fixed $\phi$, modified to be slightly simpler
to take advantage of the boundedness assumptions in Condition A and Condition
F. They hold very generally.

#### Condition D

Let $X_{1},X_{2},\ldots$ be the sequence of covariate values and $[Y_{i}\mid
X_{i},\omega_{i}]$ be independent and distributed according to Condition F.
The following conditions hold:

  * D1

We have
$\lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^{n}\omega_{i}\mathbb{E}_{\theta_{0}}\left(D\\{Y_{i};g^{-1}(X_{i}^{\top}\beta)\\}\mid
X_{i}=x\right)=D_{\infty}(\beta)$ for some $D_{\infty}(\beta)$. Furthermore,

$\displaystyle\lim_{n\to\infty}\sum_{i=1}^{n}\frac{\omega_{i}^{2}\operatorname{Var}\left(D\\{Y_{i};g^{-1}(X_{i}^{\top}\beta)\\}\right)}{n^{2}}<\infty.$

  * D2

$D(Y_{i};\mu_{i})$ is three-times continuously differentiable on the support
of $\pi(\beta)$, with the third derivative uniformly bounded.

  * D3

The matrix $\bm{X}^{\top}W\bm{X}/N$ converges to a positive definite matrix
$H(\beta_{0})$, where
$W=\operatorname{diag}\left(\frac{\omega_{i}}{g^{\prime}(\mu_{i})^{2}\,V_{\eta}(\mu_{i})}\right)$
and $\mu_{i}\equiv g^{-1}(r_{0}(X_{i}))$.

  * D4

$D_{\infty}(\beta)>D_{\infty}(\beta_{0})$ if $\beta\neq\beta_{0}$ for all
$\beta$ in the support of $\pi(\beta)$.

Finally, we assume that the moment estimator is sufficiently regular in the
sense that $\widehat{\phi}(\beta_{n})\to\phi_{0}$ provided that
$\beta_{n}\to\beta_{0}$ sufficiently fast.

#### Condition G

Let $\widehat{\phi}(\beta)$ denote the moment estimator
$\frac{1}{n}\sum_{i}\frac{\omega_{i}\\{Y_{i}-g^{-1}(X_{i}^{\top}\beta)\\}^{2}}{V_{\eta}(g^{-1}(X_{i}^{\top}\beta))}$,
and let
$\widehat{\beta}=\arg\min_{\beta}\sum_{i}\omega_{i}D\\{Y_{i};g^{-1}(X_{i}^{\top}\beta)\\}$
be the maximum quasi-likelihood estimator of $\beta$. Then, with probability
$1$, $\widehat{\phi}(\widehat{\beta}+h/\sqrt{n})\to\phi_{0}$ for all $h$.

The above conditions also imply correct Frequentist properties of the quasi-
posterior, as under standard regularity conditions (and assuming that the mean
and variance relations are correctly specified) it is well-known that
$\widehat{\beta}\stackrel{{\scriptstyle{}_{\bullet}}}{{\sim}}N(\beta_{0},\phi
N^{-1}H(\beta_{0})^{-1})$.

##  Appendix F Theoretical Results

###  F.1 Existence of the Stationary Distribution

We note that Condition A, which truncates the support of $\pi(\theta)$ to a
compact set, is overly strong and serves to eliminate technical issues in
verifying that the chain is irreducible, aperiodic, and positive recurrent,
which are standard sufficient conditions for the existence of a stationary
distribution (see Meyn and Tweedie,, 2012, Chapter 10).

###### Proof.

Let $\mathcal{S}=\\{(\beta,\phi):\|\beta\|\leq R,\phi\in[a,b]\\}$. Note that
the Markov transition function (MTF) of the two-step Gibbs sampler is given by

$\displaystyle
R\\{(\beta,\phi)\to(\beta^{\prime},\phi^{\prime})\\}=\widetilde{\pi}(\beta^{\prime}\mid\phi)\,g(\phi^{\prime}\mid\beta^{\prime}).$

From Condition A and Condition P it follows that
$\widetilde{\pi}(\beta^{\prime}\mid\phi)\,g(\phi^{\prime}\mid\beta^{\prime})$
is bounded away from $0$ on the set $\mathcal{S}\times\mathcal{S}$. We now
apply Theorem 1.3.1 (ii) of Meyn and Tweedie, (2012) to conclude the existence
of a unique stationary distribution:

  * •

The set $\mathcal{S}$ is small because the MTF is bounded away from $0$ on
$\mathcal{S}$, as this implies that the MTF is minorized by a uniform measure
on $\mathcal{S}$.

  * •

The hitting time for $\mathcal{S}$ is trivially $1$ regardless of how the
chain is initialized.

Hence Theorem 1.3.1 (i) and (ii) imply the result.

Next, we characterize the stationary distribution
$\widetilde{\pi}(\beta,\phi)$. If we sample
$(\beta,\phi)\sim\widetilde{\pi}(\beta,\phi)$ and then sample the transition
$R\\{(\beta,\phi)\to(\beta^{\prime},\phi^{\prime})\\}$ then, by definition,
$(\beta^{\prime},\phi^{\prime})\sim\widetilde{\pi}(\beta,\phi)$, but from the
structure of the Gibbs sampler it is clear that the conditional associated to
$\phi^{\prime}$ is $g(\phi^{\prime}\mid\beta^{\prime})$; hence,
$\widetilde{\pi}(\beta,\phi)=\widetilde{\pi}(\beta)\,g(\phi\mid\beta)$ for
some $\widetilde{\pi}(\beta)$.

To characterize $\widetilde{\pi}(\beta)$, we use the definition of the
stationary distribution to obtain

$\displaystyle\widetilde{\pi}(\beta^{\prime})\,g(\phi^{\prime}\mid\beta^{\prime})$
$\displaystyle=\int\widetilde{\pi}(\beta,\phi)\widetilde{\pi}(\beta^{\prime}\mid\phi)\,g(\phi^{\prime}\mid\beta^{\prime})\
d\beta d\phi$
$\displaystyle=g(\phi^{\prime}\mid\beta^{\prime})\int\widetilde{\pi}(\beta,\phi)\widetilde{\pi}(\beta^{\prime}\mid\phi)\
d\beta\ d\phi$ $\displaystyle=g(\phi^{\prime}\mid\beta^{\prime})\int
m(\phi)\,\widetilde{\pi}(\beta^{\prime}\mid\phi)\ d\phi,$

where $m(\phi)$ is the marginal distribution of $\phi$ for the stationary
distribution. Hence
$\widetilde{\pi}(\beta)=\int\widetilde{\pi}(\beta\mid\phi)\,m(\phi)\ d\phi$.

To obtain the result $m(\phi)=\int g(\phi\mid\beta)\,\widetilde{\pi}(\beta)\
d\beta$, we repeat the same logic with the roles of $\phi$ and $\beta$
reversed, and analyze the _dual chain_ that reverses the order of updating
$\beta$ and $\phi$; note that, because our sampler is not a genuine Gibbs
sampler, the dual chain does not have the same stationary distribution as the
original chain, however they do have the same marginal distributions, which is
sufficient. ∎

###  F.2 Bernstein-von Mises Theorem

We begin by proving a _uniform_ Bernstein-von Mises theorem, establishing
uniform convergence of the posterior $\widetilde{\pi}(\beta\mid\phi)$.
Throughout, we use the symbol $A\lesssim B$ to mean that $A\leq CB$ for some
positive constant $C>0$ that is allowed to depend only on
$(P,\beta_{0},a,b,H,\pi)$.

######  Lemma 1.

Suppose that Condition P, Condition A, Condition F, and D1–D4 hold. Then

$\displaystyle\sup_{\phi\in[a,b]}\int|q_{n}^{\phi}-N^{\phi}|\ dh\to
0\qquad\text{almost surely}$

where $q_{n}^{\phi}(h)$ is the quasi-posterior distribution of
$\sqrt{n}(\beta-\widehat{\beta}_{n})$, $\widehat{\beta}_{n}$ is the maximum
quasi-likelihood estimator of $\beta$, and $N^{\phi}$ is the density of a
normal distribution with mean $0$ and variance $\phi\,H^{-1}(\beta_{0})$.

###### Proof.

First, Theorem 3 of Agnoletto et al., (2023) implies the non-uniform version
of this statement: for all $\phi\in[a,b]$, we have
$\int|q_{n}^{\phi}-N^{\phi}|\ dh\to 0$ almost surely (the statement of the
result is in-probability, but as a corollary of results of Miller, (2021) the
techniques used prove almost-sure convergence). Now, fix $\delta>0$ and let
$a=\lambda_{1}<\lambda_{2}<\cdots\lambda_{K}=b$ be an evenly-space
$\delta$-cover of $[b^{-1},a^{-1}]$, i.e.,
$\lambda_{k+1}-\lambda_{k}\leq\delta$; let $\mathcal{K}$ denote this
$\delta$-cover.

For fixed $\phi$, let $\lambda_{k}=\lambda$ be such that
$\lambda_{k}\leq\phi^{-1}\leq\lambda_{k+1}$. Then

$\displaystyle\int|q^{\phi}_{n}-N^{\phi}|\
dh\leq\underbrace{\int|q_{n}^{\phi}-q_{n}^{\lambda^{-1}}|\
dh}_{(i)}+\underbrace{\int|q_{n}^{\lambda^{-1}}-N^{\lambda^{-1}}|\
dh}_{(ii)}+\underbrace{\int|N^{\lambda^{-1}}-N^{\phi}|\ dh}_{(iii)}.$

We now bound (i) – (iii) uniformly in $\phi$.

Second term: From the pointwise version of the result we have that
$\max_{\lambda\in\mathcal{K}}\int|q_{n}^{\lambda^{-1}}-N^{\lambda^{-1}}|\to 0$
as $n\to\infty$ almost surely; hence (ii) converges almost surely $0$ in
probability uniformly in $\phi$.

Third term: By Pinsker’s inequality, the third term is bounded by a multiple
of

$\displaystyle\sqrt{\log\lambda^{-1}/\phi+\phi/\lambda^{-1}-1}.$

Letting $\Delta=\frac{\lambda^{-1}}{\phi}-1$ and bounding
$\log(1+\Delta)\leq\Delta$ we have the bound

$\displaystyle\sqrt{\log(1+\Delta)-\Delta/(1+\Delta)}\leq\frac{\Delta}{\sqrt{1+\Delta}}.$

Now, $\Delta=\frac{\phi^{-1}-\lambda}{\lambda}\leq b\delta$, while
$1+\Delta=\phi^{-1}/\lambda\geq a/b$ so that we can bound this last term by
$\delta\times\sqrt{b/a}$. That is,

$\displaystyle\sup_{\phi}\int|N^{\lambda^{-1}}-N^{\phi}|\ dh\lesssim\delta.$

First term: We first bound the total variation distance with the Hellinger
distance:

$\displaystyle\int|q_{n}^{\lambda^{-1}}-q_{n}^{\phi}|\ dh\lesssim
1-\int\sqrt{q_{n}^{\phi}\,q_{n}^{\lambda^{-1}}}\
dh=1-\sqrt{\frac{[Z_{n}^{(\rho+\lambda)/2}]^{2}}{Z_{n}^{\rho}\,Z_{n}^{\lambda}}}$
(7)

where $\rho\equiv\phi^{-1}$,

$\displaystyle Z_{n}^{\rho}$ $\displaystyle=\int\exp\left\\{\log
Q_{\phi}(\widehat{\beta}+h/\sqrt{n})-\log
Q_{\phi}(\widehat{\beta})\right\\}\,\pi(\widehat{\beta}+h/\sqrt{n})\ dh$
$\displaystyle=\int\exp\left\\{\rho[\log
Q_{1}(\widehat{\beta}+h/\sqrt{n})-\log
Q_{1}(\widehat{\beta})]\right\\}\pi(\widehat{\beta}+h/\sqrt{n})\ dh$

and the equality in (7) (after routine calculation) comes from the fact that

$\displaystyle q_{n}^{\phi}(h)$ $\displaystyle\propto
Q_{\phi}(\widehat{\beta}+h/\sqrt{n})\,\pi(\widehat{\beta}+h/\sqrt{n})\qquad\text{where}$
$\displaystyle Q_{\phi}(\beta)$
$\displaystyle=\exp\left\\{-\frac{\sum_{i=1}^{n}\omega_{i}D\\{Y_{i};g^{-1}(X_{i}^{\top}\beta)\\}}{2\phi}\right\\}$

is the quasi-likelihood. Next, we observe that for all
$\lambda_{k}\leq\rho\leq\lambda_{k+1}$ we have
$(Z_{n}^{\lambda_{k}})^{\rho/\lambda_{k}}\leq
Z_{n}^{\rho}\leq(Z_{n}^{\lambda_{k+1}})^{\rho/\lambda_{k+1}}$ by Jensen’s
inequality. Accordingly, the Hellinger distance is bounded above by

$\displaystyle
1-\frac{(Z_{n}^{\lambda_{k+1}})^{(\lambda_{k}+\rho)/2{\lambda_{k+1}}}}{(Z_{n}^{\lambda_{k}})^{(\lambda_{k}+\rho)/2\lambda_{k}}}=1-\left(\frac{Z_{n}^{\lambda_{k+1}}}{Z_{n}^{\lambda_{k}}}\right)^{(\lambda+\rho)/(2\lambda_{k+1})}\times\left(Z_{n}^{\lambda_{k}}\right)^{\frac{\lambda+\rho}{2}\times\frac{\lambda_{k}-\lambda_{k+1}}{\lambda_{k}}}.$

Next, we note that

$\displaystyle\frac{\lambda_{k}+\rho}{2\lambda_{k+1}}\log\left(\frac{Z_{n}^{\lambda_{k+1}}}{Z_{n}{{}^{\lambda}_{k}}}\right)$
$\displaystyle\geq-\left|\log\left(\frac{Z_{n}^{\lambda_{k+1}}}{Z_{n}^{\lambda_{k}}}\right)\right|\equiv
A_{nk}\qquad\text{and}\qquad$
$\displaystyle-\frac{\lambda_{k}+\rho}{2}\times\frac{\lambda_{k+1}-\lambda_{k}}{\lambda_{k}}\log
Z_{n}^{\lambda_{k}}$
$\displaystyle\geq-\frac{\delta\lambda_{k+1}}{\lambda_{k}}|\log
Z_{n}^{\lambda_{k}}|\equiv B_{nk}.$

As these bounds do not depend on $\phi$, we have

$\displaystyle\sup_{\phi}\int|q_{n}^{\phi}-q_{n}^{\lambda^{-1}}|\lesssim\max_{k}1-e^{A_{nk}+B_{nk}}.$

By Theorem 4 of Miller, (2021), we have the pointwise Laplace approximation

$\displaystyle\log
Z_{n}^{\rho}\to\log\pi(\beta_{0})+\frac{P\log(\rho)}{2}+\frac{1}{2}\log|2\pi
H(\beta_{0})|\qquad\text{almost surely}.$

Hence, $A_{nk}\to\frac{P}{2}\log(\lambda_{k}/\lambda_{k+1})$ and
$B_{nk}\to\frac{-\delta\lambda_{k+1}}{\lambda_{k}}\left|\log\pi(\beta_{0})+\frac{P\log\lambda_{k}}{2}+\frac{1}{2}\log|2\pi
H(\beta_{0})|\right|$. Letting

$\displaystyle A=-\frac{P}{2}\frac{\delta}{a}\qquad\text{and}\qquad
B=-\delta\left\\{|\log\pi(\beta_{0})|+\frac{P(|\log a|+|\log
b|)}{2}+\frac{1}{2}|2\pi H(\beta_{0})|\right\\}$

we have $\lim A_{nk}\geq A$ and $\lim B_{nk}\geq B$ almost surely. Therefore

$\displaystyle\limsup\sup_{\phi}\int|q_{n}^{\phi}-q_{n}^{\lambda^{-1}}|\lesssim
1-e^{A+B}=1-e^{-\delta C}\leq\delta C$

where $C=-(A+B)/\delta>0$ does not depend on $\delta$.

Putting Everything Together: From the above arguments, we have

$\displaystyle\limsup\sup_{\phi}\int|q^{\phi}_{n}-N^{\phi}|\
dh\lesssim\delta+0+\delta.$

Because this holds for all $\delta>0$, we must have

$\displaystyle\lim_{n\to\infty}\sup_{\phi}\int|q^{\phi}_{n}-N^{\phi}|\ dh=0$

almost surely. ∎

We use this uniform Bernstein-von Mises theorem to establish the result for
the quasi-posterior defined in Section 4.

###### Proof.

We start by showing that the Bernstein-von Mises result applies to a slightly
different distribution, which is sampled from via the following algorithm:

  1. 1.

Initialize $\sigma$ with any value in $[a,b]$.

  2. 2.

Sample $\beta^{\prime}\sim\widetilde{\pi}(\beta^{\prime}\mid\sigma)$.

  3. 3.

Sample $\phi\sim g(\phi\mid\beta^{\prime})$.

  4. 4.

Sample $\beta\sim\widetilde{\pi}(\beta\mid\phi)$.

Let $m_{\sigma}$ denote the marginal distribution of $\phi$ under this
algorithm and $\widetilde{\pi}_{*}$ the marginal distribution of $\beta$ under
this algorithm. By the triangle inequality we have

$\displaystyle\int|\widetilde{\pi}_{*}-N(\cdot\mid\beta_{0},\phi_{0}N^{-1}\,H(\beta_{0})^{-1})|$
$\displaystyle\leq\int|\widetilde{\pi}_{*}-\widetilde{\pi}(\cdot\mid\phi_{0})|$
$\displaystyle\qquad+\int|\widetilde{\pi}(\cdot\mid\phi_{0})-N(\cdot\mid\beta_{0},\phi_{0}N^{-1}\,H(\beta_{0})^{-1})|.$

By Lemma 1, the second term converges to $0$ almost surely, so it suffices to
show that the first term
$A=\int|\widetilde{\pi}_{*}(\beta)-\widetilde{\pi}(\beta\mid\phi_{0})|\
d\beta$ also converges to $0$ almost surely. Next, we rewrite
$\widetilde{\pi}_{*}(\beta)$ as

$\displaystyle\int\widetilde{\pi}(\beta\mid\phi)\,m_{\sigma}(\phi)\
d\phi\qquad\text{where}\qquad m_{\sigma}(\phi)=\int
g(\phi\mid\beta^{\prime})\,\widetilde{\pi}(\beta^{\prime}\mid\sigma).$

By the triangle inequality we have

$\displaystyle A$
$\displaystyle\leq\int|\widetilde{\pi}(\beta\mid\phi)-\widetilde{\pi}(\beta\mid\phi_{0})|\,m_{\sigma}(\phi)\
d\phi\ d\beta$
$\displaystyle\leq\int|\pi(\beta\mid\phi)-N(\beta\mid\widehat{\beta},\phi\Sigma)|\,m_{\sigma}(\phi)+|N(\beta\mid\widehat{\beta},\phi\Sigma)-N(\beta\mid\widehat{\beta},\phi_{0}\Sigma)|\
m_{\sigma}(\phi)$
$\displaystyle\qquad+|\widetilde{\pi}(\beta\mid\phi_{0})-N(\beta\mid\widehat{\beta},\phi_{0}\Sigma)|m_{\sigma}(\phi)\
d\phi\ d\beta$

where $\Sigma=(NH(\beta_{0}))^{-1}$. By another application of Lemma 1, it
suffices to show that $\int|N^{\phi}-N^{\phi_{0}}|\,m_{\sigma}(\phi)\ d\phi\
d\beta\to 0$ almost surely, where $N^{\phi}$ is defined as in proof of Lemma
1. To show this, we show that $m_{\sigma}(\phi)\to\delta_{\phi_{0}}(\phi)$ in
distribution (almost surely). Letting $\psi(\phi)$ denote an arbitrary bounded
and continuous function, we have

$\displaystyle\int\psi(\phi)\,m_{\sigma}(\phi)\ d\phi$
$\displaystyle=\iint\psi(\phi)\,g(\phi\mid\widehat{\beta}+h/\sqrt{n})q_{n}^{\sigma}(h)\
d\phi\ dh$
$\displaystyle=\int\mathbb{E}_{g}\\{\psi(\phi)\mid\widehat{\beta}+h/\sqrt{n}\\}q_{n}^{\sigma}(h)\
dh$

where $q_{n}^{\sigma}$ is as defined in the proof of Lemma 1. Because
$\psi(\phi)$ is bounded, $q_{n}^{\sigma}\to N^{\sigma}$ pointwise by the
Laplace approximation, and $\int q_{n}^{\sigma}\ dh\to\int N^{\sigma}\ dh$ we
can apply the generalized dominated convergence theorem to conclude that

$\displaystyle\lim_{n\to\infty}\int\psi(\phi)\,m(\phi)\
d\phi=\int\lim_{n\to\infty}\mathbb{E}_{g}\\{\psi(\phi)\mid\widehat{\beta}+h/\sqrt{n}\\}q^{\sigma}(h)\
d\phi$

almost surely. By the truncated gamma specification of $g(\phi\mid\beta)$ we
have

$\displaystyle[\phi\mid\widehat{\beta}+h/\sqrt{n}]\sim\text{InvGam}\left\\{\frac{n}{2},\frac{n\widehat{\phi}(\widehat{\beta}+h/\sqrt{n})}{2}\right\\}\text{
truncated to }[a,b].$

By Condition G (i.e., $\widehat{\phi}(\widehat{\beta}+h/\sqrt{n})\to\phi_{0}$
for all $h$ almost surely) and the fact that $\phi_{0}\in[a,b]$ by assumption,
it is straight-forward to show that this distribution converges to a
$\delta_{\phi_{0}}$ distribution (almost surely), and hence
$\lim_{n\to\infty}\mathbb{E}_{g}\\{\psi(\phi)\mid\widehat{\beta}+h/\sqrt{n}\\}=\psi(\phi_{0})$
almost surely. Therefore

$\displaystyle\lim_{n\to\infty}\int\psi(\phi)\,m(\phi)\
d\phi=\int\psi(\phi_{0})N^{\sigma}(h)\ dh=\psi(\phi_{0}),$

so that $m_{\sigma}(\phi)\to\delta_{\phi_{0}}(\phi)$ in distribution (almost
surely). It follows that

$\displaystyle\int|N^{\phi}-N^{\phi_{0}}|m_{\sigma}(\phi)\ dh\
d\phi=\int\Xi(\phi)\,m_{\sigma}(\phi)\ d\phi\to\Xi(\phi_{0})=0$

where $\Xi(\phi)=\int|N^{\phi}-N^{\phi_{0}}|\ dh$.

Having established the result for $\widetilde{\pi}_{*}$, we now establish the
result for $\widetilde{\pi}$. By an identical argument to the one above, it
suffices to show that $m(\phi)\to\delta_{\phi_{0}}(\phi)$ in distribution
(almost surely), which is also part of the conclusion of the theorem. We use a
strategy similar to the one used to prove Lemma 1. Let
$b^{-1}=\lambda_{0}<\lambda_{1}<\cdots<\lambda_{K}=a^{-1}$ denote a
$\delta$-net of $[b^{-1},a^{-1}]$ and for each $\phi$ let
$\sigma_{\phi}=\lambda_{k}^{-1}$ such that
$\lambda_{k}\leq\phi^{-1}\leq\lambda_{k+1}$. Using the integral representation
of $m(\phi)$ we have (for an arbitrary bounded and continuous function
$\psi(\phi)$)

$\displaystyle\int\psi(\phi)m(\phi)\
d\phi=\int\psi(\phi)\,g(\phi\mid\beta)\,\widetilde{\pi}(\beta\mid\rho)\
m(\rho)\ d\rho\ d\beta\ d\phi$

Adding and subtracting $\widetilde{\pi}(\beta\mid\sigma_{\rho})$ we get

$\displaystyle\int\psi(\phi)g(\phi\mid\beta)\\{\widetilde{\pi}(\beta\mid\rho)-\widetilde{\pi}(\beta\mid\sigma_{\rho})\\}\,m(\rho)\
d\rho\ d\beta\ d\phi$
$\displaystyle\qquad+\int\psi(\phi)g(\phi\mid\beta)\,\widetilde{\pi}(\beta\mid\sigma_{\rho})\,m(\rho)\
d\rho\ d\beta\ d\phi$

We conclude by showing that the first term can be made arbitrarily small while
the second term above tends to $\psi(\phi_{0})$. For the first term we apply
the triangle inequality, boundedness of $\psi(\phi_{0})$, and change of
variables to get

$\displaystyle\left|\int\psi(\phi)g(\phi\mid\beta)\\{\widetilde{\pi}(\beta\mid\rho)-\widetilde{\pi}(\beta\mid\sigma_{\rho})\\}\,m(\rho)\
d\rho\ d\beta\ d\phi\right|$
$\displaystyle\lesssim\int|q_{n}^{\rho}-q_{n}^{\sigma_{\rho}}|\,m(\rho)\ dh\
d\rho.$

By the argument given under the heading First Term in the proof of Lemma 1, we
know that $\limsup\sup_{\rho}\int|q_{n}^{\rho}-q_{n}^{\sigma_{\rho}}|\
dh\lesssim\delta$ almost surely so that

$\displaystyle\limsup\left|\int\psi(\phi)g(\phi\mid\beta)\\{\widetilde{\pi}(\beta\mid\rho)-\widetilde{\pi}(\beta\mid\sigma_{\rho})\\}\,m(\rho)\
d\rho\ d\beta\ d\phi\right|\lesssim\delta.$

For the second term, let $m\\{(a,b)\\}$ denote the mass $m(\phi)$ assigned to
the interval $(a,b)$ and observe that

$\displaystyle\int\psi(\phi)g(\phi\mid\beta)\widetilde{\pi}(\beta\mid\sigma_{\rho})\,m(\rho)\
d\rho\ d\beta\ d\phi$
$\displaystyle=\sum_{k}m\\{(\lambda_{k},\lambda_{k+1})\\}\int\psi(\phi)\,g(\phi\mid\beta)\widetilde{\pi}(\beta\mid\sigma_{k})\
d\beta\ d\phi.$
$\displaystyle=\sum_{k}m\\{(\lambda_{k},\lambda_{k+1})\\}\int\psi(\phi)\,m_{\sigma_{k}}(\phi)\
d\phi$
$\displaystyle\to\sum_{k}m\\{(\lambda_{k},\lambda_{k+1})\\}\psi(\phi_{0})=\psi(\phi_{0}),$

with the last line following from the fact that
$m_{\sigma_{k}}(\phi)\to\delta_{\phi_{0}}(\phi)$ in distribution (almost
surely).

Summarizing, we have

$\displaystyle\limsup\left|\int\psi(\phi)\,m(\phi)-\psi(\phi_{0})\right|$
$\displaystyle\qquad\leq\limsup\left|\int\psi(\phi)g(\phi\mid\beta)\\{\widetilde{\pi}(\beta\mid\rho)-\widetilde{\pi}(\beta\mid\sigma_{\rho})\\}\,m(\rho)\
d\rho\ d\beta\ d\phi\right|$
$\displaystyle\qquad\qquad+\left|\int\psi(\phi)g(\phi\mid\beta)\,\widetilde{\pi}(\beta\mid\sigma_{\rho})\,m(\rho)\
d\rho\ d\beta\ d\phi-\psi(\phi_{0})\right|$
$\displaystyle\qquad\lesssim\delta+0,$

for all $\delta>0$, so that the limit exists and is $0$. ∎