# Identifiability and Falsifiability: Two Challenges for Bayesian Model
Expansion

Collin Cademartori

###### Abstract

We study the identifiability of model parameters and falsifiability of model
predictions under conditions of model expansion in a Bayesian setting. We
present results and examples suggesting a tendency for identifiability and
falsifiability to decrease in this context and for the severity of these
problems to trade-off against one another. Additionally, we present two
extended examples that demonstrate how these difficulties can be partially
overcome by inferential methods that leverage the joint structure of the
posterior distribution.

##  1 Introduction

In this work we connect the process of (Bayesian) model expansion to two
challenges for the interpretation and evaluation of statistical models,
namely:

  * •

identifiability \- the ability of the model to support sufficiently precise
inferences about parameters of interest, and

  * •

falsifiability \- the readiness of the model to reveal deficiencies in its fit
to the observed data.

These general concepts can be made precise in different ways. In practice,
poor identifiability can manifest as marginal posterior distributions that are
too wide to support substantively interesting conclusions about quantities of
interest. Likewise, poor falsifiability can result in reduced power for tests
of model fitness compared to nearby alternatives. In the Bayesian context, we
will argue that identifiability can be quantified using the mutual information
$\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}\right)$, and that
falsifiability can be quantified by the conditional mutual information
$\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)$
\- quantities from information theory which we will discuss in more detail in
the following sections.

Model criticism has long been recognized as an essential component of applied
statistical workflow, and this process commonly creates a need to expand our
models to capture a more diverse collection of data behaviors [4, 10, 26].
However, this process is not without challenges, as higher dimensional models
can exhibit more complex posterior distributions which frustrate simple
conclusions. Our main result quantifies two of these challenges by showing
that under appropriate classes of model expansion, there exist bounds on these
quantities which exhibit two important properties: (i) a bias towards reducing
both falsifiability and identifiability as the dimension of the expanded model
grows and (ii) a tradeoff whereby model expansions which avoid reducing one of
identifiability and falsifiability are in some sense more likely to reduce the
other.

While we do not expect the behavior of bounds (which may be quite loose) to
translate directly or universally to individual models and datasets, these two
properties of our bounds qualitatively match the patterns that we observe in
both simple cases where direct calculations are possible and in more complex
examples with simulated and real data. We thus view the main contribution of
our result as conceptually uniting and generalizing patterns observed in
particular cases. For example:

  1. 1.

The literature on Bayesian sparse regression has demonstrated that the
identification problems inherent to high-dimensional regression problems can
often be alleviated by imposing certain hierarchical priors on the
coefficients (e.g. horseshoe or normal scale-mixture priors) [20, 19, 21].
Critically, these priors work by encoding dependence between the coefficients,
and thus do not require the addition of prior information marginally to be
effective. The tradeoff between identifiability and falsifiability we observe
in our main result breaks down when the prior encodes enough dependence
between the parameters. This suggests that the kind of dependence encoding
that can resolve identification problems in regression models may be a good
strategy for addressing identifiability deficits more generally.

  2. 2.

The posterior predictive $p$-value has been criticized in the model checking
literature as being conservative or under-powered [2, 23, 31]. Because these
criticisms have hinged on frequency properties of the $p$-value, some
Bayesians have responded by pointing out that the posterior predictive
$p$-value is interpretable without reference to its distribution under
frequentist replications. We will argue that our notion of falsifiability is
directly linked to a general concept of power which does not require any
reference to frequenctist considerations. Our main result suggests that the
risk of such conservaitity problems is directly linked to the process of model
expansion, but it also motivates a generlization of the posterior predictive
$p$-value which we show is capable of resolving some of the practical problems
caused by these problems.

Our overall conclusion is thus both negative and positive. On on hand, we
believe that our main result suggests a real tension exists between some of
the basic goals of applied modeling in the context of iterative model
expansion. On the other hand, we do not believe our result militates against
successful model expansion in general. Rather, by quantifying some features of
this tradeoff, our result points towards possible tools which we believe can
form the basis of an expansion-ready statistical methodology.

We illustrate the basic shape of this tradeoff with an extremely simple
regression example to establish intuition. Suppose we have only two
observations $(y_{1},y_{2})$ and known measurement variance $\sigma^{2}=1$. In
our first model, we have one predictor $\mathbf{x}_{1}=(0,1)$ with coefficient
$\beta_{1}$. Assigning a normal prior, the resulting model is

$y_{j}\mid\beta_{1}\sim\mathrm{normal}\left(\beta_{1}x_{1j},1\right)\text{ for
}j=1,2,\quad\beta_{1}\sim\mathrm{normal}\left(0,\sigma_{b}\right),$ (1)

where the hyperparameter $\sigma_{b}$ is taken large so that the prior is
weakly informative. We then expand this model by adding a second predictor
$\mathbf{x}_{2}$ with $\|\mathbf{x}_{2}\|=1$ and with coefficient $\beta_{2}$.
Assuming $\beta_{2}$ is a priori independent of $\beta_{1}$ and assigning an
identical marginal prior, we get

$y_{j}\mid\beta_{1},\beta_{2}\sim\mathrm{normal}\left(\beta_{1}x_{1j}+\beta_{2}x_{2j},1\right)\text{
for }j=1,2,\quad\beta_{1},\beta_{2}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(0,\sigma_{b}\right).$ (2)

We consider five levels of nonnegative correlation between the predictors
$\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, which are plotted in the top row of
Figure 1. To assess the effect of adding a second predictor on
identifiability, we compare the marginal posterior of $\beta_{1}$ to its
prior, plotted in the second row of Figure 1.

As we would expect, as predictor correlation increases, the identification of
the coefficient $\beta_{1}$ decreases.

It is less obvious how we should assess the falsifiability of the model. We
argue in Section 3 that falsifiability is connected to a measure of posterior
confidence about the true data generating process, expressed as a distribution
over independent, replicated data $\mathbf{y}_{\mathrm{rep}}$.

In particular, we will argue that falsifiability tends to decrease as the
sampling distributions $p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\beta})$
and the posterior predictive distribution
$p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=\int
p\left(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\beta}\right)p\left(\boldsymbol{\beta}\mid\mathbf{y}\right)d\boldsymbol{\theta}$
become more dissimilar.

The third row of Figure 1 partially visualizes this by plotting the sampling
distributions for a replicated first observation $y^{\mathrm{rep}}_{1}$ at the
posterior means $\overline{\boldsymbol{\beta}}$ with the corresponding
posterior predictive distributions.

As the correlation between the predictors decreases, the distributions
$p(y_{1}^{\mathrm{rep}}\mid\overline{\boldsymbol{\beta}})$ and
$p(y_{1}^{\mathrm{rep}}\mid\mathbf{y})$ become less similar.

Figure 1: First row: $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, in order of
increasing correlation. Second row: the priors $p(\beta_{1})$ (blue) and the
posteriors $p(\beta_{1}\mid\mathbf{y})$ (red), both centered to allow for an
easier comparison of scales. Narrower posteriors relative to the prior
indicate better identification. Third row: the posterior predictive
$p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})$ (red) and the sampling
distributions $p(\mathbf{y}\mid\overline{\boldsymbol{\beta}})$ (blue). More
dissimilarity is connected to lower model check power on average.

The relationships displayed in the highlighted panels correspond to the
relationships between the corresponding distributions that occur in the
single-predictor model (1). These highlighted panels also correspond to the
best-case behavior, showing that the expanded model can only perform worse
than the base model on either metric. In fact, we further see that these
behaviors are inversely correlated among the expanded models, i.e. the most
precise marginal inference occurs when the sampling and posterior predictive
distributions are most dissimilar and vice versa. We will find evidence of
similar phenomena more generally in the next sections.

###  1.1 Outline

In Sections 2 and 3, we separately connect weakening of identifiability and
falsifiability to some general conditions of model expansion . Section 4
presents our main result connecting identifiability and falsifiability,
showing that expansions which decrease the severity of one challenge have an
increased risk of worsening the other (as seen in the above regression
example). In Sections 5 and 6, we demonstrate through two examples how richer
inferences which leverage the full joint structure of the posterior
distribution can alleviate some of the difficulties imposed by poor
identifiability and falsifiability.

###  1.2 Contributions

Our main result, Theorem 5 uses a generalizaton of the Bayesian Cramer-Rao
bound [1] and matrix concentration inequalities [25] to establish bounds on
the mutual and conditional mutual information which exhibit both a tradeoff
and dimension dependence. We sketch the result here. Suppose that we have a
base model and expasion thereof denoted $p_{\mathrm{base}}$ and $p$
respectively, both defined over common data $\mathbf{y}$ and with some shared
parameters $\boldsymbol{\theta}$ (a notion which we will define unambiguously
in the next section).

######  Theorem 1.

Let $\boldsymbol{\iota}_{\mathrm{base}}$ and $\boldsymbol{\iota}$ be the
eigivenvalues of the expected Fisher information matrices for
$\boldsymbol{\theta}$ in the base and expanded models $p_{\mathrm{base}}$ and
$p$, respectively. Furthermore let $\mathbf{I}_{\mathrm{base}}$ and
$\mathbf{I}$ denote the (conditional) mutual information in the base and
expanded models. Finally, let $d$ be the dimension of $\boldsymbol{\theta}$
and $d^{\mathrm{exp}}\geq d$ the dimension of the parameter space for the
expanded model. Then, under technical conditions given in the statement of
Theorem 5, we have

$\underbrace{\mathbf{I}_{\mathrm{base}}\left(\mathbf{y},\boldsymbol{\theta}\right)}_{\begin{subarray}{c}\text{base
model}\\\ \text{identifiability of
$\boldsymbol{\theta}$}\end{subarray}}\leq\Psi_{i}\left(\boldsymbol{\iota}_{\mathrm{base}}\right),\qquad\underbrace{\mathbf{I}\left(\mathbf{y},\boldsymbol{\theta}\right)}_{\begin{subarray}{c}\text{expanded
model}\\\ \text{identifiability of
$\boldsymbol{\theta}$}\end{subarray}}\leq\Psi_{i}\left(\boldsymbol{\iota}\right)-\Delta_{i}$
(3)

$\underbrace{\mathbf{I}_{\mathrm{base}}\left(\mathbf{y}_{\mathrm{rep}},\boldsymbol{\theta}\mid\mathbf{y}\right)}_{\begin{subarray}{c}\text{base
model}\\\ \text{falsifiability}\end{subarray}}\geq
c_{d}\Psi_{f}\left(\boldsymbol{\iota}_{\mathrm{base}}\right)\qquad\underbrace{\mathbf{I}\left(\mathbf{y}_{\mathrm{rep}},\boldsymbol{\theta}\mid\mathbf{y}\right)}_{\begin{subarray}{c}\text{expanded
model}\\\ \text{falsifiability}\end{subarray}}\geq
c_{d^{\mathrm{exp}}}\left[\Psi_{f}\left(\boldsymbol{\iota}\right)+\Delta_{f}\right],$
(4)

where $\Psi_{i},\Psi_{f}$ are increasing in each of the components of the
vector argument, $\Delta_{m},\Delta_{c}\geq 0$ are terms which tend to
increase in magnitude with $d^{\mathrm{exp}}$, and
$c_{d},c_{d^{\mathrm{exp}}}$ are constants depending only on $d$ and
$d^{\mathrm{exp}}$ respectively.

These inequalities are given from the adverse directions, in the sense that
smaller mutual information and larger conditional mutual information are
associated with reduced identifiability and falsifiability respectively. The
dimensional dependence enters through the $\Delta$ terms, which push the
bounds in the corresponding adverse directions. The tradeoff between these
bounds occurs through the $\Psi$ terms. For instance, if the components of
$\boldsymbol{\iota}$ are all smaller than
$\boldsymbol{\iota}_{\mathrm{base}}$, then our mutual information bound will
decrease in passing to the expanded model. In the reverse case, the
conditional mutual information bound will increase.

In conjunction with numerous examples, this result suggests that at least one
of reduced identifiability and reduced falsifiability should be expected in
the process of iterative model expansion, which is the first major
contribution of this work. The generality of this phenomenon motivates
considering methods of inference and model checking which can cope with these
conditions by extracting as much useful information from our models as
possible.

Our second contribution is to demonstrate methods of utilizing the joint
structure of the posterior in practice which may allow these challenges to be
partially overcome in many cases. We show in an extended example how the
dependence structure of the posterior can contain significant, practically
useful information even when the marginal inferences are too weak to support
strong conclusions about individual parameters. And in the context of model
checking, we provide an extension of the traditional posterior predictive
$p$-value, which we validate in a real data example, and which we argue is
often more useful and easily applied than previous solutions designed to
resolve the posterior predictive $p$-value’s claimed power deficiencies.

###  1.3 Related Work

Recently, statistical workflow has enjoyed increased attention as a discrete
topic in statistics. This literature has sought to provide a consistent
framework and practical advice for each step of a statistical analysis,
including the process of model expansion (see, e.g. [10, 26, 7]). Here we seek
to complement this perspective by studying model expansion as a distinct
regime. To this end, our main result provides interpretable bounds on the
mutual and conditional mutual information, the former of which depends
critically on Theorem 2 in [1].

Outside of the context of model expansion, the problem of weak/non-
identification has been extensively studied in the classical and Bayesian
contexts. In the Bayesian setting, methods of detecting and dealing with
identification problems have been studied in, e.g. [29, 16]. Whereas these
methods have usually been tied to particular (classes of) models, we study
this problem in a general setting of model expansion.

As we will argue in Section 3, problems of falsifiability are directly
connected to debates over the power and conservativity of the posterior
predictive $p$-value. Various forms of this problem have been described, and
possible solutions have been proposed in [3, 2, 23, 31]. We propose another
possible solution - conditional $p$-values - which differ from these previous
proposals both in their goal and method of use, and we will argue that our
approach is more practically applicable in many cases.

Our approach to studying the problems of identifiability and falsifiability
follows many previous successes in using information-theoretic tools to
understand and quantify model behaviors in great generality. We enumerate a
few connections of particular note:

  1. 1.

We quantify identifiability by thinking of the information entropy of a
posterior as representing our uncertainty about parameters of interest. This
representation of uncertainty as entropy can be traced back to Jaynes, who
used it to justify the use of maximum entropy posterior distributions [15].

  2. 2.

Information-theoretic criteria have long been used to evaluate the predictive
performance of models [28]. In the Bayesian context, the expected log
predictive density (ELPD) has been used as a flexible and model-specific
objective for model evaluation and comparison [27]. When our data consists of
a scalar quantity $y$, and the model is correctly specified, the ELPD can be
given as

$\mathrm{D}\left(p(y\mid\boldsymbol{\theta}^{*})\;||\;p(y\mid
y^{\mathrm{rep}})\right)+C,$

where $p(\mathbf{y}\mid\boldsymbol{\theta}^{*})$ is the true data generating
process, $C$ is a constant depending only on this true distribution, and
$y^{\mathrm{rep}}\sim p(y\mid\boldsymbol{\theta})$ is an independent
replication of the data. If we substitute the true value
$\boldsymbol{\theta}^{*}$ with an average over the posterior
$p(\boldsymbol{\theta}\mid\mathbf{y})$, then the first term recovers the
conditional mutual information
$\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)$
which we relate to the falsifiability of the model.

  3. 3.

The Rashomon effect, first described by Breiman in [5], is a phenomenon
whereby many models can achieve similar overall loss yet provide very
different point predictions. In our work, we find that our concept of
falsifiability is also threatened by the multiplicity of plausible sampling
distributions in a model. And indeed, our conditional mutual information rests
on a conceptually simiar KL divergence as a recently proposed metric for
quantifying the Rashomon effect, the Rashomon capacity [12].

  4. 4.

Mutual information-based quantities have also been deployed to bound measures
of other adverse model behaviors, particular bias and generalization error
[24, 30].

##  2 Weak Identifiability and Model Expansion

We start by defining the types of model expansions to which our results will
apply. We will write $p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})$ for
some base model defined over data $\mathbf{y}\in\mathbb{R}^{n}$ and parameters
$\boldsymbol{\theta}\in\mathbb{R}^{d}$. We then consider certain expansions of
this base model defined as follows.

######  Definition 1 (Model Expansion).

A model $p(\mathbf{y},\boldsymbol{\theta},\boldsymbol{\lambda})$ defined with
additional parameter $\boldsymbol{\lambda}\in\overline{\mathbb{R}}^{k}$ is an
expansion of $p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})$ if

$p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})=p(\mathbf{y},\boldsymbol{\theta}\mid\boldsymbol{\lambda}_{0})\text{
for }\boldsymbol{\lambda}_{0}\in\overline{\mathbb{R}}^{k},$ (5)

where $\overline{\mathbb{R}}=[-\infty,\infty]$.

This framework includes many common examples of model expansion:

  * •

Let $p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})$ be a generalized
linear model with response vector $\mathbf{y}$ and parameters
$\boldsymbol{\theta}$ including the coefficients and any additional
parameters. Adding a new predictor and coefficient $\lambda$ with independent
prior is then an expansion since
$p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})=p(\mathbf{y},\boldsymbol{\theta}\mid\lambda=0)$.

  * •

Let $p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})$ be an exchangeable
Poisson model over the data $y_{i}$ with $\theta$ the Poisson rate. We can
extend this with an overdispersion parameter $\lambda$ (with independent
prior). This is commonly modeled with a negative binomial distribution

$p(y\mid\theta,\lambda)=\binom{y+\lambda-1}{y}\left(\frac{\theta}{\theta+\lambda}\right)^{y}\left(\frac{\lambda}{\theta+\lambda}\right)^{\lambda}.$

Since this reduces to the Poisson as $\lambda\to\infty$, we have that
$p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})=p(\mathbf{y},\boldsymbol{\theta}\mid\lambda=\infty)$,
so this is again an expansion.

###  2.1 Weak Identification and Marginal Entropy

We now formalize our notion of identification using the information entropy.
First we establish some notation. For a joint model
$q(\boldsymbol{\theta},\mathbf{y})$, the (differential) entropy of
$q(\boldsymbol{\theta})$ is denoted
$h_{q(\boldsymbol{\theta})}(\boldsymbol{\theta})$, and the conditional entropy
of $\boldsymbol{\theta}$ given $\mathbf{y}$ is
$h_{q(\boldsymbol{\theta},\mathbf{y})}\left(\boldsymbol{\theta}\mid\mathbf{y}\right)$.
The mutual information (mi) is denoted
$\mathbf{I}_{q}\left(\boldsymbol{\theta},\mathbf{y}\right)$, which will at
times be extended to a conditional mutual information (cmi), denoted by
$\mathbf{I}_{q}\left(\boldsymbol{\theta},\mathbf{y}\mid\mathbf{x}\right)$,
when the joint model extends over an additional quantity $\mathbf{x}$. When
distributions are clear from context, we may drop subscripts from entropies
and mutual informations, writing e.g. $h(\boldsymbol{\theta})$ and
$\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}\right)$. The reader who is
unfamiliar with information theory may consult Appendix A for definitions of
these quantities and statements of the basic results that we use. With these
definitions, we can now give quantitative operational definitions of our
notions of weak marginal identification for arbitrary subsets of
$\boldsymbol{\theta}$.

######  Definition 2 ($\epsilon$-Weak Identification).

Let $I\subset[d]$. We say for any $\epsilon>0$ that
$\boldsymbol{\theta}_{I}=\left(\theta_{i}\right)_{i\in I}$ is
$\epsilon$-weakly identified for data $\mathbf{y}$ if

$h_{p(\boldsymbol{\theta}_{I}\mid\mathbf{y})}\left(\boldsymbol{\theta}_{I}\right)>h_{p(\boldsymbol{\theta}_{I})}\left(\boldsymbol{\theta}_{I}\right)-\epsilon.$
(6)

For $p_{1}\left(\boldsymbol{\theta},\mathbf{y}\right)$ and
$p_{2}(\boldsymbol{\theta},\mathbf{y})$, $\boldsymbol{\theta}_{I}$ is more
weakly identified in $p_{2}$ than $p_{1}$ if
$h_{p_{2}(\boldsymbol{\theta}_{I}\mid\mathbf{y})}\left(\boldsymbol{\theta}_{I}\right)>h_{p_{1}(\boldsymbol{\theta}_{I}\mid\mathbf{y})}\left(\boldsymbol{\theta}_{I}\right)$.

We also define weak identification for entire models (regardless of data
$\mathbf{y}$) by averaging over the prior predictive distribution.

######  Definition 3 ($\epsilon$-Weakly Identifiable Model).

Let $I\subset[d]$. We say $\boldsymbol{\theta}_{I}$ is $\epsilon$-weakly
identifiable in $p(\boldsymbol{\theta},\mathbf{y})$ if

$h\left(\boldsymbol{\theta}\mid\mathbf{y}\right)>h\left(\boldsymbol{\theta}\right)-\epsilon,$
(7)

or, equivalently, if $\mathbf{I}(\boldsymbol{\theta},\mathbf{y})<\epsilon$.

Henceforth, we will leave the $\epsilon$-dependence of this definition
implicit and simply say that a parameter is weakly identified if it is
$\epsilon$-weakly identified for an appropriate value of $\epsilon$ (which
will usually be given by domain understanding).

This operational definition of weak identification can only be interpreted
relative to the prior. In many cases, this is a natural quantity to focus on
(e.g. when we are concerned with the cost-benefit tradeoffs of data collection
or the contribution of a research finding to existing knowledge). However, if
we expand a model by adding prior information about $\boldsymbol{\theta}$,
then it possible for both the posterior entropy of $\boldsymbol{\theta}$ and
the mutual information to decrease. In other words, the identification
relative to the prior may decrease while the posterior becomes more
concentrated. This divergence between absolute and relative notions of
identification can be avoided if we exclude from consideration expansions
which decrease $h(\boldsymbol{\theta})$.

We can now show that certain model expansions tend to weaken identification in
the above sense. If $p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})$
is an expansion of $p_{\mathrm{base}}(\boldsymbol{\theta},\mathbf{y})$, then
we have the following decomposition of the mutual information:

$\mathbf{I}_{p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})}(\boldsymbol{\theta},\mathbf{y})=\mathbf{I}_{p_{\mathrm{base}}(\boldsymbol{\theta},\mathbf{y})}(\boldsymbol{\theta},\mathbf{y})+\Delta_{I}^{\mathrm{exp}}+\Delta^{\mathrm{post}}_{I},$
(8)

where we define

$\displaystyle\Delta_{I}^{\mathrm{exp}}$
$\displaystyle=\mathbf{I}_{p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})}\left(\boldsymbol{\theta},\mathbf{y}\mid\boldsymbol{\lambda}\right)-\mathbf{I}_{p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})}\left(\boldsymbol{\theta},\mathbf{y}\mid\boldsymbol{\lambda}_{0}\right)$
$\displaystyle\Delta^{\mathrm{post}}_{I}$
$\displaystyle=\mathbf{I}_{p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})}\left(\boldsymbol{\lambda},\boldsymbol{\theta}\right)-\mathbf{I}_{p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})}(\boldsymbol{\lambda},\boldsymbol{\theta}\mid\mathbf{y}).$

The $\Delta_{I}^{\mathrm{exp}}$ term is the difference in amount of
information $\mathbf{y}$ provides about $\boldsymbol{\theta}$ given
$\boldsymbol{\lambda}$ and given $\boldsymbol{\lambda}_{0}$, averaging
$\boldsymbol{\lambda}$ over the expanded model. The
$\Delta^{\mathrm{post}}_{I}$ term is the difference in the amount of
information $\boldsymbol{\lambda}$ provides about $\boldsymbol{\theta}$ before
and after observing the data $\mathbf{y}$. We regard this as a measure of the
a priori informativeness of $\boldsymbol{\lambda}$ about
$\boldsymbol{\theta}$, which is justified by the fact that we have,
$\Delta^{\mathrm{post}}_{I}\geq-\mathbf{I}\left(\boldsymbol{\lambda},\boldsymbol{\theta}\mid\mathbf{y}\right)$
with equality if and only if $\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$
are independent in the expanded model (i.e. if
$p(\boldsymbol{\theta},\boldsymbol{\lambda})=p(\boldsymbol{\theta})p(\boldsymbol{\lambda})$).
When $\Delta^{\mathrm{post}}_{I}<0$, (8) shows that the expanded model is
biased towards weaker identification of $\boldsymbol{\theta}$ compared to the
base model.

We can also use the decomposition (8) to define a concept which will be useful
in the next sections. We say the parameter $\boldsymbol{\lambda}$ dilutes the
effect of shared parameter $\boldsymbol{\theta}$ if
$h_{p(\boldsymbol{\theta},\mathbf{y}\mid\boldsymbol{\lambda})}\left(\boldsymbol{\theta}\mid\mathbf{y}\right)$
is larger than
$h_{p_{\mathrm{base}}(\boldsymbol{\theta},\mathbf{y})}(\boldsymbol{\theta}\mid\mathbf{y})$
on average over $p(\boldsymbol{\lambda})$. If this relationship is reversed,
we say that $\boldsymbol{\lambda}$ concentrates the effect of
$\boldsymbol{\theta}$. In the case that
$p(\boldsymbol{\theta}\mid\boldsymbol{\lambda})=p_{\mathrm{base}}(\boldsymbol{\theta})$
for all $\boldsymbol{\lambda}$, dilution and concentration are equivalent to
$\Delta_{I}^{\mathrm{exp}}<0$ and $\Delta_{I}^{\mathrm{exp}}>0$ respectively.

###  2.2 The Relation to Marginal Fisher Information

For model $q(\mathbf{y},\boldsymbol{\theta})$, the observed and Fisher
information matrices are defined as

$\displaystyle[\boldsymbol{\mathcal{J}}_{q}(\mathbf{y},\boldsymbol{\theta})]_{ij}=-\frac{\partial^{2}}{\partial\theta_{i}\partial\theta_{j}}\log
q(\mathbf{y}\mid\boldsymbol{\theta})\text{ for all }1\leq i,j\leq d,$
$\displaystyle\boldsymbol{\mathcal{I}}_{q}(\boldsymbol{\theta})=\mathbb{E}_{q(\mathbf{y}\mid\boldsymbol{\theta})}\boldsymbol{\mathcal{J}}_{q}(\mathbf{y},\boldsymbol{\theta}).$
(9)

We drop the subscript when the model is clear from context. We now state a
bound on the mutual information in terms of the Fisher information, which
follows directly from Theorem 2 of Aras et al. [1].

######  Theorem 2 (Mutual Information Upper Bound).

Let $q(\mathbf{y},\boldsymbol{\theta})$ be a model such that the prior
$p(\boldsymbol{\theta})$ is log-concave with covariance matrix
$\boldsymbol{\Sigma}$, then

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y})\leq
d\psi\left(\frac{1}{d}\mathrm{tr}\left(\mathbb{E}_{q(\boldsymbol{\theta})}\boldsymbol{\Sigma}^{1/2}\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\boldsymbol{\Sigma}^{1/2}\right)\right),$
(10)

where $\psi(x)$ is the concave increasing function given by

$\psi(x)=\begin{cases}\sqrt{x},&0\leq x\leq 1\\\
1+\frac{1}{2}\log(x),&x>1\end{cases}.$

Now let $v_{\mathrm{pr}}$ be the maximum eigenvalue of $\boldsymbol{\Sigma}$
(the covariance matrix over just the $\boldsymbol{\theta}$ parameters). Then
we also clearly have

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y})\leq
d\psi\left(\frac{v_{\mathrm{pr}}}{d}\mathbb{E}_{q(\boldsymbol{\theta})}\mathrm{tr}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)\right).$
(11)

If $v_{\mathrm{pr}}$ differs between a base model and expanded model, then we
can rescale the prior over $\boldsymbol{\theta}$ in the expanded model so that
they are equal. The only possible difficulty is that we may no longer have a
$\boldsymbol{\lambda}_{0}$ for which
$p_{\mathrm{base}}\left(\boldsymbol{\theta}\right)=p\left(\boldsymbol{\theta}\mid\boldsymbol{\lambda}_{0}\right)$.
Such situations can always be resolved however by passing to a larger model
which includes a prior scale hyperparameter for $\boldsymbol{\theta}$ in
$\boldsymbol{\lambda}$. With such a hyperparameter, we always have the ability
to set both the marginal prior scale and
$\boldsymbol{\lambda}_{0}$-conditional prior scale for $\boldsymbol{\theta}$
independently, allowing equality of $v_{\mathrm{pr}}$ and preservation of the
model expansion property.

We also note that rescaling $\boldsymbol{\theta}$ leaves
$\mathbf{I}(\boldsymbol{\theta},\mathbf{y})$ unchanged since the mutual
information is invariant to all inveritble transformations of
$\boldsymbol{\theta}$ and $\mathbf{y}$ separately. Thus, with loss of little
generality, we henceforth assume that $v_{\mathrm{pr}}=1$ for all models. The
weaker bound (11) will be useful for comparisons to other quantities in the
next sections, and for deriving the following further upper bound, which
applies more directly to model expansions, is easier to compute, and mirrors
the relation (8).

######  Theorem 3.

Define the partial Hessian with respect to $\boldsymbol{\lambda}$ as
$\left[\mathbf{H}\left(\boldsymbol{\lambda};\boldsymbol{\theta},\mathbf{y}\right)\right]_{jk}=-\frac{\partial^{2}}{\partial\lambda_{j}\partial\lambda_{k}}\log
p(\mathbf{y},\boldsymbol{\theta},\boldsymbol{\lambda})$ for $1\leq j,k\leq m$.
Then under the regularity conditions in Appendix E,

$\mathbb{E}\mathrm{tr}\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\leq\sum_{j=1}^{d}\left[\mathbb{E}\left\\{-\frac{\partial^{2}}{\partial\theta_{j}^{2}}\log
p(\mathbf{y}\mid\boldsymbol{\theta},\boldsymbol{\lambda})\right\\}+\Delta_{j}\right],$
(12)

where we define

$\Delta_{j}=\mathbb{E}\left\\{-\frac{\partial^{2}}{\partial\theta_{i}^{2}}\log
p\left(\boldsymbol{\lambda}\mid\boldsymbol{\theta}\right)\right\\}-\frac{\left[\sum_{j=1}^{m}\mathbb{E}\frac{\partial}{\partial\lambda_{j}}\frac{\partial}{\partial\theta_{i}}\log
p(\mathbf{y}\mid\boldsymbol{\theta},\boldsymbol{\lambda})\right]^{2}}{\mathbb{E}\|\mathbf{H}\left(\boldsymbol{\lambda};\boldsymbol{\theta},\mathbf{y}\right)\|_{\mathrm{op}}}.$

The $\Delta_{j}$ terms compare prior- and likelihood-based measures of
dependence between $\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$ and are
thus analogous to the $\Delta^{\mathrm{post}}$ term in (8). As with the
$\Delta^{\mathrm{post}}$ term, we have $\Delta_{j}\leq 0$ when
$\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$ are independent under the
prior. In particular, when
$p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})$ is an expansion of
some base model and when $\sum_{i=1}^{d}\Delta_{j}<0$, (12) again exhibits a
downward bias on the Fisher information of the expanded model compared to the
base model. We now work through two simple examples to illustrate.

  1. 1.

Take a linear regression model $p_{\mathrm{base}}$ with response
$\mathbf{y}\in\mathbb{R}^{n}$, predictors $\mathbf{X}\in\mathbb{R}^{n\times
m}$, coefficients $\boldsymbol{\beta}\in\mathbb{R}^{m}$, intercept $\alpha$,
and log noise variance $\tau$:

$\left(2\pi\exp\left(\tau\right)\right)^{-n/2}\exp\left[-\left(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}-\alpha\mathbf{1}\right)^{T}\left(\mathbf{y}-\mathbf{X}\boldsymbol{\beta}-\alpha\mathbf{1}\right)/2\exp(\tau)\right].$

We consider an expansion $p$ with additional predictor
$\mathbf{z}\in\mathbb{R}^{n}$ and coefficient $\lambda$. Suppose the
coefficients are assigned independent priors, and let
$\boldsymbol{\theta}=(\sigma,\beta_{1},\ldots,\beta_{m},\alpha)$.

We assume without loss of generality that all predictors are centered as this
does not affect the posterior entropy. We then find for all $1\leq j\leq m$,

$\displaystyle-\mathbb{E}_{p}\frac{\partial^{2}}{\partial\beta_{j}^{2}}\log
p(\mathbf{y}\mid\lambda,\boldsymbol{\theta})=\mathbb{E}_{p(\tau)}\left\\{\frac{n\mathrm{var}\left(\mathbf{x}_{j}\right)}{\exp(\tau)}\right\\}=-\mathbb{E}_{p_{\mathrm{base}}}\frac{\partial^{2}}{\partial\beta_{j}^{2}}\log
p_{\mathrm{base}}(\mathbf{y}\mid\boldsymbol{\theta}),$
$\displaystyle-\mathbb{E}_{p}\frac{\partial^{2}}{\partial\tau}\log
p(\mathbf{y}\mid\lambda,\boldsymbol{\theta})=\frac{n}{2}=-\mathbb{E}_{p_{\mathrm{base}}}\frac{\partial^{2}}{\partial\tau^{2}}\log
p_{\mathrm{base}}(\mathbf{y}\mid\boldsymbol{\theta}),\text{ and }$
$\displaystyle-\mathbb{E}_{p}\frac{\partial^{2}}{\partial\alpha^{2}}\log
p(\mathbf{y}\mid\lambda,\boldsymbol{\theta})=\mathbb{E}_{p(\tau)}\left\\{\frac{n}{\exp(\tau)}\right\\}=-\mathbb{E}_{p_{\mathrm{base}}}\frac{\partial^{2}}{\partial\alpha^{2}}\log
p_{\mathrm{base}}(\mathbf{y}\mid\lambda,\boldsymbol{\theta})$

These computations show that the first term in (12) is just
$\mathrm{tr}\left(\mathbb{E}\boldsymbol{\mathcal{I}}_{p_{\mathrm{base}}}(\boldsymbol{\theta})\right)$.
Assuming $\mathbf{z}$ is also centered, computing the second term in (12)
gives that

$\mathbb{E}_{p_{\mathrm{base}}}\mathrm{Tr}\left(\boldsymbol{\mathcal{I}}_{p_{\mathrm{base}}}\right)-\mathbb{E}_{p}\mathrm{Tr}\left(\boldsymbol{\mathcal{I}}_{p}\right)\geq
n^{2}\mathbb{E}_{p(\tau)}\left\\{\frac{1}{\exp(\tau)}\right\\}\sum_{j=1}^{m}\left[\frac{\mathrm{cov}\left(\mathbf{x}_{j},\mathbf{z}\right)}{\mathrm{var}\left(\mathbf{z}\right)}\right]^{2},$

which reflects the familiar fact that the identifiability of regression models
is reduced by significant correlation between predictors.

  2. 2.

Next consider an exchangeable Poisson base model with likelihood

$\exp\left(\mu
n\overline{\mathbf{y}}-n\exp(\mu)\right)\Big{/}\left(y_{1}!\times
y_{2}!\times\cdots\times y_{n}!\right).$

This can be expanded to a negative binomial model with likelihood

$\left[\prod_{i=1}^{n}\frac{\Gamma\left(y_{i}+\exp(\lambda)\right)}{\Gamma\left(y_{i}+1\right)\Gamma\left(\exp(\lambda)\right)}\right]\left(\frac{\exp(\mu)}{\exp(\mu)+\exp(\lambda)}\right)^{\sum_{i=1}^{n}y_{i}}\left(\frac{\exp(\lambda)}{\exp(\mu)+\exp(\lambda)}\right)^{n\exp(\lambda)}.$

This converges to the Poisson density as $\lambda\to\infty$, so the expanded
model is an expansion of the Poisson model. Next observe that the second
derivatives with respect to $\mu$ are given by

$-\frac{\partial^{2}}{\partial\mu^{2}}\log
p\left(\mathbf{y}\mid\mu,\lambda\right)=n\exp\left(\mu\right)\left[1-\frac{\exp(\mu)}{\exp(\mu)+\exp(\lambda)}\right]\left[\frac{\overline{y}+\exp(\lambda)}{\exp(\mu)+\exp(\lambda)}\right],$

which has expected value
$n\exp\left(\mu\right)\left[1-\frac{\exp(\mu)}{\exp(\mu)+\exp(\lambda)}\right]$
under $p(\mathbf{y}\mid\mu,\lambda)$. With this we can show using (12) that
the Fisher information trace must fall in passing from the base to the
expanded model:

$\mathbb{E}\mathrm{Tr}\left(\boldsymbol{\mathcal{I}}_{p}\right)\leq\mathbb{E}\left\\{n\exp\left(\mu\right)\left[\frac{\exp(\lambda)}{\exp(\mu)+\exp(\lambda)}\right]\right\\}<\mathbb{E}\left\\{n\exp(\mu)\right\\}=\mathbb{E}\mathrm{Tr}\left(\boldsymbol{\mathcal{I}}_{p_{\mathrm{base}}}\right),$

The results of Theorems 2 and 3 both connect Bayesian and classical notions of
identification and show that the marginal mi is controlled by a quantity that
is often easily approximated before fitting the expanded model. This latter
property may be useful when posterior sampling is slow, providing an
indication of difficult posterior geometry before it frustrates the sampling
algorithm.

The results of this section, particularly (12), suggest that dilution of
$\boldsymbol{\theta}$ by $\boldsymbol{\lambda}$ may be heuristically indicated
by a positive difference:

$\boldsymbol{\Delta}_{\text{dilute}}=\mathbb{E}_{p(\boldsymbol{\theta})p(\boldsymbol{\lambda})}\left\\{\boldsymbol{\mathcal{I}}_{\mathrm{base}}\left(\boldsymbol{\theta}\right)-\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}\mid\boldsymbol{\lambda})\right\\},$

where $\boldsymbol{\mathcal{I}}_{\mathrm{base}}$ is the Fisher information of
the base model, and
$\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}\mid\boldsymbol{\lambda})$ is the
Fisher information of the expanded model conditional on
$\boldsymbol{\lambda}$, i.e. the principal submatrix of the full Fisher
information matrix
$\boldsymbol{\mathcal{I}}(\boldsymbol{\theta},\boldsymbol{\lambda})$ obtained
by deleting those twos and columns involving derivatives with respect to the
components of $\boldsymbol{\lambda}$. We will also say that the effect of
$\boldsymbol{\lambda}$ is totally diluting/concentrating of
$\boldsymbol{\theta}$ if $\boldsymbol{\Delta}_{\text{dilute}}$ is
positive/negative semidefinite, respectively.

##  3 Weak Falsifiability and Model Expansion

We now turn to the behavior of the posterior predictive distribution (ppd)
$p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})$ under model expansion. For joint
model $p(\mathbf{y},\boldsymbol{\theta})$, the ppd is

$p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=\int
p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})p(\boldsymbol{\theta}\mid\mathbf{y})d\boldsymbol{\theta}.$
(13)

Comparisons between posterior predictive samples and observed data are
commonly used to check Bayesian models. It is often convenient to formalize
these checks as posterior predictive $p$-values (ppp-vs).

######  Definition 4 (Posterior Predictive $p$-Value).

For observed data $\mathbf{y}$, joint model
$p(\mathbf{y},\boldsymbol{\theta})$, and real-valued test statistic $T$, the
right-tailed ppp-v for $T$ is

$p_{T}=\int_{\\{T(\mathbf{y}_{\mathrm{rep}})\geq
T(\mathbf{y})\\}}p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})d\mathbf{y}_{\mathrm{rep}}.$
(14)

The left-tailed and two-tailed $p$-values are defined analogously.

It will be useful to set ppp-vs in a general framework of model evaluations.

######  Definition 5 (Data Distribution Evaluation).

Given $p(\mathbf{y},\boldsymbol{\theta})$, let $\mathcal{Y}$ be the (common)
support of the densities $p(\mathbf{y}\mid\boldsymbol{\theta})$,
$P(\mathcal{Y})$ be the space of all densities on $\mathcal{Y}$, and
$\mathcal{E}\subset\mathbb{R}$. Then for any data $\mathbf{y}$, an evaluation
is a (measurable) map

$e_{\mathbf{y}}:P(\mathcal{Y})\to\mathcal{E}$

For posterior sample $\\{\boldsymbol{\theta}_{(s)}\\}_{s=1}^{S}$ and
evaluation $e_{\mathbf{y}}$, the modeler has data

$\left\\{\boldsymbol{\theta}_{(s)},e_{\mathbf{y}}\left(p(\mathbf{y}\mid\boldsymbol{\theta}_{(s)})\right)\right\\}_{s=1}^{S}.$
(15)

with which to evaluate the model. Since (15) can be complex and high-
dimensional, it may not proffer easy conclusions about overall model fitness.
Posterior predictive checks solve this by providing simple summaries of (15).
For statistic $T$, we define conditional ppp-vs $p_{T}(\boldsymbol{\theta})$
as the evaluations
$e_{\mathbf{y}}\left(p(\cdot\mid\boldsymbol{\theta})\right)$ for the map
$q(\cdot)\to\int_{\\{T(\mathbf{y}_{\mathrm{rep}})\geq
T(\mathbf{y})\\}}q(\mathbf{y}_{\mathrm{rep}})d\mathbf{y}_{\mathrm{rep}}$. The
usual ppp-v is then just the average:

$p_{T}=\int
e_{\mathbf{y}}\left(p(\cdot\mid\boldsymbol{\theta})\right)p(\boldsymbol{\theta}\mid\mathbf{y})d\boldsymbol{\theta},$
(16)

which is naturally estimated by
$\frac{1}{S}\sum_{s=1}^{S}e_{\mathbf{y}}\left(p(\cdot\mid\boldsymbol{\theta}_{(s)})\right)$.
Thus, the ppp-v may be limited if relevant information in (15) is lost in
(16). If
$p(\mathbf{y}\mid\boldsymbol{\theta}_{(s)})=p(\mathbf{y}\mid\boldsymbol{\theta}_{(t)})$
for all $\mathbf{y}$ and $1\leq s,t\leq S$, the ppd reduces to this one
distribution, and no information is lost in (16). But generally the ppd will
not be able to totally summarize all of the sampling distributions which are
plausible under the posterior.

We quantify this loss of information with the Kullback-Leibler (KL)
divergence, which is given for densities $p$ and $q$ over common support as
$D\left(p(\mathbf{y})||q(\mathbf{y})\right)=\mathbb{E}_{p(\mathbf{y})}\log\left(\frac{p(\mathbf{y})}{q(\mathbf{y})}\right)$.
This is a measure of discrepancy between distributions, and with it, we define
a metric for the average discrepancy between distributions
$p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})$ drawn from the posterior
and the ppd.

######  Definition 6 (Posterior Sampling Divergence).

For data $\mathbf{y}$ and model $p(\boldsymbol{\theta},\mathbf{y})$, the
posterior sampling divergence is

$\mathsf{psd}\left(\mathbf{y}\right)=\mathbb{E}_{p(\boldsymbol{\theta}\mid\mathbf{y})}D\left(p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})||p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})\right).$
(17)

Using the Donsker-Varadhan representation and Jensen’s inequality, we get

$\mathsf{psd}(\mathbf{y})\leq\mathbb{E}_{p(\boldsymbol{\theta}\mid\mathbf{y})}\left\\{\sup_{T:\mathcal{Y}\to\mathbb{R}}\left\lvert\mathbb{E}_{p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})}T\left(\mathbf{y}_{\mathrm{rep}}\right)-\mathbb{E}_{p(\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})}T\left(\mathbf{y}_{\mathrm{rep}}\right)\right\rvert\right\\}$

Figure 2: Left: The posterior of $\theta$. Right: The posterior predictive
distribution of $(y_{1},y_{2})$.

In words, the psd lower bounds the degree to which typical sampling
distributions $p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})$ (with
respect to the posterior) and the ppd can be distinguished by a statistic $T$.
A large psd thus indicates increased risk of information loss when using a
ppp-v compared to (15). An example shows how this information loss can be
relevant for practical model evaluation by making it difficult to falsify the
model using a ppp-v. Let $\mathbf{y}=(-10,10)$, with model

$y_{1},y_{2}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathsf{student\text{-}t}\left(\theta,1,10\right),\qquad\theta\sim\mathsf{uniform}\left(-15,15\right),$
(18)

where $\mathsf{student\text{-}t}(\mu,\sigma,d)$ is the t distribution with
location $\mu$, scale $\sigma$, and $d$ degrees of freedom. This results in a
multimodal posterior, plotted in the left panel of Figure 2. The right panel
plots joint samples from the ppd, which is also bimodal despite the
unimodality of the individual sampling distributions.

Consider the test statistics $T_{1}(y_{1},y_{2})=-y_{1}$ and
$T_{2}(y_{1},y_{2})=y_{2}$, and let $p_{T_{1}}(\theta)$ and
$p_{T_{2}}(\theta)$ be the corresponding conditional ppp-vs for
$p(y_{1},y_{2}\mid\theta)$. Figure 3 plots these against $\theta$. The ppp-vs
$p_{T_{1}}$ and $p_{T_{2}}$ are $\approx 0.165$, above usual thresholds for
rejection and thus insufficient for falsification. However, the conditional
$p$-values are vanishingly small over the bulk of the posterior support,
suggesting that the model may be improved by introducing a scale parameter, or
allowing the means to differ for the two observations, for example.

This example points towards a notion of power which does not make reference to
the frequency properties of the ppp-v. Specifically, we will consider a model
assessment to be underpowered if there is additional data (e.g. that contained
in (15)) which would lead us to consider the model fitness deficient (with
respect to the data feature we are testing) despite the particular model
assessment passing (i.e. indicating acceptable compatibility between data
feature and model). Therefore, in light of the above, we view an increasing
psd as increasing the risk of our chosen model assessments suffering power
deficits. In these cases, we have to work harder to find strong evidence for
the falsity of the model (e.g. by examining the conditional $p$-value plots in
Figure 3), and in this sense the model exhibits weaker falsifiability.

Figure 3: Conditional posterior predictive $p$-values for $T_{1}$ (left panel)
and $T_{2}$ (right panel), evaluated at and plotted against posterior draws of
$\theta$.

###  3.1 Posterior Sampling Divergence and Model Expansion

Since $\mathbf{y}$ and $\mathbf{y}_{\mathrm{rep}}$ are conditionally
independent given $\boldsymbol{\theta}$, it follows that

$\mathbb{E}_{p(\mathbf{y})}\mathsf{psd}\left(\mathbf{y}\right)=\mathbf{I}_{p(\boldsymbol{\theta},\mathbf{y},\mathbf{y}_{\mathrm{rep}})}\left(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right).$
(19)

Then for base model $p_{\mathrm{base}}(\boldsymbol{\theta},\mathbf{y})$ and
expansion $p(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y})$, define

$\Delta_{I}=\mathbb{E}_{p(\boldsymbol{\lambda})}\left[\mathbf{I}_{p(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}},\mathbf{y}\mid\boldsymbol{\lambda})}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})-\mathbf{I}_{p(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}},\mathbf{y}\mid\boldsymbol{\lambda}_{0})}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})\right].$

Using the chain rule for cmi, we have that

$\mathbf{I}_{p}\left((\boldsymbol{\theta},\boldsymbol{\lambda}),\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)=\mathbf{I}_{p_{\mathrm{base}}}\left(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)+\Delta_{I}+\mathbf{I}_{p}\left(\boldsymbol{\lambda},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right).$
(20)

As before, the nonnegative term
$\mathbf{I}_{p}\left(\boldsymbol{\lambda},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)$
creates an upward bias for the overall cmi. In some simple examples, the cmi
can be computed exactly:

  1. 1.

For base model $y\mid\theta\sim\mathrm{normal}\left(\theta,1\right)$ and
$\theta\sim\mathrm{normal}\left(0,1\right)$, we get

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=h\left(\boldsymbol{\theta}\mid\mathbf{y}\right)-h\left(\boldsymbol{\theta}\mid\mathbf{y},\mathbf{y}_{\mathrm{rep}}\right)=\log\left(3/2\right)/2.$

Now we add a redundant location parameter, so
$\boldsymbol{\theta}=(\theta_{1},\theta_{2})$, and

$y\mid\boldsymbol{\theta}\sim\mathrm{normal}\left((\theta_{1}+\theta_{2})/\sqrt{2},1\right),\qquad\theta_{1},\theta_{2}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(0,1\right).$

After an invertible reparametrization
$(\mu_{1},\mu_{2})=\phi(\theta_{1},\theta_{2})$, this model is

$y\mid\boldsymbol{\theta}\sim\mathrm{normal}\left(\mu_{1},1\right),\qquad\mu_{1},\mu_{2}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(0,1\right).$

By invariance of the cmi under invertible reparametrization, we have

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=\mathbf{I}(\boldsymbol{\mu},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=h\left(\mu_{1}\mid\mathbf{y}\right)-h\left(\mu_{1}\mid\mathbf{y},\mathbf{y}_{\mathrm{rep}}\right)=\log\left(3/2\right)/2.$

  2. 2.

Now consider a normal location model with data $\mathbf{y}\in\mathbb{R}^{2n}$
for $n\geq 1$:

$\mathbf{y}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\theta,1\right)\text{ for }1\leq i\leq
2n,\qquad\theta\sim\mathrm{normal}\left(0,1\right).$

We expand this model by dividing $\mathbf{y}$ as
$\mathbf{y}=\left(\mathbf{y}^{1},\mathbf{y}^{2}\right)$ with
$\mathbf{y}^{1},\mathbf{y}^{2}\in\mathbb{R}^{n}$ and introducing separate
means $\theta_{1},\theta_{2}$, arriving at:

$\mathbf{y}^{j}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\theta_{j},1\right)\text{ for }1\leq i\leq
n\text{ and }j=1,2,\qquad\theta_{1},\theta_{2}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(0,1\right).$

Now the cmi of the base model is
$\mathrm{CMI}_{\mathrm{base}}(n)=\frac{1}{2}\log\left(\frac{4n+1}{2n+1}\right)$,
whereas the cmi of the expanded model is
$\mathrm{CMI}_{\mathrm{exp}}(n)=\log\left(\frac{2n+1}{n+1}\right)$. Figure 4
plots
$(\mathrm{CMI}_{\mathrm{exp}}(n)-\mathrm{CMI}_{\mathrm{base}}(n))/\mathrm{CMI}_{\mathrm{base}}(n)$
against $n$. Clearly
$\mathrm{CMI}_{\mathrm{exp}}(n)>\mathrm{CMI}_{\mathrm{base}}(n)$ for all $n$,
and $\mathrm{CMI}_{\mathrm{exp}}(n)\to 2\mathrm{CMI}_{\mathrm{base}}(n)$ as
$n\to\infty$.

Figure 4: Percent change in cmi against $1\leq n\leq 30$.

The change in cmi can be separated into two pieces. First, by splitting
$\mathbf{y}$, we reduce the data we have to estimate each of the means
$\theta_{1}$ and $\theta_{2}$. This is reflected in the inequality
$\log\left(\frac{4n+1}{2n+1}\right)>\log\left(\frac{2n+1}{n+1}\right)$. But
parametrizing with two independent means adds a degree of freedom in the
sampling distribution of the expanded model, doubling the constant factor,
which dominates the comparison. However, the latter effect will not always
determine the change in cmi between models, as the next example shows.

  3. 3.

We take the base model from the last example with
$\mathbf{y}\in\mathbb{R}^{n}$ and expand it by adding a precision parameter
and using a jointly normal-gamma prior:

$\mathbf{y}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\theta_{1},\theta_{2}^{-1/2}\right)\text{
for }1\leq i\leq
n,\qquad(\theta_{1},\theta_{2})\sim\mathrm{NG}\left(0,\mu_{\theta_{2}^{-1}},2,\mu_{\theta_{2}^{-1}}\right)$

Here, $\mu_{\theta_{2}^{-1}}>0$ is the prior mean of the variance
$\theta_{2}^{-1}$. The marginal prior on $\theta_{1}$ is
$\mathrm{normal}\left(0,1\right)$, matching the prior in the base model.
Figure 5 shows estimated percentage changes in cmi from base to expanded model
against $n$ for a range of noise levels $r=\mu_{\theta_{2}^{-1}}/n$.

As before, increasing $n$ makes it easier to distinguish sampling
distributions. Similarly, the added degree of freedom introduced by the
precision parameter pushes the cmi larger. Hence, the majority of points in
Figure 5 lie above $0$. However, unlike the last example, the cmi can decrease
in the expanded model if $r$ is sufficiently large. This is because large
values of $r$ create priors that favor sampling distributions with large
scales that are correspondingly harder to distinguish. Nevertheless, the
effect of the added degree of freedom dominates this comparison. For example,
for $n=2$, the noise level in the base model is $r=0.5$. Unless the prior
average noise level in the expanded model is more than double that of the base
model (i.e. $r=1$), we can see that the cmi will increase.

Figure 5: Percentage change in cmi from the base model to the expanded model
against data size $n$ for a range of noise levels $r=\mu_{\theta_{2}^{-1}}$.

  4. 4.

We now vary the prior scale in the model of the first example. Specifically,
take $y\mid\theta\sim\mathrm{normal}\left(\theta,1\right)$ and
$\theta\sim\mathrm{normal}\left(0,\sigma_{p}\right)$. This has cmi given by

$\frac{1}{2}\log\left(\frac{2\sigma^{2}_{p}+1}{\sigma^{2}_{p}+1}\right).$

This is increasing in $\sigma_{p}$, and converges to $0$ as $\sigma_{p}\to 0$
and to $\frac{1}{2}\log(2)$ as $\sigma_{p}\to\infty$. In this case, there is
no change in data size or degrees of freedom in the likelihood, and the cmi
changes only because of the prior.

In these examples we derived simple expressions for the cmi that depended on
sample size, sampling variance, and prior variance. In most cases, these
expressions increased to a finite upper bound in relevant limits, despite the
fact that the cmi is unbounded above in general. The following lower bound in
terms of the Fisher information demonstrates that this self-limiting behavior,
as well as the dominance of the parameter dimension in driving increases the
cmi, is not limited to these simple examples.

######  Theorem 4 (Conditional Mutual Information Lower Bound).

For $M\geq 1$, define the $M$-replicated model:

$p\left(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)},\boldsymbol{\theta}\right)=p(\boldsymbol{\theta})\prod_{i=1}^{M}p\left(\mathbf{y}^{(i)}\mid\boldsymbol{\theta}\right).$

Suppose for $M$ sufficiently large, we have that

  * •

the posterior distributions
$p\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
are normal,

  * •

the observed information matrix of $p(\mathbf{y},\boldsymbol{\theta})$ is
$\gamma$-subexponential for some $\gamma>0$ (i.e. the observed information
does not have heavy tails),

  * •

$\mathbb{E}\lambda_{d}(\boldsymbol{\Sigma})$,
$\mathbb{E}\lambda^{-1}_{1}(\boldsymbol{\Sigma})$,
$\mathbb{E}\lambda^{2}_{d}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)$,
and
$\mathbb{E}\lambda_{1}^{-1}(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})))$
are bounded by some $B>0$ where
$\boldsymbol{\Sigma}=\mathrm{Cov}\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
(i.e. the posterior covariance and Fisher information are neither too small
nor too large on average).

Then for $C$ a constant depending on $\gamma$ and $B$, we have

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})\geq\frac{C}{\log
d}\mathrm{tr}\left(\mathbb{E}_{p(\boldsymbol{\theta},\mathbf{y})}\boldsymbol{\Sigma}^{1/2}_{\mathbf{y}}\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\boldsymbol{\Sigma}_{\mathbf{y}}^{1/2}\right),$
(21)

Remarks:

  * •

The assumption that
$\mathbb{E}\lambda^{-1}_{1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)<B$
rules out singular models with
$\lambda_{1}(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}))=0$. However, such
models can often be reparametrized as
$\boldsymbol{\theta}^{\prime}=\Psi(\boldsymbol{\theta})$ using some
$\Psi:\mathbb{R}^{d}\to\mathbb{R}^{r}$ with $r<d$ such that
$\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}^{\prime})$ is nonsingular.
Applying the result to such a parametrization gives a lower bound for the
original cmi
$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})$.

  * •

Normality of
$p\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
for $M\geq 1$ sufficiently large is almost certainly not satisfied unless it
is satisfied for $M=1$. Nevertheless, if the Bernstein-von-Mises theorem
holds, we would expect
$p\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
to be nearly normal for large $M$ even if it is far from normal for $M=1$. We
thus conjecture a similar bound for more general posteriors.

We note that our bound depends on the parameter dimension through the number
of terms in the trace, and the other influences on the cmi observed in the
above examples enter through the magnitudes of these terms. the self-limiting
phenomenon can be seen in this bound through this multiplication of covariance
and information matrices. For instance, in a Bernstein-von-Mises type limit,
the posterior concentrates around the true parameter
$\boldsymbol{\theta}_{0}$, we get
$\boldsymbol{\Sigma}_{\mathbf{y}}\approx\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}_{0})$,
and the bound becomes $Cd/\log d$, which is independent of all non-dimensional
factors.

The dependence on the model dimension in this bound comes directly through the
dimensions of the covariance and information matrices. The sampling and prior
effects appear through the individual eigenvalues of the posterior-normalized
Fisher information matrix. We also note that the self-limiting phenomenon can
be seen in this bound through this multiplication of covariance and
information matrices. For instance, in a Bernstein-von-Mises type limit, the
posterior concentrates around the true parameter $\boldsymbol{\theta}_{0}$, we
get
$\boldsymbol{\Sigma}_{\mathbf{y}}\approx\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}_{0})$,
and the bound becomes $Cd/\log d$.

##  4 Marginal Entropy and Sampling Divergence

Let $p_{\mathrm{base}}(\mathbf{y},\boldsymbol{\theta})$ and
$p(\mathbf{y},\boldsymbol{\theta},\boldsymbol{\lambda})$ be a base model and
expansion. Consider the following two extreme scenarios:

  * •

Let $p_{\mathrm{base}}(\mathbf{y}\mid\boldsymbol{\theta})=q(\mathbf{y})$ for a
density $q$. Then the likelihood is constant, and
$h_{p_{\mathrm{base}}}(\boldsymbol{\theta})=h_{p_{\mathrm{base}}}(\boldsymbol{\theta}\mid\mathbf{y},\mathbf{y}_{\mathrm{rep}})$
for all $\mathbf{y}$ and $\mathbf{y}_{\mathrm{rep}}$. It follows immediately
that
$\mathbf{I}(\boldsymbol{\theta},\mathbf{y})=\mathbf{I}(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y})=0$,
and so any nontrivial expansion of $p_{\mathrm{base}}$ with
$\boldsymbol{\theta}$ and $\boldsymbol{\lambda}$ independent must decrease the
marginal posterior entropy of $\boldsymbol{\theta}$ and increase the cmi.

  * •

Let
$p_{\mathrm{base}}(\mathbf{y}\mid\theta)=\mathrm{normal}(\mathbf{y}\mid\theta,1)$,
and then take the expansion
$p(\mathbf{y}\mid\theta,\lambda)=p_{\mathrm{base}}(\mathbf{y}\mid\theta+\lambda)$
with priors $\theta\sim\mathrm{normal}(0,\sigma_{\theta}^{2})$ and
$\lambda\sim\mathrm{normal}\left(0,\sigma_{\lambda}^{2}\right)$. We then have
for the expanded model that

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y})=\frac{1}{2}\log\left(1+\frac{\sigma^{2}_{\theta}}{1+\sigma^{2}_{\lambda}}\right),\quad\mathbf{I}\left((\boldsymbol{\theta},\boldsymbol{\lambda}),\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)=\frac{1}{2}\log\left(\frac{1+2(\sigma^{2}_{\theta}+\sigma^{2}_{\lambda})}{1+\sigma^{2}_{\theta}+\sigma^{2}_{\lambda}}\right).$

Similarly, the cmi in the base model is

$\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}_{\mathrm{rep}}\mid\mathbf{y}\right)=\frac{1}{2}\log\left(\frac{1+2\sigma^{2}_{\theta}}{1+\sigma^{2}_{\theta}}\right).$

Both cmi expressions are bounded above by $\frac{1}{2}\log(2)$. Furthermore,
by taking $\sigma^{2}_{\theta}$ large, we can ensure that the cmi for the base
model is arbitrarily close to this limit. Then the cmi can increase only
negligibly in the expanded model regardless of $\sigma_{\lambda}^{2}$.
However, fixing $\sigma^{2}_{\theta}$, the mi in the expanded model tends to
$0$ as $\sigma^{2}_{\lambda}\to\infty$.

In each case, one of the models has singular Fisher information, and either
the mi or cmi degrades (i.e. decreases or increases respectively) while the
other changes negligibly. When both the base model and expanded model have
nonsingular Fisher information, we can show that in certain cases there is a
strict trade-off asymptotically between worsening (i.e. decreasing) mi and
worsening (i.e. increasing) cmi.

######  Theorem 5.

Let $p\left(\boldsymbol{\theta},\boldsymbol{\lambda},\mathbf{y}\right)$ be an
expansion of $p_{\mathrm{base}}(\boldsymbol{\theta},\mathbf{y})$, and let
$\\{\iota_{i}\\}_{i=1}^{d}$ and $\\{\iota^{\mathrm{cond}}_{i}\\}_{i=1}^{d}$ be
the eigenvalues of
$\mathbb{E}\boldsymbol{\mathcal{I}}_{\mathrm{base}}(\boldsymbol{\theta})$ and
$\mathbb{E}\boldsymbol{\mathcal{I}}(\boldsymbol{\theta}\mid\boldsymbol{\lambda})$
respectively.

Furthermore, suppose that the conditions of Theorems 2, 3, and 4 hold. In
particular this requires that:

  1. 1.

the marginal priors $p_{\mathrm{base}}(\boldsymbol{\theta})$ and
$p(\boldsymbol{\theta})$ are log concave.

  2. 2.

the posteriors $p_{\mathrm{base}}(\boldsymbol{\theta}\mid\mathbf{y}_{M})$ and
$p(\boldsymbol{\theta},\boldsymbol{\lambda}\mid\mathbf{y}_{M})$ are normal for
all $\mathbf{y}_{M}$ for $M$ large enough, where $\mathbf{y}_{M}$ is a vector
of $M$ i.i.d. replicated draws from
$p_{\mathrm{base}}(\mathbf{y}\mid\boldsymbol{\theta})$ and
$p(\mathbf{y}\mid\boldsymbol{\theta},\boldsymbol{\lambda})$ respectively,

  3. 3.

the expected spectra of the Fisher information matrices
$\boldsymbol{\mathcal{I}}_{\mathrm{base}}(\boldsymbol{\theta})$ and
$\boldsymbol{\mathcal{I}}(\boldsymbol{\theta},\boldsymbol{\lambda})$ are
bounded above and below by universal constants, and similarly for the expected
spectra of the posterior covariance matrices,

  4. 4.

and a few other regularity conditions on the smoothness of the relevant
densities and tails of the observed information.

Additionally, we further assume that the distributions of the Fisher
information matrices and posterior covariance matrices are not too skewed in
the sense of Lemma 10 in Appendix F. This ensures that the (random) Fisher
information and covariance matrices are sufficiently well-summarized by their
means.

Under these conditions, for increasing functions $\psi_{1}$ and $\psi_{2}$, we
have the following inequalities:

$\mathbf{I}_{\mathrm{base}}\left(\mathbf{y},\boldsymbol{\theta}\right)\leq\psi_{1}\left(\sum_{j=1}^{d}\iota_{j}\right),\qquad\mathbf{I}_{\mathrm{base}}\left(\mathbf{y}_{\mathrm{rep}},\boldsymbol{\theta}\mid\mathbf{y}\right)\geq\frac{1}{\log
d}\sum_{i=j}^{d}\psi_{2}\left(\iota_{j}\right)$ (22)
$\displaystyle\mathbf{I}\left(\mathbf{y},\boldsymbol{\theta}\right)\leq\psi_{1}\left(\sum_{j=1}^{d}\iota^{\mathrm{cond}}_{j}\right)-\Delta_{i},$
(23)
$\displaystyle\mathbf{I}\left(\mathbf{y}_{\mathrm{rep}},(\boldsymbol{\theta},\boldsymbol{\lambda})\mid\mathbf{y}\right)\geq\frac{1}{\log
d^{\mathrm{exp}}}\left[\sum_{j=1}^{d}\psi_{2}\left(\iota^{\mathrm{cond}}_{j}\right)+\Delta_{f}\right],$
(24)

where $\Delta_{f}\geq 0$, and $\Delta_{i}\geq 0$ so long as knowledge of
$\mathbf{y}$ does not decrease the information that $\boldsymbol{\lambda}$
provides about $\boldsymbol{\theta}$ in the sense that
$\sum_{j=1}^{d}\Delta_{j}\leq 0$ (where the $\Delta_{j}$ are defined in
Theorem 3). Furthermore, if $p$ is a totally diluting expansion of
$p_{\mathrm{base}}$, then we have

$\psi_{i}\left(\sum_{j=1}^{d}\iota^{\mathrm{cond}}_{j}\right)\leq\psi_{i}\left(\sum_{j=1}^{d}\iota_{j}\right),$
(25)

and if $p$ is a totally nondiluting expansion of $p$, then we have

$\sum_{j=1}^{d}\psi_{f}\left(\iota^{\mathrm{cond}}_{j}\right)\geq\sum_{j=1}^{d}\psi_{f}\left(\iota_{j}\right).$
(26)

###### Proof.

See Appendix F. ∎

This result substantially generalizes the pattern we observed in our
introductory regression example, where predictor correlation structures that
offered better identification created greater posterior uncertainty about the
sampling distribution and vice versa. We can interpret this theorem in a
positive and negative light. Negatively, when our prior information is
relatively unstructured, the process of model expansion may force us to
confront either weakly identified marginal inferences or a large posterior
sampling divergence. Positively, a prior with sufficient dependence between
the parameters may allow us to avoid these difficulties, even if this prior is
weak in the sense of carrying relatively little marginal information about any
particular parameter. In the next sections, we explore examples where these
two difficulties occur and demonstrate how they can be partially overcome with
sufficiently rich posterior summaries.

##  5 Example: Inference Under Poor Identification

The next two sections explore methodological implications of the above results
with concrete examples. Here we consider an example in which we get ‘stuck’
between an implausibly simple base model and an expanded model where the
parameter of interest is too weakly identified to allow strong conclusions. We
then show how the expanded model can still support nontrivial inferences
which, in this case, inform how we should collect future data. In so doing, we
hope to demonstrate (i) how the concerns raised in Sections 2 and 4 motivate
looking beyond standard marginal posterior summaries, and (ii) that weak
identification need not be an inferential dead end for a statistical analysis.

###  5.1 Two Models for Grouped Data

We first define three simulated data sets, each generated by drawing random
samples of $M$ measurements from each of $L$ subpopulations. We take $M=2$ and
$L=20$ for all data sets, but vary the ratio of within- and between-
subpopulation variances between them. Appendix B provides a complete
description of the data generating process. Figure 6 plots the data, with a
row for each subpopulation, a column for each data set, and dots for the
individual measurements.

Figure 6: Columns: the three data sets. Rows: the 20 subpopulations. Cells:
the two data points drawn from each subpopulation, connected with a horizontal
line to show their range.

The unobserved grand (i.e. superpopulation-level) mean will be our quantity of
interest. We also take more positive values of the grand mean to represent
“better” outcomes (i.e. more desirable from the researcher’s standpoint). We
think about the identification of the grand mean in two ways:

  1. 1.

In terms of the reduction in entropy from marginal prior to posterior or the
ratio of their standard deviations (which is just a monotonic transform of the
former when both distributions are normal). This notion is general, but being
unit-free, it is somewhat unnatural for drawing practical conclusions.

  2. 2.

The posterior probability that the grand mean is below a threshold of
practical significance. Specifically, we consider the grand mean to be
practically significant only if it exceeds $1$.

While these notions are distinct, greater entropy will tend to be associated
with larger probability of a not-practically-significant effect (so long as
the prior is sufficiently constraining in the extremes).

####  5.1.1 An Oversimplified Initial Model

In our base model, we assume the subpopulation-level distributions are
identical (so the subpopulation means equal the grand mean). Letting
$\mathbf{y}_{ml}$ be the data with $1\leq m\leq M$ indexing measurements and
$1\leq l\leq L$ indexing subpopulations, the sampling distribution is
$\mathbf{y}_{ml}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\mu,\sigma_{*}\right)$.

Here $\mu$ represents the grand mean, to which we assign prior
$\mathrm{normal}\left(0,\left\lvert\mu_{0}\right\rvert\right)$, where
$\mu_{0}$ is the “true” value used in simulating the three data sets. We treat
$\sigma_{*}$ as a hyperparameter which we set to a prior guess.

Figure 7 shows histograms of draws from the posteriors for each data set along
with the prior density. For each data set, the posterior standard deviation is
approximately a quarter of the prior standard deviation, and the probabilities
of a practically insignificant effect (i.e. $\mu<1$) are between $2\%$ and
$10\%$. Overall, we have improved on our prior knowledge, but the posterior
cannot fully rule out the possibility of an insignificant effect.

Figure 7: Histograms of samples from the posterior distributions of $\mu$
under the base model fit to each data set. The red curve shows the density of
the prior on $\mu$.

####  5.1.2 A More Plausible Expanded Model

Two features of the base model stand out for criticism:

  1. 1.

We usually cannot confidently guess the scale $\sigma_{*}$ with just prior
information.

  2. 2.

If we regard the subpopulations as distinct for data collection purposes, then
we likely have reason to believe their distributions could be distinct.

Our expanded model thus adds a parameter for the subpopulation scale and
allows the subpopulations to have distinct means (drawn from some
superpopulation). In symbols,

$\begin{gathered}\mathbf{y}_{ml}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\theta_{l},\sigma\right),\qquad\theta_{l}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\mu,\tau\right),\\\
\mu\sim\mathrm{normal}\left(0,\lvert\mu_{0}\rvert\right),\qquad\sigma\sim\mathrm{gamma}\left(a_{\sigma},b_{\sigma}\right),\qquad\tau\sim\mathrm{gamma}\left(a_{\tau},b_{\tau}\right).\end{gathered}$
(27)

In the above, $a_{\sigma},b_{\sigma},a_{\tau},b_{\tau}$ are hyperparameters of
the model, and $\mu_{0}$ is unchanged from the base model. The hyperparameters
$a_{\sigma}$ and $b_{\sigma}$ are taken such that the prior mode for $\sigma$
equals our prior guess $\sigma_{*}$. The prior on $\mu$ is unchanged from the
base model, and $\mu$ again corresponds to the grand mean. Setting
$\sigma=\sigma_{*}$ and $\tau=0$ recovers the base model as a special case.

Figure 8 shows the posterior inferences for $\mu$ in the expanded model. In
the first and third data sets, the posterior standard deviation is $\approx
1/3$ larger than in the base model. For the second data set, it is more than
twice as large as in the base model. In all cases, the standard deviation is
at least a third that of the prior. Similarly, the probabilities of a
practically insignificant effect are all at least three times larger than in
the base model, ranging from $\approx 6\%$ in the best case to over $30\%$.

Figure 8: Histograms of samples from the marginal posteriors of $\mu$ under
the expanded model fit to each data set. The red curve shows the density of
the prior on $\mu$.

We might now return to the base model or focus on subpopulation means to get
better identification. However, the former is inadvisable since the base model
was implausible, and the latter will likely fail since the subpopulation
sample sizes are tiny. A less convenient but much safer approach is simply to
gather more data. In fact, inferences from the expanded model can yield
relevant information for future data collection despite the weak
identification.

###  5.2 A Bootstrap Comparison of Two Sampling Schemes

It has been recognized that bootstrapping new datasets, fitting a Bayesian
model to each, and then aggregating the results can provide more information
than the fit to the observed data alone [6]. This idea has been used, e.g.,
for model criticism and selection tasks [14, 13]. Here, we demonstrate that a
similar procedure with a Bayesian parametric bootstrap of future data can help
us learn about how to sample such data despite the underwhelming marginal
identification of the posterior.

First, we must identify the candidate sampling schemes. We may sample the same
subpopulations that we sampled originally, or we may sample from new
subpopulations, or some combination of these. If the within-subpopulation
variance is low, then sampling the same subpopulations may yield little
improvement for identification. If instead the between-subpopulation variance
is low, sampling the same subpopulations may provide enough information about
the individual $\theta_{l}$ to identify $\mu$ well. The issue may be further
complicated by cost differences. For our purposes, we assume that a sample
from a new subpopulation is four times the cost of a sample from an existing
subpopulation.

To simplify our comparison, we focus on whether it is more efficient to sample
exclusively from existing subpopulations or from new subpopulations. To answer
this, we use the joint posterior over parameters and replicated data to
simulate sampling new data with these two schemes. We then compare posteriors
on $\mu$ after refitting the model to these enlarged data sets. Replicating
many times then allows us to assess the relative risks of each approach.
Figure 9 diagrams these resampling schemes with full pseudocode in Appendix B.

Figure 9: Schematic representation of our two resampling schemes. Dashed
arrows represent sampling the output variable conditional on the input. Solid
arrows represent that the output variable is formed by evaluating the circled
function on the input variables.

We take $M^{\mathrm{new}}=8$ additional samples from each of the $L=20$
existing subpopulations in the first scheme and $M=2$ samples from
$L^{\mathrm{new}}=20$ new subpopulations in the second scheme. These are
performed at equivalent overall cost, hence the difference in total sample
size. We replicate each scheme $R=500$ times. We also repeat this process with
the prior replacing the posterior at each step in order to demonstrate that
the posterior inferences differ from prior inferences despite the marginal
weak identification.

For each of the $r=1,\ldots,R$ replications and each method, define the ratio
$\rho_{r}=\sigma^{(r)}/\sigma_{\mathrm{obs}}$, where $\sigma^{(r)}$ is the
posterior standard deviation of $\mu$ for the $r^{\mathrm{th}}$ expanded data
set simulated from the given method, and where $\sigma_{\mathrm{obs}}$ is the
posterior standard deviation of $\mu$ given just the observed data set. Table
10 presents the averages of these ratios over all replications
$\overline{\rho}=\frac{1}{R}\sum_{r=1}^{R}\rho_{r}$ for each method and each
of our observed data sets.

Despite weak identification, the posterior and prior columns differ,
indicating that our simulated data reflects information learned from the
observed data. Furthermore, the average improvement to identification for each
scheme depends substantially on the data set, with each scheme winning in one
data set and tying in the third.

Figure 11 gives a finer-grained comparison, plotting histograms of the
$\rho_{r}$ for the two schemes. While the first data set yields a tie in the
comparison of averages, the distribution of $\rho_{r}$ is wider when sampling
the same subpopulations than when sampling new subpopulations. This suggests a
risk-reward trade-off: the latter scheme gives a more predictable reduction in
uncertainty while the former carries the possibility of a greater reduction.
Together this demonstrates that we can get nontrivial inferences which depend
strongly on the particular data observed despite weak identification.

|  Same subpops |  New subpops |  Same subpops (prior) |  New subpops (prior)
---|---|---|---|---
Data set #1 | 0.74 | 0.74 | 1.18 | 1.04
Data set #2 | 0.69 | 0.8 | 0.74 | 0.77
Data set #3 | 0.94 | 0.71 | 1.16 | 1.03
Figure 10:  Figure 11: Histograms of estimated posterior standard deviations
of $\mu$ for 500 simulations of future data under each sampling scheme and
each data set.

##  6 Example: Conditional Model Checking

We next consider a model where the conditional ppp-vs contain substantially
more information than the marginal $p$-value (i.e. the regular ppp-v) for a
relevant test statistic. This conditional information motivates a modification
which further improves model fitness and resolves a problem that was masked by
the marginal $p$-value.

We take as our base model the election forecasting model of [11]. We pay
particular attention to herding, a process in which pollsters systematically
augment their raw data in such a way that their published poll numbers are
closer to the existing consensus of recent polls than they would otherwise be.
That the variance among presidential election polls fell well below the
minimal expected variance (assuming independence between polls) just before
the election has been used as evidence for the existence of herding in 2016
[22]. As this model does not explicitly account for herding, this will be
starting point of our model checking.

###  6.1 The Base Forecasting Model and a Check for Herding

A more complete discussion of the model specification is given in Appendix C,
and full details can be found in [11]. The primary purpose of the model is to
infer the level of support for the Democratic candidate over time and across
states. This level of support is represented by a matrix parameter
$\boldsymbol{\mu}\in\mathbb{R}^{S\times T}$ with rows representing the $S=51$
states (including Washington DC) and columns representing the $T$ days from
the start of measurement until election day. This parameter is assigned a time
series prior:

$\boldsymbol{\mu}_{t}\mid\boldsymbol{\mu}_{t+1}\sim\mathrm{normal}\left(\boldsymbol{\mu}_{t+1},\boldsymbol{\Sigma}^{\mu}\right)\text{
for }1\leq t\leq T-1,\text{ and
}\boldsymbol{\mu}_{T}\sim\mathrm{normal}\left(\mathbf{m}^{\mathrm{f}},\mathbf{S}^{\mathrm{f}}\right).$
(28)

Here $\boldsymbol{\Sigma}^{\mu}\in\mathbb{R}^{S\times S}$ is a hyperparameter
encoding correlation between states and variation over time, constructed using
demographic data, polling from previous elections, and domain knowledge.
Likewise, $\mathbf{m}^{f}\in\mathbb{R}^{S}$ and
$\mathbf{S}^{f}\in\mathbb{R}^{S\times S}$ are hyperparameters set using a
‘fundamentals forecast’ derived from variables known in political science to
be good predictors of U.S. election outcomes.

Results of state and national polls are modeled with a binomial distribution
that combines $\boldsymbol{\mu}$ with terms representing sources of polling
bias. Letting $i=1,\ldots,N_{\mathrm{state}}$ index state polls and $y_{i}$
denote the number of respondents supporting the Democratic candidate out of
$n_{i}$ respondents in poll $i$, we have

$y_{i}\sim\mathrm{binomial}\left(\mathsf{logit}^{-1}\left(\boldsymbol{\mu}_{s_{i},t_{i}}+\beta_{i}\right),n_{i}\right),$
(29)

where $s_{i}$, $t_{i}$ denote the state and day for poll $i$ and $\beta_{i}$
models various sources of bias (see Appendix C for further discussion).

National polls are modeled similarly, except state-level terms including
$\boldsymbol{\mu}_{st}$ are averaged with weights accounting for each state’s
share of the national vote in the previous election.

Since most national polls lie in $(0.4,0.6)$ in any close race, the binomial
sampling model (29) effectively places a lower bound of $\sqrt{0.6\times
0.4/n_{i}}\approx 0.49/\sqrt{n_{i}}$ on the standard deviation of national
poll $i$. Thus, this model will be incompatible with sufficiently low poll
variance over any interval with enough published polls. Figure 12 compares the
observed variation in the last ten days with posterior predictive replications
and shows that the observed polling variance is substantially smaller than in
posterior simulations.

Figure 12: The observed data is highlighted in red. Left: box plots of the
last tend days of polls for the observed data and fifty posterior predictive
replications, ordered by range. Right: histogram of standard deviations of the
last ten days of polls for the observed data and 6000 posterior predictive
replications.

###  6.2 A Simple Model of Herding

By allowing dependence among polls, adding a herding mechanism may resolve
this tension between model and data. First define
$\theta^{\mathrm{state}}_{i}=\mathsf{logit}^{-1}\left(\boldsymbol{\mu}_{s_{i},t_{i}}+\beta_{i}\right)$
and
$\theta^{\mathrm{nat}}_{j}=\mathsf{logit}^{-1}\left(\overline{\mu}_{t_{j}}+\overline{\beta}_{j}\right)$
for $1\leq i\leq N^{\mathrm{state}}$ and $1\leq j\leq N^{\mathrm{nat}}$, where
$\overline{\mu}_{t_{j}}$ and $\overline{\beta}_{j}$ are weighted averages of
the means and biases over the states. The $\theta$ represent the expected
average support for the Democratic candidate in random samples drawn from the
sampling frame of the corresponding poll. We model the sampling process of
national poll $j$ as follows.

  1. 1.

The pollster samples their sampling frame, observing unherded result:

$y_{j}^{\mathrm{nat}}\mid\theta^{\mathrm{nat}}_{j}\sim\mathrm{binomial}\left(\theta_{j}^{\mathrm{nat}},n_{j}^{\mathrm{nat}}\right).$

  2. 2.

The pollster calculates a herding target $\mu^{\mathrm{herd}}_{j}$
representing their quantification of the consensus of polls to which they
compare their $y^{\mathrm{nat}}_{j}$.

  3. 3.

The pollster herds $y_{j}^{\mathrm{nat}}$ by some fraction
$\lambda_{j}^{\mathrm{herd}}\in(0,1)$ towards $\mu^{\mathrm{herd}}_{j}$. They
publish this herded result:

$p_{j}^{\mathrm{nat}}=\left(1-\lambda^{\mathrm{herd}}_{j}\right)\frac{y_{j}^{\mathrm{nat}}}{n^{\mathrm{nat}}_{j}}+\lambda^{\mathrm{nat}}_{j}\mu_{j}^{\mathrm{herd}}.$

To simplify the implementation of this scheme, we use the normal approximation
to the binomial, arriving at

$\displaystyle\frac{y_{j}^{\mathrm{nat}}}{n^{\mathrm{nat}}_{j}}\;\Bigg{|}\;\theta^{\mathrm{nat}}_{j}$
$\displaystyle\sim\mathrm{normal}\left(\theta_{j}^{\mathrm{nat}},\sqrt{\theta_{j}^{\mathrm{nat}}\left(1-\theta_{j}^{\mathrm{nat}}\right)\Big{/}n^{\mathrm{nat}}_{j}}\right),$
$\displaystyle
p_{j}^{\mathrm{nat}}\mid\theta_{j}^{\mathrm{nat}},\mu^{\mathrm{herd}}_{j},\lambda_{j}^{\mathrm{herd}}$
$\displaystyle\sim\mathrm{normal}\Big{(}\left(1-\lambda^{\mathrm{herd}}_{j}\right)\theta^{\mathrm{nat}}_{j}+\lambda^{\mathrm{nat}}_{j}\mu_{j}^{\mathrm{herd}},$
(30) $\displaystyle\hskip
59.75095pt\left(1-\lambda^{\mathrm{herd}}_{j}\right)\sqrt{\theta_{j}^{\mathrm{nat}}\left(1-\theta_{j}^{\mathrm{nat}}\right)\Big{/}n^{\mathrm{nat}}_{j}}\Big{)}.$
(31)

The herding model for the $1\leq j\leq N^{\mathrm{state}}$ polls of states
$1\leq s\leq 51$ uses analogous targets $\mu^{\mathrm{herd}}_{s_{j},j}$ and
percentages $\lambda^{\mathrm{herd}}_{s_{j},j}$, where $s_{j}$ is the state in
which poll $j$ was conducted. Now let $m^{\mathrm{nat}}$, $s^{\mathrm{nat}}$,
$m^{\mathrm{state}}_{s}$ and $s^{\mathrm{state}}_{s}$ be the sample averages
and standard deviations of the national polls and polls of state $1\leq s\leq
51$ respectively. We set the following priors for these parameters.

$\mu^{\mathrm{herd}}_{i}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(m^{\mathrm{nat}},\frac{2}{3}s^{\mathrm{nat}}\right),\quad\mu^{\mathrm{herd}}_{s_{j},j}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(m^{\mathrm{state}}_{s_{j}},\frac{2}{3}s^{\mathrm{state}}_{s_{j}}\right).$
(32)

These essentially just constrain the herding targets to lie in $\pm 2$
standard deviations of the overall polling means, keeping them away from the
extremes of our data. Next let $\mathcal{C}=\\{1\leq i\leq
N^{\mathrm{nat}}\mid t_{i}\geq T-10\\}$ be the indices of national polls
conducted in the last ten days prior to the election. We then set priors:

$\displaystyle\lambda^{\mathrm{nat}}_{i},\lambda^{\mathrm{state}}_{j,s_{j}}$
$\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{Beta}\left(1,9\right),\text{ for all }1\leq j\leq
N^{\mathrm{state}}\text{ and }i\in\mathcal{C}^{c},$
$\displaystyle\lambda^{\mathrm{nat}}_{i}\Bigg{|}\mu^{\mathrm{last}}_{\lambda},k^{\mathrm{last}}_{\lambda}$
$\displaystyle\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{Beta}\left(\mu^{\mathrm{last}}_{\lambda}k^{\mathrm{last}}_{\lambda},\left(1-\mu^{\mathrm{last}}_{\lambda}\right)k^{\mathrm{last}}_{\lambda}\right),\text{
for all }i\in\mathcal{C},$ (33)

where we use hyperpriors
$\mu^{\mathrm{last}}_{\lambda}\sim\mathrm{uniform}\left[0,1\right]$ and
$k^{\mathrm{last}}_{\lambda}\sim\mathrm{normal}\left(200,120\right)$. We also
constrain the $\lambda$ explicitly to $[0,0.9]$ since values too close to $1$
create adverse geometry that frustrates the sampler, and since our prior
beliefs rules out such values. The hierarchical prior on the last ten days of
national herding parameters serves two purposes. Since the polls in this
period are our primary evidence for herding, inference for
$\mu^{\mathrm{last}}_{\lambda}$ and $k^{\mathrm{last}}_{\lambda}$ is of
substantive interest. Also, since data is especially dense in this period, we
can better estimate the herding parameters for these polls (and can thus
afford the weaker marginal priors implied by the hierarchical structure).

Figure 13 displays the posterior predictive check of polling variation in the
last ten days for the expanded model. While the observed variation is still
relatively small, it is no longer implausible. We may now be tempted to stop
and declare our model good enough, as it is not obvious what further
improvements to make. We can, however, extract more information than is
revealed by this (marginal) posterior predictive check.

Figure 13: The observed data is highlighted in red. Left: box plots of the
last tend days of polls for the observed data and fifty posterior predictive
replications, ordered by total range. Right: histogram of the standard
deviation of the last ten days of national polls for the observed data and six
thousand posterior predictive replications.

###  6.3 A More Informative Conditional Model Check

Recall that the conditional posterior predictive $p$-values (cppp-vs) are
defined for test statistic $T$ as
$p_{T}(\boldsymbol{\theta})=\mathbb{E}_{p(\mathbf{y}_{\mathrm{rep}}\mid\boldsymbol{\theta})}\mathbbm{1}\left\\{T(\mathbf{y}_{\mathrm{rep}})\geq
T(\mathbf{y})\right\\}$. The top panel of Figure 14 plots the cppp-vs with $T$
equal to the standard deviation of the last ten days of national polls against
posterior samples of $\mu^{\mathrm{last}}_{\lambda}$, the average herding
percentage in the last ten days of polling. The marginal distributions
displayed on the axes show that the distribution of cppp-vs has a heavy right
tail with the preponderance of sampled $p$-values less than the average of
$0.05$. Furthermore, the fit of the model appears much better for larger
values of $\mu^{\mathrm{last}}_{\lambda}$, suggesting that model fitness may
improve if the posterior favored more herding. But it is unclear how we should
achieve this as using a more informative prior to favor higher values of
$\mu^{\mathrm{last}}_{\lambda}$ would be inconsistent with our (lack of) prior
knowledge.

Figure 14: Conditional cppp-vs for the standard deviation of the last ten days
of national polls in the expanded model plotted against the population average
herding percentage over that time frame (top) and the standard deviation of
the herding targets over that time frame (bottom). Black lines indicate the
marginal ppp-v. Estimated marginal distributions are displayed on the margins.

Another comparison tells us more. The bottom panel of Figure 14 plots the same
cccp-vs against the standard deviation of the herding targets
$\mu^{\mathrm{herd}}_{i}$ in the last ten days before the election. Less
variability in the herding targets is associated with better fit, indicating
another avenue for improvement.

By leveraging only the gross features of the poll results, the prior on
$\mu^{\mathrm{herd}}_{i}$ was designed to be weakly informative. However, this
prior only enforces that $\mu^{\mathrm{herd}}_{i}$ not be extreme compared to
all polls in the series. But clearly $\mu^{\mathrm{herd}}_{i}$ should also be
non-extreme compared to polls around the specific time that poll $i$ was
conducted. Including this prior information should reduce the posterior
standard deviation of herding targets and thus hopefully improve the fit.

###  6.4 Using the Conditional Check to Improve the Model

To construct a new prior on the targets $\mu^{\mathrm{herd}}_{i}$, we
calculate the trailing averages

$c^{i}_{d}=\frac{1}{\lvert\mathcal{C}_{t_{i}}^{d}\rvert}\sum_{k\in\mathcal{C}^{i}_{d}}\frac{y_{k}^{\mathrm{nat}}}{n_{k}^{\mathrm{nat}}}$
(34)

for $d=2,\ldots,7$ days before each national poll $i$, where
$\mathcal{C}^{d}_{t}=\left\\{1\leq i\leq N^{\mathrm{nat}}\mid t-d\leq
t_{i}\leq t\right\\}$.

We then define the average and spread of these:

$m_{c}^{i}=\frac{1}{6}\sum_{d=2}^{7}c^{i}_{d},\text{ and
}s_{c}^{i}=\sqrt{\frac{1}{6}\sum_{d=2}^{7}\left(c^{i}_{d}-m_{c}^{i}\right)^{2}}+0.0002.$
(35)

The addition of a small constant $2^{-4}$ ensures the spread cannot shrink to
$0$ when data is sparse. We then define the improved prior
$\mu_{i}^{\mathrm{herd}}\sim\mathrm{normal}\left(m_{c}^{i},s_{c}^{i}\right)$,
with priors for the state herding targets constructed similarly from the
state-level time series. This prior now encodes the idea that
$\mu_{i}^{\mathrm{herd}}$ should look like a consensus of recent polls.

Figure 15: Posterior samples under the first expanded model of the average
herding percentage for national polls against the standard deviation of
herding targets for national polls both over the last ten days of polls.
Estimated marginals are given on the margins.

Figure 16 displays the same pair of conditional posterior predictive checks
using this more informative prior. The (marginal) ppp-v is nearly $0.3$, and
the unpleasant ‘spike and slab’ shape of the cppp-vs that appeared in Figure
14 has been attenuated. Furthermore, we no longer find intervals of large
posterior probability where the cppp-vs are vanishingly small in either plot.
Thus, the marginal check looks good, and the conditional check no longer
indicates obvious directions for improvement.

The improvement in cppp-vs versus the standard deviation of
$\mu_{i}^{\mathrm{herd}}$ is expected, but the improvement in the comparison
with the mean herding percentage (achieved by higher mean herding in the
posterior) may be somewhat surprising. The source of the improvement is
revealed in Figure 15, which displays a strong negative association between
the standard deviations of the $\mu_{i}^{\mathrm{herd}}$ and the
$\mu^{\mathrm{last}}_{\lambda}$ in the last ten days for the first expanded
model. Thus, solving the problem for the herding targets also solved the
herding percentage problem for free.

Figure 16: cppp-vs for the standard deviation of national polls under the
second expanded model, plotted against the population average herding
percentage (top) and the standard deviation of the herding targets (bottom),
all computed over the last ten days of polls. Black lines indicate marginal
ppp-vs. Estimated marginals are given on the margins.

This example shows that the concern of Section 3 that marginalized model
checks like the ppp-v could obscure information useful for assessing model
fitness is more than purely theoretical. More positively, it also clearly
demonstrates how the cppp-v can be a powerful tool for motivating specific
model improvements, especially when the marginal ppp-v is neither implausibly
small nor reassuringly large.

##  7 Conclusions

When constructing a model for a given data analysis, a statistician should
balance various desiderata, including:

  * •

model predictions compatible with what is known about the world;

  * •

inferences sufficiently well identified to support nontrivial conclusions;

  * •

model checks powerful enough to reveal frictions between model and data.

When model checks reveal deficiencies, the first item is no longer satisfied,
and a better model should be sought. In practice, this is often an expansion
of the previous model. If such an expansions are not accompanied by
sufficiently strong prior information (in the form of prior dependence of
parameters, not the marginal scales), then our results demonstrate that a
tension may easily arise within these three goals as the model dimension
grows. Insofar as the first desideratum is most essential, this motivates
methods that can extract useful information even when one of the latter two
desiderata are not satisfied. One avenue is to pursue richer inferential
summaries. As demonstrated by the last two examples, the full posterior often
contains such rich inferential data. This data is both accessible (by looking
beyond e.g. marginal means, standard deviations, and $p$-values) and capable
of supporting nontrivial conclusions (which are concealed by the common
marginal summaries).

Many directions for future work remain. Determining whether the trade-off
observed in Lemma 11, which depended on Fisher information-based bounds, could
be strengthened (e.g. to a relation between information theoretic quantities
directly) would be of particular interest. We would also like to have better
tools for extracting joint inferences from the posterior. For example, the
posterior predictive resampling we performed in Section 5 could become
prohibitively expensive in large models, and methods to rapidly approximate
such results would be useful. More broadly, there is at the time of writing no
canonical method for estimating conditional analogs of the usual marginal
summaries (e.g. means and variances), which would aid the study of joint
inferences in practice.

## References

  * [1] Efe Aras, Kuan-Yun Lee, Ashwin Pananjady, and Thomas A. Courtade.  A family of bayesian cramér-rao bounds, and consequences for log-concave priors.  CoRR, abs/1902.08582, 2019.
  * [2] M. J. Bayarri and James O. Berger.  P values for composite null models.  Journal of the American Statistical Association, 95(452):1127–1142, 2000.
  * [3] M.J. Bayarri and J.O. Berger.  Quantifying surprise in data and model verification.  In J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith, editors, Bayesian Statistics 6, pages 53–82. Oxford University Press, 1999\.
  * [4] David M. Blei.  Build, compute, critique, repeat: Data analysis with latent variable models.  Annual Review of Statistics and Its Application, 1(1):203–232, 2014\.
  * [5] Leo Breiman.  Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author).  Statistical Science, 16(3):199 – 231, 2001.
  * [6] Peter Bühlmann.  Discussion of big bayes stories and BayesBag.  Statistical Science, 29(1), feb 2014.
  * [7] Jonah Gabry, Daniel Simpson, Aki Vehtari, Michael Betancourt, and Andrew Gelman.  Visualization in bayesian workflow.  Journal of the Royal Statistical Society Series A: Statistics in Societ, 182(2):389–402, jan 2019.
  * [8] Andrew Gelman.  Comment: Bayesian Checking of the Second Levels of Hierarchical Models.  Statistical Science, 22(3):349 – 352, 2007.
  * [9] Andrew Gelman, Xiao-Li Meng, and Hal Stern.  Posterior predictive assessment of model fitness via realized discrepancies.  Statistica Sinica, 6(4):733–760, 1996.
  * [10] Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C. Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák.  Bayesian workflow, 2020.
  * [11] Merlin Heidemanns, Andrew Gelman, and G. Elliott Morris.  An Updated Dynamic Bayesian Forecasting Model for the US Presidential Election.  Harvard Data Science Review, 2(4), oct 27 2020.  https://hdsr.mitpress.mit.edu/pub/nw1dzd02.
  * [12] Hsiang Hsu and Flavio Calmon.  Rashomon capacity: A metric for predictive multiplicity in classification.  In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh, editors, Advances in Neural Information Processing Systems, volume 35, pages 28988–29000. Curran Associates, Inc., 2022.
  * [13] Jonathan H. Huggins and Jeffrey W. Miller.  Robust inference and model criticism using bagged posteriors, 2020.
  * [14] Jonathan H. Huggins and Jeffrey W. Miller.  Reproducible Model Selection Using Bagged Posteriors.  Bayesian Analysis, 18(1):79 – 104, 2023.
  * [15] E. T. Jaynes.  Information theory and statistical mechanics.  Phys. Rev., 106:620–630, May 1957.
  * [16] Subhash Lele, Khurram Nadeem, and Byron Schmuland.  Estimability and likelihood inference for generalized linear mixed models using data cloning.  Journal of the American Statistical Association, 105:1617–1625, 12 2010.
  * [17] Arnaud Marsiglietti and Victoria Kostina.  A lower bound on the differential entropy of log-concave random vectors with applications.  Entropy, 20, 04 2017.
  * [18] Xiao-Li Meng.  Posterior predictive p-values.  The Annals of Statistics, 22(3):1142–1160, 1994.
  * [19] T. J. Mitchell and J. J. Beauchamp.  Bayesian variable selection in linear regression.  Journal of the American Statistical Association, 83(404):1023–1032, 1988.
  * [20] Juho Piironen and Aki Vehtari.  Sparsity information and regularization in the horseshoe and other shrinkage priors.  Electronic Journal of Statistics, 11(2):5018 – 5051, 2017.
  * [21] Nicholas G. Polson and James G. Scott.  501Shrink Globally, Act Locally: Sparse Bayesian Regularization and Prediction.  In Bayesian Statistics 9. Oxford University Press, 10 2011.
  * [22] Christopher Prosser and Jonathan Mellon.  The twilight of the polls? a review of trends in polling accuracy and the causes of polling misses.  Government and Opposition, 53(4):757–790, 2018.
  * [23] James M. Robins, Aad van der Vaart, and Valerie Ventura.  Asymptotic distribution of p values in composite null models.  Journal of the American Statistical Association, 95(452):1143–1156, 2000.
  * [24] Daniel Russo and James Zou.  How much does your data exploration overfit? controlling bias via information usage.  IEEE Trans. Inf. Theor., 66(1):302–323, jan 2020.
  * [25] Joel A. Tropp.  User-friendly tail bounds for sums of random matrices.  Foundations of Computational Mathematics, 12(4):389–434, aug 2011\.
  * [26] Rens van de Schoot, Sarah Depaoli, Andrew Gelman, Ruth King, Bianca Kramer, Kaspar Märtens, Mahlet G. Tadesse, Marina Vannucci, Joukje Willemsen, and Christopher Yau.  Bayesian statistics and modelling.  Nature Reviews Methods Primers, 1, January 2021.
  * [27] Aki Vehtari, Andrew Gelman, and Jonah Gabry.  Practical bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing, 27(5):1413–1432, aug 2016.
  * [28] Sumio Watanabe and Manfred Opper.  Asymptotic equivalence of bayes cross validation and widely applicable information criterion in singular learning theory.  Journal of machine learning research, 11(12), 2010.
  * [29] Yang Xie and Bradley Carlin.  Measures of bayesian learning and identifiability in hierarchical models.  Journal of Statistical Planning and Inference, 136:3458–3477, 10 2006.
  * [30] Aolin Xu and Maxim Raginsky.  Information-theoretic analysis of generalization capability of learning algorithms.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 2521–2530, Red Hook, NY, USA, 2017. Curran Associates Inc.
  * [31] Yoshitaka Yano, Stuart Beal, and Lewis Sheiner.  Evaluating pharmacokinetic/pharmacodynamic models using the predictive check.  Journal of pharmacokinetics and pharmacodynamics, 28:171–92, 04 2001.

##  Appendix A Basic Quantities and Relations from Information Theory

In this section, we provide statements of the basic results from information
theory that we make use of throughout this paper. Proofs of these results can
be found in any introductory course on information theory. We state all
results in terms of conditional entropies and mutual informations when
appropriate since these contain the non-conditional statements as special
cases. First we review relevant definitions.

######  Definition 7 (Basic Quantities of Information Theory).

Let $q(\boldsymbol{\theta},\mathbf{y})$ be some joint model. Then the entropy
of $\boldsymbol{\theta}$ is defined as

$h_{q(\boldsymbol{\theta})}(\boldsymbol{\theta})=-\mathbb{E}_{q(\boldsymbol{\theta})}\log
q(\boldsymbol{\theta}).$ (36)

The conditional entropy of $\boldsymbol{\theta}$ given $\mathbf{y}$ is just
the average entropy of the conditional distributions:

$h_{q(\boldsymbol{\theta},\mathbf{y})}\left(\boldsymbol{\theta}\mid\mathbf{y}\right)=\mathbb{E}_{q(\mathbf{y})}h_{q(\boldsymbol{\theta}\mid\mathbf{y})}\left(\boldsymbol{\theta}\right)=-\mathbb{E}_{q(\boldsymbol{\theta},\mathbf{y})}\log
q(\boldsymbol{\theta}\mid\mathbf{y}).$ (37)

The mutual information between $\boldsymbol{\theta}$ and $\mathbf{y}$ is the
amount by which entropy is expected to decrease after conditioning
$\mathbf{y}$:

$\mathbf{I}_{q}\left(\boldsymbol{\theta},\mathbf{y}\right)=h_{q(\boldsymbol{\theta})}\left(\boldsymbol{\theta}\right)-h_{q(\boldsymbol{\theta},\mathbf{y})}\left(\boldsymbol{\theta}\mid\mathbf{y}\right).$
(38)

Finally, if we extend our joint model to
$q(\boldsymbol{\theta},\mathbf{y},\mathbf{x})$ where $\mathbf{x}$ is any
additional quantity, then the conditional mutual information given
$\mathbf{x}$ is just the difference of the corresponding conditional
entropies:

$\mathbf{I}_{q}\left(\boldsymbol{\theta},\mathbf{y}\mid\mathbf{x}\right)=h_{q(\boldsymbol{\theta},\mathbf{x})}\left(\boldsymbol{\theta}\mid\mathbf{x}\right)-h_{q(\boldsymbol{\theta},\mathbf{y},\mathbf{x})}\left(\boldsymbol{\theta}\mid\mathbf{y},\mathbf{x}\right).$
(39)

The first important result allows us to break up an entropy or mutual
information expression additive over the components of vector arguments.

######  Lemma 1 (Chain Rule for Entropy and Mutual Information).

Let $q(\mathbf{x},\mathbf{y},\mathbf{z})$ be a joint model and suppose that
$\mathbf{x}$ can be partitioned into sub-vectors
$\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{m}\right)$ for some $m\geq 1$. Then
we have that

$h_{q}\left(\mathbf{x}\mid\mathbf{y}\right)=\sum_{i=1}^{m}h\left(\mathbf{x}_{i}\mid\mathbf{x}_{<i},\mathbf{y}\right),$

where $\mathbf{x}_{<i}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{i-1}\right)$
for $i\geq 2$, and $\mathbf{x}_{<1}=\\{\\}$. Furthermore, we have that

$\mathbf{I}_{q}\left(\mathbf{y},\mathbf{x}\mid\mathbf{z}\right)=\sum_{i=1}^{m}\mathbf{I}\left(\mathbf{y},\mathbf{x}_{i}\mid\mathbf{z},\mathbf{x}_{<i}\right).$

Next, it can be useful to express the (conditional) mutual information in
terms of the KL divergence. This can be done in two different ways.

######  Lemma 2 (Mutual Information as KL Divergence).

Let $q(\mathbf{x},\mathbf{y},\mathbf{z})$ be a joint model. Then we have

$\mathbf{I}(\mathbf{y},\mathbf{x}\mid\mathbf{z})=\mathbb{E}_{q(\mathbf{z})}D\left(q(\mathbf{x},\mathbf{y}\mid\mathbf{z})||q(\mathbf{x}\mid\mathbf{z})q(\mathbf{y}\mid\mathbf{z})\right)=\mathbb{E}_{q(\mathbf{y},\mathbf{z})}D\left(q(\mathbf{x}\mid\mathbf{y},\mathbf{z})||q(\mathbf{x}\mid\mathbf{z})\right).$

It is of fundamental importance that the KL divergence is always nonnegative,
which follows by an application of Jensen’s inequality.

######  Lemma 3 (Nonnegativity of the KL Divergence).

For any densities $p(\mathbf{y})$ and $q(\mathbf{y})$, we have

$D\left(p(\mathbf{y})||q(\mathbf{y})\right)\geq 0$

with equality if and only if $p(\mathbf{y})=q(\mathbf{y})$ $p$-almost surely.

This immediately implies nonnegativity of the mutual information, and in turn
the fact that

$h_{q(\mathbf{x},\mathbf{y})}\left(\mathbf{y}\mid\mathbf{x}\right)\leq
h_{q(\mathbf{y})}\left(\mathbf{y}\right)$

for any joint distribution $q(\mathbf{x},\mathbf{y})$.

It is often useful to know how these quantities operate under certain
transformations of the random quantities in terms of which they are defined.
This is characterized by the following result.

######  Lemma 4 (Entropy and Mutual Information Under Transformation).

Let $\mathbf{A}\in\mathbb{R}^{d\times d}$ be any invertible matrix and let
$\mathbf{y}^{\prime}=\mathbf{A}\mathbf{y}$. Then we have

$h_{q(\mathbf{y}^{\prime})}\left(\mathbf{y}^{\prime}\right)=h_{q(\mathbf{y})}\left(\mathbf{y}\right)+\log\left\lvert\det\mathbf{A}\right\rvert.$

Furthermore, if $\mathbf{y}^{\prime}=\mathbf{y}+c$ for any $c\in\mathbb{R}$,
then $h\left(\mathbf{y}\right)=h\left(\mathbf{y}^{\prime}\right)$. Thus, the
entropy is invariant under translations and orthogonal transformations. The
mutual information satisfies the stronger property of invariance under
arbitrary smooth reparametrizations of the individual arguments. Specifically,
let $\phi,\psi$ be smooth, invertible maps, and define
$\mathbf{y}^{\prime}=\phi(\mathbf{y})$ and
$\mathbf{x}^{\prime}=\psi(\mathbf{x})$. Then we have that

$\mathbf{I}(\mathbf{y}^{\prime},\mathbf{x}^{\prime}\mid\mathbf{z})=\mathbf{I}(\mathbf{y},\mathbf{x}\mid\mathbf{z}).$

The general behavior of the mutual information under potentially noninvertible
transformations is characterized by the data processing inequality.

######  Lemma 5 (Data Processing Inquality).

Let $q(\mathbf{x},\mathbf{y},\mathbf{z})$ be any distribution, and suppose
that $\mathbf{x}$ and $\mathbf{z}$ are conditionally independent given
$\mathbf{y}$. Then we have that

$\mathbf{I}(\mathbf{x},\mathbf{y})\geq\mathbf{I}(\mathbf{x},\mathbf{z}).$

In particular, the above inequality holds if $\mathbf{z}=\psi(\mathbf{y})$ for
any function $\psi$.

Finally, certain distributions maximize the entropy under certain conditions.
For our purposes, it suffices to note that normal distributions on
$\mathbb{R}^{d}$ maximize the entropy among all distributions with fixed
covariance matrix and support equal to $\mathbb{R}^{d}$.

######  Lemma 6 (Maximum Entropy of Normal).

Let $q(\mathbf{y})$ be any probability distribution supported on
$\mathbb{R}^{d}$, and let $p(\mathbf{y})$ be a normal distribution with any
mean and covariance matrix equal to the covariance $\boldsymbol{\Sigma}_{q}$
of $q(\mathbf{y})$. Then we have

$h_{q(\mathbf{y})}\left(\mathbf{y}\right)\leq
h_{p(\mathbf{y})}\left(\mathbf{y}\right)=\frac{1}{2}\log\left(\det\left(2\pi
e\boldsymbol{\Sigma}_{q}\right)\right).$

##  Appendix B Weak Identification Simulation Details

Algorithm 1 Posterior Bootstrap for Sampling Existing Subpopulations

1: Observed data $\mathbf{y}$, # of replications $R$, # of new samples per
subpopulation $M^{\mathrm{new}}$, # of posterior samples $S$.

2:for $r\leftarrow 1,\ldots,R$ do

3: Sample $(\theta_{1},\ldots\theta_{L},\sigma)\sim
p(\boldsymbol{\theta},\sigma\mid\mathbf{y})$.

4: for $l\leftarrow 1,\ldots L$ do

5: Sample
$\mathbf{y}^{\mathrm{new}}_{1l},\ldots,\mathbf{y}^{\mathrm{new}}_{M^{\mathrm{new}}l}\stackrel{{\scriptstyle
iid}}{{\sim}}p(y\mid\theta_{l},\sigma)$.

6: end for

7: $\mathbf{y}^{\mathrm{rep}}\leftarrow\left[\begin{smallmatrix}\mathbf{y}\\\
\mathbf{y}^{\mathrm{new}}\end{smallmatrix}\right]$.

8: Sample $\mu_{(1)},\ldots,\mu_{(S)}\stackrel{{\scriptstyle
iid}}{{\sim}}p\left(\mu\mid\mathbf{y}^{\mathrm{rep}}\right)$

9:
$\hat{\sigma}_{r}\leftarrow\mathsf{sd}\left(\mu_{(1)},\ldots,\mu_{(S)}\right)$

10:end for

11:return $\\{\hat{\sigma}_{1},\ldots,\hat{\sigma}_{R}\\}$.

Algorithm 2 Posterior Bootstrap for Sampling New Subpopulations

1:Observed data $\mathbf{y}$, # of replications $R$, # of new subpopulations
$L^{\mathrm{new}}$, # of posterior samples $S$.

2:for $r\leftarrow 1,\ldots,R$ do

3: Sample $(\mu,\tau,\sigma)\sim p(\mu,\tau,\sigma\mid\mathbf{y})$.

4: Sample $(\theta_{1},\ldots\theta_{L^{\mathrm{new}}})\sim
p(\boldsymbol{\theta}\mid\mu,\tau)$.

5: for $l\leftarrow 1,\ldots L^{\mathrm{new}}$ do

6: Sample
$\mathbf{y}^{\mathrm{new}}_{1l},\ldots,\mathbf{y}^{\mathrm{new}}_{Ml}\stackrel{{\scriptstyle
iid}}{{\sim}}p(y\mid\theta_{l},\sigma)$.

7: end for

8:
$\mathbf{y}^{\mathrm{rep}}\leftarrow\left[\begin{smallmatrix}\mathbf{y}&\mathbf{y}^{\mathrm{new}}\end{smallmatrix}\right]$.

9: Sample $\mu_{(1)},\ldots,\mu_{(S)}\stackrel{{\scriptstyle
iid}}{{\sim}}p\left(\mu\mid\mathbf{y}^{\mathrm{rep}}\right)$

10:
$\hat{\sigma}_{r}\leftarrow\mathsf{sd}\left(\mu_{(1)},\ldots,\mu_{(S)}\right)$

11:end for

12:return $\\{\hat{\sigma}_{1},\ldots,\hat{\sigma}_{R}\\}$.

We simulated our three data sets in Section 5 from the following model.

$\begin{gathered}\mathbf{y}_{ml}\mid\boldsymbol{\theta}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\theta_{l},\sigma_{*}\right),\qquad\boldsymbol{\theta}_{l}\stackrel{{\scriptstyle
iid}}{{\sim}}\mathrm{normal}\left(\mu_{*},\tau_{*}\right),\\\ \end{gathered}$
(40)

The $\sigma_{*}$, $\mu_{*}$, and $\tau_{*}$ are all hyperparameters. We fixed
$\mu_{*}$ and $\tau_{*}$ for all data sets, but we varied $\sigma_{*}$ from
$\tau_{*}/2$ to $2\tau_{*}$. For our purposes, the key differentiator of the
data sets is the ratio of the variance of the subpopulation-level (sample)
means to the overall (sample) variance. When this statistic is closer to $1$,
the between-subpopulation variation swamps the within-population variation,
and the reverse is true when this statistic is close to $0$. For the three
resulting data sets, this statistic was approximately $0.45$, $0.6$, and
$0.95$ respectively. The complete data sets along with the R and Stan code
used to generate and analyze them are available on the companion github
repository https://github.com/collin-cademartori/BayesianModelExpansionPaper/.

Algorithms 1 and 2 give detailed pseudo-code for the two resampling procedures
we used to simulate the process of sampling data from (a) the same
subpopulations from which the original data was sampled and (b) new
subpopulations within the larger superpopulation. Each of these algorithms was
run with $R=500$ replication of the resampling scheme, $S=2000$ samples drawn
from each of the resulting posterior distributions, and with
$M^{\mathrm{new}}=8$ and $L^{\mathrm{new}}=20$ respectively.

##  Appendix C Details of the Election Forecasting Model

The bias terms $\beta_{i}$ in (29) are decomposed into several further terms
which are designed to capture the following sources of polling bias:

  * •

Pollster-level “house” effects, i.e. the observed phenomenon that most
pollsters exhibit a nonzero, temporally stable bias toward one party or the
other.

  * •

Poll mode effects reflecting the observed phenomenon that the averages of
phone- and internet-based polls tend to differ.

  * •

Poll population effects reflecting the observed phenomenon that the averages
of “registered voter” and “likely voter” polls tend to differ.

  * •

Partisanship nonresponse effects that account for the fact that differential
rates of response between members of the two parties can alter poll results.
This term is only included for polls that do not attempt to adjust their
results for the partisan composition of their sample. Unlike the previous
terms, this term varies through time (and is assigned a time-series prior) in
order to account for the fact that differential partisan nonresponse has been
observed to vary over the course of an election cycle.

  * •

State-level measurement error effects reflecting the observation that polls of
some states have historically been more accurate than others.

  * •

Poll-level random measurement error, i.e. a catch-all term representing any
additional variance in the average support for the Democrat among the polls’
sampling frames unaccounted for by the above terms.

A complete description of the polling model specification can be found in
[11], and the code and data for the original model can be accessed at
https://github.com/TheEconomist/us-potus-model. The code for our expanded
models and our conditional posterior checks can be found at
https://github.com/collin-cademartori/BayesianModelExpansionPaper.

##  Appendix D Conditional Mutual Information Bounds

In this section we state and prove Theorem 4, our lower bound on the
conditional mutual information. To establish the lower bound, we will need to
use the following key fact on the monotonicity of the conditional mutual
information under conditioning on successive replications.

######  Lemma 7 (Decreasing Conditional Mutual Information).

Let $p(\boldsymbol{\theta},\mathbf{y})$ be a joint model of parameters and
data. For any $M,R\geq 1$, let

$(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M+R)})\stackrel{{\scriptstyle
iid}}{{\sim}}p(\mathbf{y}\mid\boldsymbol{\theta}).$

Then we have that
$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}^{(M)}\mid\mathbf{y},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M-1)})$
is decreasing in $M$. Furthermore, we have that

$\mathbf{I}(\boldsymbol{\theta},\mathbf{y}^{(1)}\mid\mathbf{y})\geq\frac{1}{M}\mathbf{I}(\boldsymbol{\theta},(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)})\mid\mathbf{y}^{(M+1)},\ldots\mathbf{y}^{(M+R)})$

for all $M\geq 1$.

###### Proof.

Letting
$\overline{\mathbf{y}}_{M}=\left(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
and
$\widetilde{\mathbf{y}}_{R}=\left(\mathbf{y}^{(M+1)},\ldots,\mathbf{y}^{(M+R)}\right)$,
we have that

$\displaystyle\mathbf{I}\left(\boldsymbol{\theta},\overline{\mathbf{y}}_{M}\mid\widetilde{\mathbf{y}}_{R}\right)$
$\displaystyle=h\left(\overline{\mathbf{y}}_{M}\mid\widetilde{\mathbf{y}}_{R}\right)-h\left(\overline{\mathbf{y}}_{M}\mid\widetilde{\mathbf{y}}_{R},\boldsymbol{\theta}\right)$
$\displaystyle=h\left(\overline{\mathbf{y}}_{M}\mid\widetilde{\mathbf{y}}_{R}\right)-h\left(\overline{\mathbf{y}}_{M}\mid\boldsymbol{\theta}\right)$
(41) $\displaystyle\leq
h\left(\overline{\mathbf{y}}_{M}\mid\mathbf{y}^{(M+1)}\right)-h\left(\overline{\mathbf{y}}_{M}\mid\boldsymbol{\theta}\right)$
(42)
$\displaystyle=h\left(\overline{\mathbf{y}}_{M}\mid\mathbf{y}^{(M+1)}\right)-h\left(\overline{\mathbf{y}}_{M}\mid\boldsymbol{\theta},\mathbf{y}^{(M+1)}\right)$
(43)
$\displaystyle=\mathbf{I}\left(\boldsymbol{\theta},\overline{\mathbf{y}}_{M}\mid\mathbf{y}^{(M+1)}\right),$

where (41) follows from the conditional independence of the $\mathbf{y}^{(i)}$
given $\boldsymbol{\theta}$, (42) follows from the fact that conditioning
decreases entropy, and (43) again follows in the same way as (41).

Now, by the chain rule for mutual information and the same argument as above,
we can bound this latter conditional mutual information as

$\displaystyle\mathbf{I}\left(\boldsymbol{\theta},\overline{\mathbf{y}}_{M}\mid\mathbf{y}^{(M+1)}\right)$
$\displaystyle=\sum_{j=1}^{M}\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}^{(j)}\mid\mathbf{y}^{(M+1)},\overline{\mathbf{y}}_{j-1}\right)$
$\displaystyle\leq\sum_{j=1}^{M}\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}^{(j)}\mid\mathbf{y}^{(M+1)}\right)$
$\displaystyle=M\mathbf{I}\left(\boldsymbol{\theta},\mathbf{y}^{(1)}\mid\mathbf{y}\right),$

where $\overline{\mathbf{y}}_{0}=\left\\{\right\\}$, and where the last
inequality follows from the fact that the distributions of
$\left(\boldsymbol{\theta},\mathbf{y}^{(j)},\mathbf{y}^{(M+1)}\right)$ are all
equal to the distribution of
$\left(\boldsymbol{\theta},\mathbf{y}^{(1)},\mathbf{y}\right)$. Now dividing
both sides by $M$ completes the proof. ∎

Before stating our main lemma, we first give a definition that extends the
notion of a subexponential distribution to random symmetric matrices.

######  Definition 8 (Subexponential Random Matrix).

A random symmetric matrix $\mathbf{M}$ is said to be
$(\alpha,\beta)$-subexponential if, for all $k\geq 2$, we have

$\mathbb{E}\mathbf{M}^{k}\prec\alpha\frac{k!}{2}\beta^{k-2}\mathbf{I},$

where the inequality denotes the Loewner order.

We will also say that a random matrix is $\gamma$-subexponential if it is
$(\alpha,\beta)$-subexponential for some $\alpha,\beta$ such that
$\gamma=\alpha+\beta$.

######  Lemma 8 (Fisher Information Lower Bound).

For $M,R\geq 1$, define the $(M+R)$-replicated model:

$p\left(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M+R)},\boldsymbol{\theta}\right)=p(\boldsymbol{\theta})\prod_{i=1}^{M+R}p\left(\mathbf{y}^{(i)}\mid\boldsymbol{\theta}\right).$

Suppose for $R$ sufficiently large, we have that

  * •

the posterior distributions
$p\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(R)}\right)$
are normal,

  * •

the observed information matrix of $p(\mathbf{y},\boldsymbol{\theta})$ is
$(\alpha,\beta)$-subexponential for some $\alpha,\beta>0$ (i.e. the observed
information does not have heavy tails),

  * •

$\mathbb{E}\lambda_{d}(\boldsymbol{\Sigma})$,
$\mathbb{E}\lambda^{-1}_{1}(\boldsymbol{\Sigma})$,
$\mathbb{E}\lambda^{2}_{d}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)$,
and
$\mathbb{E}\lambda_{1}^{-1}(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})))$
are bounded by some $B>0$ where
$\boldsymbol{\Sigma}=\mathrm{Cov}\left(\boldsymbol{\theta}\mid\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$
(i.e. the posterior covariance and Fisher information are neither too small
nor too large on average).

Then for $C$ a constant depending on $\gamma=\alpha+\beta$ and $B$, we have
for $M$ sufficiently large that

$\frac{1}{M}\mathbf{I}(\boldsymbol{\theta},\overline{\mathbf{y}}_{M}\mid\widetilde{\mathbf{y}}_{R})\geq\frac{C}{\log
d}\mathrm{tr}\left(\mathbb{E}_{p(\boldsymbol{\theta},\mathbf{y})}\boldsymbol{\Sigma}^{1/2}_{\overline{\mathbf{y}}_{R}}\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\boldsymbol{\Sigma}_{\overline{\mathbf{y}}_{R}}^{1/2}\right).$
(44)

We give the proof in the case $R=1$ and write $\mathbf{y}$ for
$\mathbf{y}^{(M+1)}$ for simplicity, but the proof for $R>1$ follows in the
exact same way with
$\left(\mathbf{y}^{(M+1)},\ldots,\mathbf{y}^{(M+R)}\right)$ replacing
$\mathbf{y}$ throughout.

###### Proof.

By the invariance of the mutual information under invertible transformations,
if we define
$\widetilde{\boldsymbol{\theta}}=\boldsymbol{\Sigma}^{-1/2}_{\mathbf{y}}\boldsymbol{\theta}$,
we get

$\mathbf{I}(\boldsymbol{\theta},(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)})\mid\mathbf{y})=\mathbf{I}(\widetilde{\boldsymbol{\theta}},(\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)})\mid\mathbf{y})\geq\mathbb{E}\log\left[\det\widetilde{\boldsymbol{\Sigma}}_{M}^{-1}\right],$
(45)

where the last inequality follows by expressing the conditional mutual
information as the entropy difference
$h\left(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}\right)-h\left(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\right)$,
plugging in the expressions for the entropy of multivariate normal
distributions, and simplifying.

Now for $1\leq i\leq d$, let $\lambda^{d}_{i}(\cdot):H_{d}\to\mathbb{R}$ be
the map from $d\times d$ Hermitian matrices to their $i^{\mathrm{th}}$
eigenvalue under increasing order. We will usually write $\lambda_{i}$ for
$\lambda_{i}^{d}$ when the matrix dimension is clear from context. We will
henceforth write $\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}$ for
$\widetilde{\boldsymbol{\Sigma}}_{\mathbf{y},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}}^{-1}$
for convenience.

Furthermore, define
$\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})=\frac{1}{M}\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})$
where $\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})=-\log
p(\mathbf{y}^{(i)}\mid\widetilde{\boldsymbol{\theta}})$, and

$\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}=\left\\{\left\|\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})-\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right\|_{\mathrm{op}}<\frac{\delta}{2}\right\\}.$
(46)

Finally, define the event

$\mathcal{G}=\left\\{\lambda_{1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)\geq\delta\right\\}.$

The proof will now proceed by analyzing the decomposition

$\displaystyle\mathbb{E}\log\left[\det\widetilde{\boldsymbol{\Sigma}}_{M}^{-1}\right]$
$\displaystyle=\sum_{i=1}^{d}\mathbb{E}\left[\log\left(\lambda_{i}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]+\sum_{i=1}^{d}\mathbb{E}\left[\log\left(\lambda_{i}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{E}^{c}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]$
$\displaystyle\hskip
28.45274pt+\sum_{i=1}^{d}\mathbb{E}\left[\log\left(\lambda_{i}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{G}^{c}}\right],$
(47)

where, on the right hand side, we now take the expectations over
$\left(\mathbf{y}^{(0)},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)},\widetilde{\boldsymbol{\theta}}\right)$.
Label these terms T1 - T3.

Lower bound for T3.

First we note that

$\displaystyle\mathbb{E}\left[\log\left(\lambda_{i}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{G}^{c}}\right]$
$\displaystyle\geq\mathbb{E}\left[\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{G}^{c}}\right]$
$\displaystyle=-\mathbb{E}\left[\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)\mathbbm{1}_{\\{\lambda_{1}^{-1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)>\delta^{-1}\\}}\right]$
$\displaystyle\geq-\sqrt{\mathbb{E}\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)^{2}}\sqrt{\mathbb{P}\left(\lambda_{1}^{-1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)\geq\delta^{-1}\right)}$
$\displaystyle\geq-\delta^{1/2}\sqrt{\mathbb{E}\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)^{2}\mathbb{E}\lambda^{-1}_{1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)}$
$\displaystyle\geq-\delta^{1/2}\sqrt{\left(\mathbb{E}\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)+\mathbb{E}\lambda^{-1}_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)\mathbb{E}\lambda^{-1}_{1}\left(\boldsymbol{\mathcal{I}}(\boldsymbol{\theta})\right)}$
$\displaystyle\geq-\sqrt{2}\delta^{1/2}B.$

So T3 is lower bounded by $-\sqrt{2}d\delta^{1/2}B$.

Lower Bound for T2.

We now lower bound T2. First observe that

$\displaystyle\mathbb{E}\left[\log\left(\lambda_{i}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{E}^{c}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]$
$\displaystyle\geq\mathbb{E}\left[\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}^{-1}_{M}\right)\right)\mathbbm{1}_{\mathcal{E}^{c}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]$
$\displaystyle=-\mathbb{E}\left[\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)\mathbbm{1}_{\mathcal{E}^{c}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]$
$\displaystyle\geq-\sqrt{\mathbb{E}\left[\log\left(\lambda_{1}\left(\widetilde{\boldsymbol{\Sigma}}_{M}\right)\right)^{2}\right]}\sqrt{\mathbb{P}\left(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\cap\mathcal{G}\right)}$

Now we observe that we can upper bound the probability factor as

$\mathbb{P}\left(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\cap\mathcal{G}\right)=\mathbb{E}_{p(\boldsymbol{\theta},\mathbf{y})}\mathbb{P}\left(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\cap\mathcal{G}\mid\widetilde{\boldsymbol{\theta}},\mathbf{y}\right)\leq\mathbb{E}_{p(\boldsymbol{\theta},\mathbf{y})}\mathbb{P}\left(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\mid\widetilde{\boldsymbol{\theta}},\mathbf{y}\right).$

Since we assume that the observed information matrices
$\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})$
are sub-exponential in the sense of (8), we can apply the matrix Bernstein
inequality conditional on $(\widetilde{\boldsymbol{\theta}},\mathbf{y})$ with
rank of the Fisher information
$\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})$ to see that, if
$\delta\leq 1$,

$\mathbb{P}\left(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\mid\widetilde{\boldsymbol{\theta}}\right)\leq
2d\exp\left(-M\frac{\delta^{2}}{4\alpha+2\beta\delta}\right)\leq
2d\exp\left(-M\frac{\delta^{2}}{4\gamma}\right),$ (48)

where $\alpha,\beta$ are the sub-exponential parameters controlling the tails
of the observed information, $\gamma=\alpha+\beta$. Since the right hand side
is free of $(\widetilde{\boldsymbol{\theta}},\mathbf{y})$, we obtain the same
upper bound for the marginal probability
$\mathbb{P}(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c}\cap\mathcal{G})$.
This gives an overall bound for T2 of
$-\sqrt{2}Bd\exp\left(-M\frac{\delta^{2}}{8\gamma}\right)$.

Lower Bound for T1.

Now, for any Gaussian distribution, the Hessian of the potential function is
exactly the precision matrix. Combining this with the fact that derivatives of
the log normalizing constant vanish and the fact that
$\mathrm{Cov}\left(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}\right)=\mathbf{\mathsf{I}}$,
we have that, for all $\widetilde{\boldsymbol{\theta}}\in\mathbb{R}^{d}$,

$\displaystyle\widetilde{\boldsymbol{\Sigma}}^{-1}_{\mathbf{y}^{(0)},\ldots,\mathbf{y}^{(M)}}$
$\displaystyle=\mathbf{H}_{-\log
p(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(0)},\ldots,\mathbf{y}^{(M)})}(\widetilde{\boldsymbol{\theta}})$
$\displaystyle=\mathbf{H}_{-\log(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y})}(\widetilde{\boldsymbol{\theta}})+\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})$
$\displaystyle=\mathbf{\mathsf{I}}+\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}}).$
(49)

Now note that

$\lambda_{i}\left(\mathbf{\mathsf{I}}+\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})\right)=1+\lambda_{i}\left(\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})\right)$
(50)

Combining the observations of (50) and (D) with the identity (45), we obtain

$\displaystyle\mathbf{I}(\widetilde{\boldsymbol{\theta}},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\mid\mathbf{y})$
$\displaystyle=\frac{1}{2}\mathbb{E}\sum_{i=1}^{d}\log\left(1+\lambda_{i}\left(\sum_{i=1}^{M}\mathbf{H}_{\phi(\widetilde{\boldsymbol{\theta}}\mid\mathbf{y}^{(i)})}(\widetilde{\boldsymbol{\theta}})\right)\right)$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{d}\mathbb{E}\log\left(1+M\lambda_{i}\left(\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})\right)\right).$
(51)

Now Weyl’s inequalities imply that

$\left\lvert\lambda_{i}\left(\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})\right)-\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\right\rvert\leq\left\lvert\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})-\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right\rvert_{\mathrm{op}}.$
(52)

Using this, we have for all $1\leq i\leq M$ that

$\displaystyle\mathbb{E}_{p(\overline{\mathbf{y}}\mid\widetilde{\boldsymbol{\theta}},\mathbf{y})}$
$\displaystyle\left[\log\left(1+M\lambda_{i}\left(\mathbf{I}_{\mathrm{obs}}(\overline{\mathbf{y}},\widetilde{\boldsymbol{\theta}})\right)\right)\mathbbm{1}_{\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}\cap\mathcal{G}}\right]$
$\displaystyle\hskip
28.45274pt\geq\log\left(1+\frac{M}{2}\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\right)\mathbbm{1}_{\mathcal{G}}(\widetilde{\boldsymbol{\theta}})\left(1-\mathbb{P}_{p(\overline{\mathbf{y}}\mid\widetilde{\boldsymbol{\theta}},\mathbf{y})}(\mathcal{E}_{\widetilde{\boldsymbol{\theta}}}^{c})\right)$
(53)

where the expectations are over
$\overline{\mathbf{y}}=\left(\mathbf{y}^{(1)},\ldots\mathbf{y}^{(M)}\right)$
conditional on $(\widetilde{\boldsymbol{\theta}},\mathbf{y})$.

Now combining (D), and (51), and letting

$M=8k\frac{\gamma}{\delta^{2}}\log\left(2d\right),$

we obtain that

$\displaystyle\mathbf{I}$
$\displaystyle(\widetilde{\boldsymbol{\theta}},\mathbf{y}^{(1)},\ldots,\mathbf{y}^{(M)}\mid\mathbf{y})$
$\displaystyle\hskip
28.45274pt\geq\frac{1}{4}\sum_{i=1}^{d}\mathbb{E}_{p(\widetilde{\boldsymbol{\theta}},\mathbf{y})}\left[\log\left(1+\frac{4k\gamma}{\delta^{2}}\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\right)\mathbbm{1}_{\mathcal{G}}(\widetilde{\boldsymbol{\theta}})\right]-cd\left[\delta^{1/2}+2^{-k}\right].$
(54)

If we now define the function

$\sigma^{\tau}(x)=\min\left(x,\tau\right)$ (55)

to be the $\tau$-truncation of $x$, then we have for $\delta^{-1}$ and $\tau$
sufficiently large that the first term above is lower bounded by

$\displaystyle\mathbb{E}_{p(\widetilde{\boldsymbol{\theta}},\mathbf{y})}$
$\displaystyle\sum_{i=1}^{d}\log\left[1+\frac{4k\gamma}{\delta^{2}}\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\right]\mathbbm{1}_{\mathcal{G}}(\widetilde{\boldsymbol{\theta}})$
$\displaystyle\geq\mathbb{E}_{p(\widetilde{\boldsymbol{\theta}},\mathbf{y})}\sum_{i=1}^{d}\log\left[1+\frac{4k\gamma}{\delta^{2}}\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\right]\mathbbm{1}_{\\{\delta<\lambda_{1}(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}}))\leq\lambda_{d}(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}}))<\tau\\}}$
$\displaystyle\geq\frac{\log\left(1+4k\gamma\tau/\delta^{2}\right)}{\tau}\sum_{i=1}^{d}\mathbb{E}_{p(\widetilde{\boldsymbol{\theta}},\mathbf{y})}\lambda_{i}\left(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}})\right)\mathbbm{1}_{\\{\delta<\lambda_{1}(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}}))\leq\lambda_{d}(\boldsymbol{\mathcal{I}}(\widetilde{\boldsymbol{\theta}}))<\tau\\}}$