[3]

# Population Predictive Checks

Rajesh Ranganath  Department of Computer Science and Center for Data Science;
New York University  David M. Blei  Department of Computer Science and
Department of Statistics; Columbia University

###### Abstract

Bayesian modeling has become a staple for researchers analyzing data. Thanks
to recent developments in approximate posterior inference, modern researchers
can easily build, use, and revise complicated Bayesian models for large and
rich data. These new abilities, however, bring into focus the problem of model
assessment. Researchers need tools to diagnose the fitness of their models, to
understand where a model falls short, and to guide its revision. In this paper
we develop a new method for Bayesian model checking, the population predictive
check (pop-pc). pop-pcs are built on posterior predictive checks (ppcs), a
seminal method that checks a model by assessing the posterior predictive
distribution on the observed data. Though powerful, ppcs use the data
twice—both to calculate the posterior predictive and to evaluate it—which can
lead to overconfident assessments. pop-pcs, in contrast, compare the posterior
predictive distribution to the population distribution of the data. This
strategy blends Bayesian modeling with frequentist assessment, leading to a
robust check that validates the model on its generalization. Of course the
population distribution is not usually available; thus we use tools like the
bootstrap and cross validation to estimate the pop-pc. Further, we extend pop-
pcs to hierarchical models. We study pop-pcs on classical regression and a
hierarchical model of text. We show that pop-pcs are robust to overfitting and
can be easily deployed on a broad family of models.

##  1 Introduction

Bayesian modeling has become an important tool for the modern applied
statistician. In Bayesian modeling, we cast our knowledge and assumptions into
a probabilistic model of hidden and observed variables. Then we condition on a
dataset and perform inference of the hidden variables by calculating the
posterior distribution. The posterior is key both for making predictions about
the future and for using the model to understand more qualitative properties
of the data.

Calculating the posterior has traditionally been the bottleneck to practical
Bayesian analysis. However, Bayesian methods have been made more practical
thanks to new methods for approximate posterior inference (Beaumont et al.,,
2002; Hoffman et al.,, 2013; Hoffman and Gelman,, 2014; Ranganath et al.,,
2014; Angelino et al.,, 2016; Blei et al.,, 2017; Srivastava et al.,, 2018)
and corresponding software for easily devising and using Bayesian models (Stan
Development Team,, 2015; Mansinghka et al.,, 2014; Wood et al.,, 2014; Goodman
and Stuhlmuller,, 2015; Tran et al., 2016a, ). With these innovations,
Bayesian statisticians can draw from an extensive catalog of models to develop
tailored data analyses for specific applied problems.

This freedom—to be able to develop and compose many models—has turned the
practice of model building into a cycle (Blei,, 2014; Gelman et al.,, 1995).
We iteratively construct more complex models from simpler ones, using our
assessments of earlier “drafts” to guide the structure of subsequent
revisions. Currently, this process is largely driven by creativity, intuition,
and experience. But it brings into sharp focus the problem of model
diagnostics. How do we navigate the space of models we can use for a problem?
How do we decide when a model needs to change? In what ways should we change
it? Are there statistical tools that can guide this craft? This paper develops
a new way to diagnose Bayesian models.

One tool that aims to provide flexible model diagnostics is the posterior
predictive check (ppc), the result of a line of thinking that sought to adapt
classical goodness-of-fit tests to Bayesian statistics (Guttman,, 1967;
Rubin,, 1984). A ppc formalizes the following check: “If my model is good (in
a relevant way) then its posterior predictive distribution will generate data
that looks like my observations (in that same way).” The main idea is to
locate the observed data in a predictive distribution of data sets; observed
data that is unlikely under the posterior predictive distribution indicates
model misfit.

But there is a crucial issue with the ppc—it uses the data twice. First, the
data are used to construct the posterior predictive distribution of data sets.
Second, the data are used as the “reference,” i.e., the observations that
should be well-modeled by the posterior predictive distribution. Thus the ppc
checks the fit of the same data that the model is built to predict, which
makes it prone to overconfidence about the model (Gelfand et al.,, 1992;
Chatfield,, 1995; Bayarri and Berger,, 2000; Robins et al.,, 2000). A ppc can
lead the modeler to undesirable conclusions and, within the
model/diagnose/revise cycle, can send the modeler to overfit their data.

This paper introduces a solution to this problem. We develop the population
predictive check (pop-pc), a new method for Bayesian model checking. The pop-
pc combines the Bayesian ppc with the frequentist idea of the population
distribution. The result changes the basic intuition of a diagnostic: “If my
model is good (in a relevant way) then data drawn from the posterior
predictive distribution will look like data from the true population (in that
way).” This alternative intuition leads to model checks that do not overfit.

As a practical matter, the population distribution is unknown. To approximate
it, the pop-pc appeals to familiar frequentist concepts such as the bootstrap
(Efron,, 1979), plug-in principle (Efron and Tibshirani,, 1993), and cross
validation (Stone,, 1974; Geisser,, 1975). This perspective leads to a
practical method for diagnosing Bayesian models, one that naturally evokes
ideas around held-out data and out-of-sample generalization. (It also
generalizes and justifies a research literature about using data-splitting to
assess models; we review this literature below.)

Concretely, the recipe for computing a pop-pc on a fixed data set is as
follows. First, sample a new data set from the empirical distribution. Next,
use this data set to construct a posterior predictive distribution. Finally,
evaluate the posterior predictive probability of the observed data that were
omitted from the resampled data, i.e., held-out data. The model is misfit when
the the posterior predictive does not model the held-out data well. We can
repeat this procedure across many resampled data sets.

This paper is organized as follows. Section 2 discusses the historical
development of Bayesian model evaluation, and where this paper fits in.
Section 3.1 reviews ppc s and discusses the double-use of data; Section 3.2
develops pop-pc s; Section 3.3 provides a simple example to demonstrate both;
and Section 3.4 extends the ideas to realized discrepancies, those that
include latent variables, and to multi-level models. Section 4 discusses
computation, providing algorithms to approximate pop-pc s by blending
approximate Bayesian inference and frequentist methods for estimating the
population distribution. Section 5 demonstrates the algorithms empirically on
two types of models, a regression model and a hierarchical document model.
Section 6 concludes the paper.

##  2 Related work

This work builds on posterior predictive check (ppc) s, which were developed
in several seminal papers. First, inspired by earlier ideas (Geisser,, 1975),
Box, (1980) used the prior predictive distribution as the distribution of
data; this is a prior predictive check. Later, Rubin, (1984) mimicked Box’s
framework, but replaced the prior predictive with the posterior predictive;
this is both a more practical strategy for diagnosing models and one that is,
in Rubin’s language, “Bayesianly justifiable.” (Independently, Guttman, (1967)
proposed a similar approach.) Finally, Gelman et al., (1996) showed how to
develop diagnostic functions of the data—termed discrepancy functions—that
depend on both a data set and the hidden variables. This made the ppc even
more practical, allowing investigators to define focused diagnostic checks for
examining the performance of their models.

We first describe Box, (1980) and Rubin, (1984), and how this work fits into
their thinking. We then describe subsequent related work.

Box, (1980) posits two essential activities for the applied statistician:
inference of unknown quantities and criticism of a model. His position is that
inference should be Bayesian—we condition on our data and calculate posterior
quantities of unknowns—but that criticism is based on sampling theory, a
frequentist perspective.

Figure 1: A schematic diagram that relates Box, (1980), Rubin, (1984), and
this paper. This diagram posits the data come from an unknown population
distribution. Box, (1980) uses the data as reference in the marginal
distribution induced by the model (top); Rubin, (1984) uses the data as
reference in the posterior predictive distribution induced by the model
(middle); our method uses the population distribution as a reference in the
posterior predictive (bottom).

Set a model $p(\bm{\theta},\mathbf{y})$, where $\bm{\theta}$ is unknown and
$\mathbf{y}$ is the data. We observe data $\mathbf{y}^{\textrm{obs}}$ and
calculate the posterior $p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})$—this
tells us about the unknown quantities that help explain the observations. But
this inference is subject to the “correctness” of the model, where a correct
model provides a distribution of data similar to their true distribution. To
assess that correctness we criticize the model, asking if the observations
$\mathbf{y}^{\textrm{obs}}$ could have plausibly come from it. For Box, this
involves forming the marginal distribution $p(\mathbf{y})=\int
p(\bm{\theta},\mathbf{y})\mathop{}\\!\mathrm{d}\bm{\theta}$ and then locating
the $\mathbf{y}^{\textrm{obs}}$ within it; this strategy is depicted in Figure
1 (top). Box argues that if the $\mathbf{y}^{\textrm{obs}}$ could have come
from the marginal then we can use
$p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})$ as an inference of the unknown
variables that explain them.

Box’s perspective clarifies the roles for Bayesian inference and frequentist
theory, but it misses how modern statisticians use Bayesian models. Box seeks
joint distributions whose marginals capture the data, that is, models we can
believe in. In contrast, statisticians use Bayesian models to express a rich
space of possible explanations and then use data (through the posterior) to
focus on the relevant portion of that space. In real-world applications,
realistic marginal distributions are less important than realistic predictive
distributions, where posterior inference points the model towards a good
distribution of future data.

It is in this spirit that Rubin, (1984) adapts Box’s setup to a more practical
end. He sets up a model as above and does inference with the posterior. To
criticize the model, however, he calculates the posterior predictive
distribution $p(\mathbf{y}\,|\,\mathbf{y}^{\textrm{obs}})=\int
p(\mathbf{y}\,|\,\theta)p(\theta\,|\,\mathbf{y}^{\textrm{obs}})$ and locates
the observations within it; see Figure 1 (middle). The essence of the ppc is
that a good model will produce a posterior predictive that captures the
observed data, and this justifies interpreting the posterior and using it in
other downstream tasks.

Rubin’s motivation was philosophical as well as practical. He emphasized that
the ppc is “Bayesianly justifiable,” in that it always conditions on the
observed data. (This is in contrast to Box, who checks the model with
$p(\mathbf{y})$ rather than $p(\mathbf{y}\,|\,\mathbf{y}^{\textrm{obs}})$.)
However, the philosophical advantage of the ppc is also its flaw: it asks the
posterior predictive to capture the same data that it has already conditioned
on. While this does not violate any formal principles of Bayesian statistics,
it can lead to diagnostics that are overly optimistic about the quality of the
posterior predictive. The ppc has been criticized for this; see the excellent
discussion and references in Bayarri and Castellanos, (2007).

In this paper, we circle back to Box. Like Rubin, we build diagnostics that
stem from the posterior predictive. But, like Box, we maintain the idea that
model criticism is a question of sampling theory. Ideally, we ask if the
posterior predictive distribution captures the population distribution of the
data; see Figure 1 (bottom). When such an ideal check passes, we can interpret
and use Bayesian inference of the latent variables—their posterior
distribution provides a plausible mechanism for generating data from the
population. This is the population predictive check (pop-pc).

Other related work. We described how pop-pcs assess Bayesian models by
comparing the posterior predictive distribution to the population
distribution. With fixed data, we approximate them using frequentist concepts
like the bootstrap and cross validation. pop-pcs improve on ppcs (Rubin,,
1984), which can overfit on account of the double use of the data. Here we
discuss alternative methods for addressing this issue.

Gelfand et al., (1992) developed cross-validated checks. They iteratively hold
out each data point, condition on the remaining data, and compare samples from
the corresponding posterior predictive distribution to the held-out point.
This is a type of pop-pc, one that uses a leave-one-out cross validation
estimator (see Section 4). Similar strategies are discussed in Draper, (1996);
O’Hagan, (2003); Dahl et al., (2007); pop-pc’s generalize and justify these
approaches.

Bayarri and Berger, (2000) propose partial posterior predictive and
conditional predictive checks. For a given discrepancy, partial posterior
predictive checks use simulations from the posterior predictive distribution
formed with a likelihood that is conditioned on the discrepancy. Similarly,
conditional predictive checks form simulations from a distribution conditioned
on an alternative statistic, chosen to not have information about the
discrepancy. Both strategies avoid the double use of the data and enjoy
asymptotically uniform p-values (Robins et al.,, 2000). These approaches are
limited, however, in that information needed to build a meaningful posterior
may also be found in the discrepancy. Further, developing new conditional
posterior approximations for each new discrepancy can be computationally
prohibitive.

Robins et al., (2000) propose calibrating ppc p-values by computing the
empirical distribution using the bootstrap; Hjort et al., (2006) propose a
similar calibration. These approaches apply to $p$-values and tail
probabilities. But the full ppc framework and pop-pc framework allow for other
types of checks, such as those based on graphs (Gelman et al.,, 1996) or other
measures of surprise (Bayarri and Morales,, 2003).

Finally, Gelman et al., (1996) define mixed predictive checks. To form
replicated data in hierarchical models, they use the prior of group-specific
latent variables and the posterior for latent variables shared across groups.
This mitigates overfitting, but there are no guarantees. If a model has rich
enough latent variables shared across data, a mixed predictive check will
still overfit. Other predictive checks for hierarchical models include Dey et
al., (1998); Marshall and Spiegelhalter, (2003); Bayarri and Castellanos,
(2007).

##  3 Population Predictive Checks

We develop population predictive checks (pop-pcs). First, we set up notation
for Bayesian model checking and review the posterior predictive check (ppc).

###  3.1 Posterior Predictive Checks

A ppc formalizes the following intuition: “If my model is good (in a relevant
way) then its posterior predictive distribution will generate data that looks
like my observations (in that same way).” Rubin, (1984) designed the ppc as a
procedure to support the art of iterative model building, which requires tools
for checking and criticizing a model (Box,, 1980; Rubin,, 1984; Meng,, 1994;
Gelman et al.,, 1996).

Here is the recipe for a ppc. We have an observed data set and a model. We
simulate replicated data from the posterior predictive distribution, compute a
discrepancy that measures the “relevant ways” we want our model to capture the
data, and then compare the discrepancy on the observed data to the
distribution of the discrepancy on replicated data. (Proponents of the ppc s
prefer the term “discrepancy” to “test statistic” to veer away from language
that implies the model might be true.) One way to make such a comparison, for
example, is to estimate a tail probability, the probability that the
discrepancy of the replicated data is greater than the discrepancy of the
observed data.

More formally, consider a Bayesian model with hidden variables $\bm{\theta}$
and data $\mathbf{y}$. It defines a joint distribution as a product of the
prior and the likelihood,

$\displaystyle
p(\bm{\theta},\mathbf{y})=p(\bm{\theta})p(\mathbf{y}\,|\,\bm{\theta}).$ (1)

There are two types of data. The first is observed data
$\mathbf{y}^{\textrm{obs}}$. The second is replicated data
$\mathbf{y}^{\textrm{rep}}$. Replicated data are drawn from the posterior
predictive distribution,

$\displaystyle
p(\mathbf{y}^{\textrm{rep}}\,|\,\mathbf{y}^{\textrm{obs}})=\int_{\bm{\theta}}p(\mathbf{y}^{\textrm{rep}}\,|\,\bm{\theta})p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}}).$
(2)

The discrepancy $d(\mathbf{y})$ is a function of the data that measures
something we want to see replicated in the posterior predictive distribution.
For example, one discrepancy is the empirical mean,

$\displaystyle d(\mathbf{y})=\frac{1}{n}\sum_{i=1}^{n}y_{i}.$ (3)

With it, we can check whether the empirical mean of the observed data is
“typical” for data that come from the posterior predictive.

The ppc compares the distribution of discrepancies
$d(\mathbf{Y}^{\textrm{rep}})$ with the observed discrepancy
$d(\mathbf{y}^{\textrm{obs}})$. The original posterior predictive p-value is a
tail probability (Rubin,, 1984),

$\displaystyle
p_{\textsc{ppc}}=p(d(\mathbf{y}^{\textrm{rep}})>d(\mathbf{y}^{\textrm{obs}})\,|\,\mathbf{y}^{\textrm{obs}}),$
(4)

where the distribution is defined via the posterior predictive
$p(\mathbf{y}^{\textrm{rep}}\,|\,\mathbf{y}^{\textrm{obs}})$. In this paper,
we generalize this idea with some additional notation. Let $g(\cdot,\cdot)$ be
a distance function between discrepancies, also called a measure of surprise
(Bayarri and Morales,, 2003). A ppc is a posterior predictive expectation

$\displaystyle\textsc{ppc}(\mathbf{y}^{\textrm{obs}};g,d)=\mathrm{E}\left[g(d(\mathbf{Y}^{\textrm{rep}}),d(\mathbf{y}^{\textrm{obs}}))\,|\,\mathbf{y}^{\textrm{obs}}\right].$
(5)

The expectation is with respect to the distribution of
$\mathbf{Y}^{\textrm{rep}}$, which is the posterior predictive distribution.
In practice, this expectation can be approximated with Monte Carlo, where MCMC
or a variational approximation provides samples of
$\mathbf{Y}^{\textrm{rep}}$. For example, we recover the posterior predictive
p-value when the distance function is an indicator
$g(a,b)=\mathbb{1}\left[a>b\right]$. This generalization lets us consider
alternative distance functions and more flexible checks, such as those that
involve multivariate discrepancies.

An important innovation in Bayesian model checking was the realized
discrepancy, which is a function of both data and hidden variables
$d(\mathbf{y},\bm{\theta})$ (Meng,, 1994; Gelman et al.,, 1996). The realized
discrepancy assesses aspects of the relationship between data and the latent
variables. For example, one realized discrepancy is the log likelihood,

$\displaystyle d(\mathbf{y},\bm{\theta})=\sum_{i=1}^{n}\log
p(y_{i}\,|\,\bm{\theta}).$

Using a realized discrepancy, a ppc is

$\displaystyle\textsc{ppc}(\mathbf{y}^{\textrm{obs}};g,d)=\mathrm{E}\left[g(d(\mathbf{Y}^{\textrm{rep}},\bm{\theta}),d(\mathbf{y}^{\textrm{obs}},\bm{\theta}))\,|\,\mathbf{y}^{\textrm{obs}}\right].$
(6)

Here the expectation is with respect to a joint distribution of the latent
variables and replicated data
$p(\bm{\theta},\mathbf{y}^{\textrm{rep}}\,|\,\mathbf{y}^{\textrm{obs}})=p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})p(\mathbf{y}^{\textrm{rep}}\,|\,\bm{\theta})$.
Note the observed discrepancy $d(\mathbf{y}^{\textrm{obs}},\bm{\theta})$ is
random in its second argument; the replicated discrepancy
$d(\mathbf{Y}^{\textrm{rep}},\bm{\theta})$ is random in both arguments.

It is in the discrepancy, whether simple $d(\mathbf{y})$ or realized
$d(\mathbf{y},\bm{\theta})$, that the investigator defines what is important
about the model. Simple discrepancies connect the ppc to classical goodness-
of-fit tests (Gelman et al.,, 1996). Other discrepancies are defined by the
scientific community that use the models, e.g., in population genetics (Mimno
et al.,, 2015). Still more elaborate discrepancies can be multivariate (Lewis
and Raftery,, 1996). For example, the discrepancy might use the latent
variables to partition the observations and then measure a scalar discrepancy
on each group (Mimno and Blei,, 2011; Mimno et al.,, 2015).

We have described the ppc as an expectation, the calculation of a single
value. An alternative approach, and possibly one that is more widely used, is
to form visualizations that express the differences and similarities between
the observed and replicated data. With a simple discrepancy, one practice is
to plot the distribution of $d(\mathbf{y}^{\textrm{rep}})$ and then locate
$d(\mathbf{y}^{\textrm{obs}})$ within it. This plot is a nice companion to the
posterior predictive p-value; see, e.g., Belin and Rubin, (1995) for an
example from medical data. Realized discrepancies of both data and latent
variables require more intricate plots (Gelman et al.,, 1996; Mimno and Blei,,
2011; Mimno et al.,, 2015).

For $p$-values to be interpretable, their reference distribution should be
uniform (Robins et al.,, 2000). Posterior predictive $p$-values are not
interpretable in this way. They tend to place more mass around one-half than a
uniformly distributed variable, because simulations from the posterior
predictive distribution tend to be similar to the observations.

Double use of data. One of the central criticisms of posterior predictive
checks (ppcs) is that they use the data twice. The ppc first conditions on the
data to find the posterior predictive distribution and then uses the same data
to calculate the observed discrepancy; this is explicit in Equation 5 and
Equation 6. Though this property might be unsettling, especially to an
orthodox Bayesian, the dual use of the data is not technically wrong: the ppc
is a well-defined expectation that happens to use the data in two ways.

That said, an unfortunate consequence of this property is that it can lead to
a type of “overfitting” when we use a ppc to help build and refine a model.
The issue is that the ppc asks discrepancies of replicated data to be close to
the discrepancy of the same data that we condition on. To illustrate the
problem, consider a model whose posterior predictive distribution simply
memorizes the empirical distribution of the data. (Below we show that the
Dirichlet process can provide such an idiosyncratic model.) In this setting,
the observed discrepancy $d(\mathbf{y}^{\textrm{obs}})$ will be close to the
replicated discrepancy $d(\mathbf{y}^{\textrm{rep}})$, and so it will pass
most checks. But this goes against the spirit of expressive Bayesian modeling.
An investigator does not build a sequence of Bayesian models only to finally
decide on the empirical distribution of the data.

Of course, the overfitting problem does not render ppc’s useless. When used
responsibly—by maintaining simple models and using sensible discrepancies—they
are an effective tool for practical model building. But the double use of the
data can be a conceptual and practical pitfall of the method. To address this,
we expand the ppc methodology to use the idea of an unknown population
distribution, i.e., the true underlying distribution of the data. The
resulting procedure, which we call a pop-pc, provides a framework for model
checking that is inherently protected from overfitting.

###  3.2 Population Predictive Checks

The pop-pc changes the basic intuition of a model check to one that uses the
idea of a hypothetical population distribution: “If my model is good (in a
relevant way) then data drawn from the posterior predictive distribution will
look like data from the true population (in that way).” This shift leads to
techniques that naturally combine the ideas behind ppc’s with estimates of
out-of-sample measurements, such as through cross validation (Geisser,, 1975)
or the bootstrap (Efron,, 1979). We will show that pop-pc’s are more robust to
overfitting than ppc’s.

In addition to the usual ingredients of a ppc—the observed data
$\mathbf{y}^{\textrm{obs}}$ and replicated data
$\mathbf{y}^{\textrm{rep}}$—the pop-pc also uses new data
$\mathbf{y}^{\mathrm{new}}$, drawn from the true population distribution $F$.
(We do not assume the population is related to the considered model.) A pop-pc
uses the new data $\mathbf{y}^{\mathrm{new}}$ to check if the predictive
distribution induced by the observed data $\mathbf{y}^{\textrm{obs}}$ leads to
a distribution that is close to the population $F$. If so, this indicates that
the model’s posterior predictive captures future data well and gives reason to
investigate the posterior of the latent variables to help understand,
simplify, and summarize the observations.

As for a ppc, each check involves a discrepancy function $d(\cdot)$ and a
distance function $g(\cdot,\cdot)$. Recall that that distance function $g$
measures how close two discrepancies are. For now, we consider simple
discrepancies that only depend on data $d(\mathbf{y})$. A simple discrepancy
can apply to any type of data: observed data, replicated data from the
posterior predictive, or new data from the population.

With these quantities in place, we first define an ideal pop-pc. It compares
replicated data from the posterior predictive to new data from the true
population,

$\displaystyle\textsc{pop-
pc}_{\textrm{ideal}}(\mathbf{y}^{\textrm{obs}},F;d,g)=\mathrm{E}\left[g(d(\mathbf{Y}^{\textrm{rep}}),d(\mathbf{Y}^{\mathrm{new}}))\,|\,\mathbf{y}^{\textrm{obs}},F\right].$
(7)

The expectation is taken with respect to the joint distribution of replicated
data and future data,

$\displaystyle
p(\mathbf{y}^{\textrm{rep}},\mathbf{y}^{\mathrm{new}}\,|\,\mathbf{y}^{\textrm{obs}},F)=p(\mathbf{y}^{\textrm{rep}}\,|\,\mathbf{y}^{\textrm{obs}})p(\mathbf{y}^{\mathrm{new}}\,|\,F).$
(8)

Of course, we do not have access to the true population distribution or to new
data from it. In Section 4 we discuss procedures to approximate it in various
ways, connecting pop-pc’s to ideas such as the bootstrap (Efron,, 1979) and
cross validation (Stone,, 1974; Geisser,, 1975). But first we consider some of
the properties of an ideal pop-pc.

As an example, consider again the posterior predictive p-value. It is the tail
probability of Equation 4 and corresponds to a ppc with a particular
definition of $g(\cdot,\cdot)$. Like all ppc’s, it is widely criticized for
using the data twice (Gelfand et al.,, 1992; Bayarri and Berger,, 2000; Robins
et al.,, 2000) once to compare discrepancies and once to determine the
posterior predictive distribution. Again we emphasize that this “double
dipping” is not a mathematical error. Rather, it is a natural consequence of
the question being asked.

The analogous population predictive p-value asks a different question. It is

$\displaystyle
p_{\textrm{ideal}}(\mathbf{y}^{\textrm{obs}},F)=p(d(\mathbf{y}^{\textrm{rep}})>d(\mathbf{y}^{\mathrm{new}})\,|\,\mathbf{y}^{\textrm{obs}},F).$
(9)

The joint distribution of
$[\mathbf{y}^{\textrm{rep}},\mathbf{y}^{\mathrm{new}}]$ is in Equation 8. This
quantity no longer uses the data twice. It asks how often the discrepancy from
the posterior predictive is in the tail of the discrepancy from the true
population distribution.

To place these ideas in context, Figure 1 shows graphical diagram for
population predictive checks from this paper, posterior predictive checks
(Rubin,, 1984; Gelman et al.,, 1996), and prior predictive checks (Box,,
1980). Unlike prior predictive checks, posterior predictive checks and
population predictive checks condition on the observed data to generate
replicated data. The key difference between posterior and population checks is
in how each constructs the reference distribution, i.e., to what we compare
the posterior predictive distribution. Posterior predictive checks use the
observed data as reference; this is the root of the overfitting issue.
Population checks use new data from the true population distribution. This
mitigates the possibility of overfitting because data from the population is
different data from those that are conditioned on in the posterior predictive.

We now describe more details about the pop-pc. We present a toy example to
build intuition, describe how to incorporate realized discrepancies into this
framework, and describe how to check hierarchical Bayesian models with a pop-
pc.

###  3.3 A simple example of ppc s and pop-pc s

We demonstrate the pop-pc with a simple example, a Bayesian model for which
the posterior predictive distribution is effectively the empirical
distribution of the observations. Our example illuminates how a ppc can be led
astray by its double use of the data, and how the pop-pc is robust to such
issues.

Begin by defining the population distribution $F$ to be a Normal with mean 5
and variance 2. The population is unknown, but we observe a sequence of $n$
examples, $y_{i}\sim F$, $i\in\\{1,\ldots,n\\}$.

Now define an incorrect model. Our model uses the Dirichlet process (dp), a
Bayesian nonparametric prior on discrete random distributions (Ferguson,,
1973). Let $G$ be a random distribution drawn from a dp. The dp has two
parameters, a scalar $\alpha$ and the base measure $H$. In our example, we fix
$H$ to be the true population distribution. The draw $G$ will be a random
discrete distribution with the same support as $H$. Its atoms are distributed
by $H$ and it’s “peakiness” is determined by $\alpha$.

Using a dp, we assume

$\displaystyle G$ $\displaystyle\sim\mathrm{DP}(\alpha H)$ $\displaystyle
y_{i}$ $\displaystyle\sim G.$

For our purposes here, the details of the dp are not important. What is
important is that the corresponding predictive distribution is a mixture of
$H$ and the empirical distribution,

$\displaystyle
y\,|\,\mathbf{y}^{\textrm{obs}}\sim\frac{\alpha}{\alpha+n}H(y)+\sum_{i=1}^{n}\frac{1}{\alpha+n}\delta_{y^{\textrm{obs}}_{i}}(y).$
(10)

For small values of $\alpha$ (or large values of $n$), the posterior
predictive distribution is essentially the empirical distribution. In this
example, large values of $\alpha$ lead to a more correct model— the posterior
predictive more often draws from the true population $H$.

Suppose we want to check if replications from model have similar likelihood to
the data. The discrepancy is

$\displaystyle d(\mathbf{y})\triangleq\frac{1}{n}\sum_{i=1}^{n}\log
p(y_{i}\,|\,\mathbf{y}^{\textrm{obs}}).$ (11)

In the context of a ppc, this discrepancy is discussed in Lewis and Raftery,
(1996). Further, suppose we are interested in a $p$-value. We set $g$ to be
the indicator $g(a,b)=\mathbb{1}\left[a>b\right]$. For the ppc we compute
Equation 4; for the pop-pc we compute Equation 9.

Figure 2: The log-probability discrepancy as a function of the Dirichlet
process concentration parameter for both ppcs and pop-pcs. ppcs suggest the
most overfit model works best.

In the model, we set the base distribution to be the true population
distribution, a Gaussian, and observe a data set of ten points. We then vary
the scaling parameter $\alpha$ so that the posterior predictive ranges from
the empirical distribution ($\alpha$ low) to the true population distribution
($\alpha$ high); see Equation 10. Figure 2 compares posterior predictive
checks and (ideal) population predictive checks. On the x-axis is the log-
concentration parameter; on the y-axis is the p-value.

The behavior of the ideal pop-pc is how we would hope model-checking performs:
the p-value is extreme (zero) when there is model mismatch; it becomes less
extreme as the posterior predictive matches the data-generating distribution.
The ppc does not behave as well. It does not produce extreme p-values when the
posterior predictive is far from the generating distribution and, if compared
quantitatively, even slightly prefers the empirical distribution to the true
population. The ppc also has a bump when $\alpha$ is close to the observed
data size. The bump occurs because the replicated data samples from a balanced
mixture of the empirical and true population distribution which has lower
likelihood than only the observed data.

The reason the ppc fails is the double use of the data. When $\alpha$ is small
the posterior predictive is close to the empirical distribution. It
“memorizes” $\mathbf{y}^{\textrm{obs}}$ and then the ppc compares samples from
the memorized distribution to the same $\mathbf{y}^{\textrm{obs}}$; these two
sets of samples will always look similar. In contrast, the ideal pop-pc
compares its samples from the memorized distribution to data from the
population; it does not overfit.

This example illustrates the underlying issue with the ppc in the idealized
setup where we have access to new samples from the population distribution. In
real-world settings, we need to approximate sampling from the population
distribution using ideas like the bootstrap (Efron,, 1979). We discuss these
approximations in Section 4.

###  3.4 pop-pc s with realized discrepancies

We developed pop-pc s with a simple discrepancy, a function only of data.
However, many ppc s are built around a realized discrepancies (Gelman et al.,,
1996), functions of both latent variables and data. A realized discrepancy
measures the strength of the connection between latent variables and a data
set.

Consider a latent variable $\bm{\theta}$ and a data set $\mathbf{y}$. As one
example, Gelman et al., (1996) define the chi-squared discrepancy

$\displaystyle
d(\mathbf{y},\bm{\theta})=\sum_{i=1}^{N}\frac{(y_{i}-\mathrm{E}[Y_{i}|\bm{\theta}])^{2}}{\textrm{Var}(Y_{i}|\bm{\theta})}.$
(12)

Realized discrepancies admit more focused model checks than simple
discrepancies. For example, they can be used to check individual topics in a
topic model (Mimno and Blei,, 2011). In this section, we extend pop-pc s to
realized discrepancies.

We distinguish between global and local latent variables. A global latent
variable is one that helps govern the entire data set; this is the type of
latent variable we usually see in classical Bayesian statistics, such as a
parameter to a likelihood. A local latent variable is one that governs a
subset of data points; one example is a data point’s cluster assignment in a
mixture model. Figure 3 shows the graphical models for both a global latent
variable model and a global and local latent variable model. For now we focus
on global latent variables. We will discuss local latent variables below.

Figure 3: The left figure shows the graphical model with only a global latent
variable $\bm{\theta}$. The right model shows the graphical model with both a
global latent variable $\bm{\theta}$ and local latent variables $z_{i}$.

.

Realized discrepancies with global latent variables. We consider a generic
Bayesian model that contains global latent variables,

$\displaystyle\bm{\theta}$ $\displaystyle\sim p(\bm{\theta})$ $\displaystyle
y_{i}\,|\,\bm{\theta}$ $\displaystyle\sim p(y\,|\,\bm{\theta}).$

The prior $p(\bm{\theta})$ defines a distribution on the latent parameter; the
likelihood $p(\mathbf{y}\,|\,\bm{\theta})$ defines the conditional probability
of a data point.

To include global variables, we extend the ideal check of Section 3.2 to the
generalized ideal check. We define the discrepancy $d(\mathbf{y},\bm{\theta})$
as a function of both data and the latent variable (e.g., Equation 16). We
compute the corresponding check,

$\displaystyle\textsc{pop-
pc}_{\textrm{ideal}}(\mathbf{y}^{\textrm{obs}},F;d,g)=\mathrm{E}\left[g(d(\mathbf{Y}^{\textrm{rep}},\bm{\theta}),d(\mathbf{Y}^{\mathrm{new}},\bm{\theta}))\,|\,\mathbf{y}^{\textrm{obs}},F\right].$
(13)

Here, the sampling distribution is
$p(\mathbf{y}^{\textrm{rep}},\mathbf{y}^{\mathrm{new}},\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}},F)$.
We draw the global variable from the posterior, the replicated data from the
likelihood given the global variable, and the new data from the population,

$\displaystyle\bm{\theta}$ $\displaystyle\sim
p(\bm{\theta}|\mathbf{y}^{\textrm{obs}})$
$\displaystyle\mathbf{y}^{\textrm{rep}}$ $\displaystyle\sim
p(\mathbf{y}\,|\,\bm{\theta})$ $\displaystyle\mathbf{y}^{\mathrm{new}}$
$\displaystyle\sim F.$

This generalizes the ideal checks of Section 3.2, where the sampling
distribution (Equation 8) implicitly marginalizes out $\bm{\theta}$ because
the discrepancy does not depend on it.

The discrepancy $d(\mathbf{y},\bm{\theta})$ is meant to measure how well a
latent variable $\bm{\theta}$ “explains” a data set $\mathbf{y}$. Consider a
latent variable $\bm{\theta}$ from the posterior
$p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})$. Replicated data
$\mathbf{y}^{\textrm{rep}}$ is drawn directly from $\bm{\theta}$, which should
be explained well by $\bm{\theta}$. However, new data
$\mathbf{y}^{\mathrm{new}}$ is drawn from the true population. It is only
well-explained by $\bm{\theta}$ when the model is good, i.e., when the
posterior distribution produces latent variables that govern realistic data.
When a model passes such an ideal check, we are justified in interpreting the
latent variable as being indicative of a type of structure in the true
population.

Realized discrepancies with global and local latent variables. Along with
global variables, many models also contain local variables that help govern
individual groups of data. When we use a simple discrepancy (Equation 7) or a
realized discrepancy that only involves global variables (Equation 13), the
local variables can be marginalized out of the likelihood to form the sampling
distribution. Some forms of model criticism, however, use discrepancies that
depend on local variables too. As examples, Mimno and Blei, (2011) and Mimno
et al., (2015) use local variables to partition data and compute a multi-
dimensional discrepancy.

As an example of local variables, consider a model of grouped data. Denote the
data to be $y_{ij}$, the $i$th data point in the $j$th group. It is governed
in part by the global latent variable $\bm{\theta}$ and in part by its group-
wide local variable $z_{j}$. The model is

$\displaystyle\bm{\theta}$ $\displaystyle\sim p(\bm{\theta})$
$\displaystyle\mathbf{z}_{j}\,|\,\bm{\theta}$ $\displaystyle\sim
p(\mathbf{z}\,|\,\bm{\theta})$ $\displaystyle y_{ij}\,|\,z_{j},\bm{\theta}$
$\displaystyle\sim p(y\,|\,\mathbf{z}_{j},\bm{\theta})$

This outlines a large class of models, including probabilistic matrix
factorization (Salakhutdinov and Mnih,, 2008), mixed-membership models
(Airoldi et al.,, 2014), hierarchical regression models (Gelman and Hill,,
2007), and many Bayesian nonparametric models (Hjort et al.,, 2010). In these
models, all data depend on $\theta$ and each data group also depends on its
local variable.

With grouped data, we can imagine a per-group check: “If my model is good (in
a relevant way) then data drawn from the posterior predictive distribution of
a group will look like data from the true population of that group.” The
corresponding discrepancy is $d(\mathbf{y}_{j},\mathbf{z}_{j},\bm{\theta})$.
Notice that using this discrepancy requires new local data from the population
distribution of the $j$th group $F_{j}$. Per-group pop-pc s check how well new
data in a group match replications from the posterior predictive distribution
of that group.

Working with local variables and grouped data is conceptually the same as in
the simpler generalized ideal check of the previous section. However, the per-
group check requires two new quantities, $\mathbf{z}^{\textrm{rep}}_{j}$ and
$\mathbf{z}^{\mathrm{new}}_{j}$, which are the local variables attached to the
replicated data and new data for the $j$th group. The check is

$\displaystyle\textsc{pop-
pc}_{\textrm{ideal}}^{j}(\mathbf{y}^{\textrm{obs}},F_{j};d,g)=\mathrm{E}[g(d(\mathbf{Y}^{\textrm{rep}}_{j},\mathbf{Z}^{\textrm{rep}}_{j},\bm{\theta}),d(\mathbf{Y}^{\mathrm{new}}_{j},\mathbf{Z}^{\mathrm{new}}_{j},\bm{\theta}))|F_{j},\mathbf{y}^{\textrm{obs}}].$
(14)

There is a nuance to defining this check, which points to the essential
difference between a ppc and a pop-pc. Equation 14 requires a distribution of
$[\bm{\theta},\mathbf{z}^{\textrm{rep}}_{j},\mathbf{y}^{\textrm{rep}}_{j},\mathbf{y}^{\mathrm{new}}_{j},\mathbf{z}^{\mathrm{new}}_{j}]$.
Intuitively, the replicated local variables $\mathbf{z}^{\textrm{rep}}$ should
help explain the replicated data $\mathbf{y}^{\textrm{rep}}_{j}$ and the new
local variables $\mathbf{z}^{\mathrm{new}}_{j}$ should help explain the new
data $\mathbf{y}^{\mathrm{new}}_{j}$. (This mirrors the role of the global
latent variables in the discrepancy.) This is straightforward for the
replicated data set—it comes from the (posterior) model and thus naturally
includes $\mathbf{z}^{\textrm{rep}}_{j}$ in its generative process. But the
new data comes from the population distribution, which does not involve any
local latent variables, and we cannot assume that the local latent variables
for $\mathbf{y}^{\textrm{rep}}_{j}$ also explain the data from the population
distribution.

Thus we use the following sampling distribution

$\displaystyle\bm{\theta}$ $\displaystyle\sim
p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})$
$\displaystyle\mathbf{z}^{\textrm{rep}}_{j}$ $\displaystyle\sim
p(\mathbf{z}\,|\,\bm{\theta})$ $\displaystyle\mathbf{z}^{\mathrm{new}}_{j}$
$\displaystyle\sim
p(\mathbf{z}\,|\,\mathbf{y}^{\textrm{obs}}_{j},\bm{\theta})$
$\displaystyle\mathbf{y}^{\textrm{rep}}_{j}$ $\displaystyle\sim
p(\mathbf{y}\,|\,\mathbf{z}^{\textrm{rep}}_{j},\bm{\theta})$
$\displaystyle\mathbf{y}^{\textrm{obs}}_{j}$ $\displaystyle\sim F_{j}$
$\displaystyle\mathbf{y}^{\mathrm{new}}_{j}$ $\displaystyle\sim F_{j}$

The distribution of $\mathbf{z}^{\mathrm{new}}$ comes from a posterior
$p(\mathbf{z}\,|\,\mathbf{y}^{\textrm{obs}}_{j},\bm{\theta})$ conditioned on
the $j$th group. (Because of the structure of the model, conditional on
$\bm{\theta}$, the posterior of $\mathbf{z}^{\mathrm{new}}_{j}$ is independent
of the other data.) This choice captures that in a good model
$\mathbf{z}^{\mathrm{new}}_{j}$, drawn from the posterior, should help explain
$\mathbf{y}^{\mathrm{new}}_{j}$, drawn from the population. In fact, this can
be seen as a generalized ideal check conditional on $\bm{\theta}$. Notice that
we can use this check to criticize entirely new groups (where there is no
$\mathbf{y}^{\textrm{obs}}_{j})$ as well as “future data” from existing
groups.

To use a per-group check in an omnibus measure, we can average across multiple
groups drawn from the population. Imagine a hypothetical super population
distribution $F$ from which the subpopulations are drawn $F_{j}\sim F$, then
we can define an omnibus check as

$\displaystyle\textsc{pop-
pc}_{\textrm{ideal}}(\mathbf{y}^{\textrm{obs}};d,g)=\mathrm{E}_{F_{j}\sim
F}[\textsc{pop-pc}_{\textrm{ideal}}^{j}|\mathbf{y}^{\textrm{obs}}].$ (15)

Alternatively, we can define other discrepancies that involve multiple groups.

Note this construction extends naturally to hierarchical models with more than
two levels and to other conditional probability structures, including those
not organized in a tree.

##  4 Computing Population Predictive Checks

Population predictive checks require both an expectation with respect to the
posterior distribution and an expectation with respect to the true population
distribution. Recall the form of the ideal check with realized discrepancies,

$\displaystyle\textsc{pop-
pc}_{\textrm{ideal}}(\mathbf{y}^{\textrm{obs}},F;d,g)=\mathrm{E}\left[g(d(\mathbf{Y}^{\textrm{rep}},\bm{\theta}),d(\mathbf{Y}^{\mathrm{new}},\bm{\theta}))\,|\,\mathbf{y}^{\textrm{obs}},F\right].$

Approximating this expectation requires samples from the posterior
distribution, the posterior predictive distribution, and the population
distribution.

Sampling from the posterior and posterior predictive is straightforward. First
sample a realization of the latent variables $\bm{\theta}$ from an
approximation to the posterior $p(\bm{\theta}\,|\,\mathbf{y}^{\textrm{obs}})$,
e.g., by using Markov chain Monte Carlo or variational inference. Then sample
replicated data $\mathbf{y}^{\textrm{rep}}$ from
$p(\mathbf{y}\,|\,\bm{\theta})$. This procedure is the same both for ppc s and
pop-pc s.

What differentiates pop-pc s from ppc s is that pop-pc s also require
$\mathbf{y}^{\mathrm{new}}$, which are samples from the population
distribution. We will develop various strategies for obtaining such samples.

In some circumstances, a statistical inference will involve a stream of data
where new observations from the population distribution continuously appear.
Examples include user behavior data from online services (Twitter, Facebook,
Google) or purchase data from a market place. In these settings, we can
approximate the pop-pc directly with Monte Carlo. This also extends to
scientific applications, where we can obtain new samples from running more
replications of the same experiment. These scenarios allow use of the ideal
check.

In most applications of Bayesian analysis, however, new samples from the
population distribution are not readily available; we only have access to a
fixed data set $\mathbf{y}$. For these applications we will study different
ways to approximate the population expectation required for a pop-pc. At a
high level, all of our approaches approximate the population distribution with
a form of an empirical distribution, thus connecting Bayesian model checking
to traditional frequentist ideas like the bootstrap and the plug-in principle
(Efron,, 1979; Efron and Tibshirani,, 1993). In a pop-pc, these methods
produce both the data we condition on $\mathbf{y}^{\textrm{obs}}$ and the
“future data” $\mathbf{y}^{\mathrm{new}}$ as functions of $\mathbf{y}$.

Algorithm 1 outlines the general procedure for approximating a pop-pc. This
algorithm requires methods for sampling $\mathbf{y}^{\mathrm{new}}$ and
$\mathbf{y}^{\textrm{obs}}$ from $\mathbf{y}$, the fixed pool of observed data
points. Using $\mathbf{y}^{\textrm{obs}}$, it uses a posterior approximation
to provide samples of $\mathbf{y}^{\textrm{rep}}$. Finally, it computes the
distance between the discrepancy on $\mathbf{y}^{\textrm{rep}}$ and
$\mathbf{y}^{\mathrm{new}}$. We now detail several strategies for sampling
$\mathbf{y}^{\mathrm{new}}$ and $\mathbf{y}^{\textrm{obs}}$ from a fixed pool
of observations $\mathbf{y}$.

Cross validation. One strategy is to approximate the distribution of future
data with the knowledge that we will not see the same points as we did in the
observed dataset. This line of thinking naturally leads to cross validation
(Geisser,, 1975).

Let $\textrm{Emp}(\mathbf{y})=\frac{1}{n}\sum_{i=1}^{n}\delta_{y_{i}}$ be the
empirical distribution on the points in $\mathbf{y}$ and let
$\textrm{Emp}^{m}_{\textrm{no-replace}}(\mathbf{y})$ be the distribution of
$m$ points from the empirical distribution without replacement and define
$\mathcal{A}\setminus\mathcal{B}$ to be the data points in $\mathcal{A}$ not
in $\mathcal{B}$. In a cross-validated check, we sample
$\mathbf{y}^{\textrm{obs}}$ and $\mathbf{y}^{\mathrm{new}}$ as

$\displaystyle\mathbf{y}^{\textrm{obs}}_{\rm{cv}}$
$\displaystyle\sim\textrm{Emp}^{m}_{\rm{no-replace}}(\mathbf{y})$
$\displaystyle\mathbf{y}^{\mathrm{new}}_{\rm{cv}}$
$\displaystyle=\mathbf{y}\setminus\mathbf{y}^{\textrm{obs}}_{\rm{cv}}.$

Here $m$ (less than $n$) points are chosen directly for
$\mathbf{y}^{\textrm{obs}}$ without duplication and
$\mathbf{y}^{\mathrm{new}}$ takes the remainder. In this setting, pop-pcs are
a form of cross validated check (Gelfand et al.,, 1992; Marshall and
Spiegelhalter,, 2003).

Bootstrap and out-of-bag samples. The observations $\mathbf{y}^{\textrm{obs}}$
in cross-validation must contain fewer points than the given data because
cross-validation samples without replacement. Alternatively, we can sample
with replacement—this is a bootstrap sample. Let
$\textrm{Emp}^{m}(\mathbf{y})$ be the distribution of $m$ iid samples from the
empirical distribution. We obtain $\mathbf{y}^{\textrm{obs}}$ and
$\mathbf{y}^{\mathrm{new}}$ as

$\displaystyle\mathbf{y}^{\textrm{obs}}_{\rm{boot}}$
$\displaystyle\sim\textrm{Emp}^{n}(\mathbf{y})$
$\displaystyle\mathbf{y}^{\mathrm{new}}_{\rm{oobag}}$
$\displaystyle=\mathbf{y}\setminus\mathbf{y}^{\textrm{obs}}_{\rm{boot}}.$

As for cross-validation, the selection of $\mathbf{y}^{\mathrm{new}}$ is
deterministic given the unique points in $\mathbf{y}^{\textrm{obs}}$—it is
called the “out-of-bag” sample (Breiman,, 1996). On average
$\mathbf{y}^{\mathrm{new}}_{\rm{oobag}}$ will have 37% of the points of
$\mathbf{y}$. Alternatively, $\mathbf{y}^{\mathrm{new}}$ may also be randomly
sampled from the empirical distribution of the points not in
$\mathbf{y}^{\textrm{obs}}$.

The double bootstrap. Both cross-validation and the out-of-bag bootstrap
enforce that $\mathbf{y}^{\textrm{obs}}$ and $\mathbf{y}^{\mathrm{new}}$ be
disjoint. As with other bootstrap-based estimations of generalization error,
this can be too conservative. It overestimates the error or model
misspecification both because the new data is “too different” from the
observed data and because we have removed some of the information in the
observed data to create the new data (Efron and Tibshirani,, 1997).

The double bootstap method draws both $\mathbf{y}^{\textrm{obs}}$ and
$\mathbf{y}^{\mathrm{new}}$ with replacement from the empirical distribution
of $\mathbf{y}$. This is a direct application of the plug-in principle,

$\displaystyle\mathbf{y}^{\textrm{obs}}_{\rm{boot}}$
$\displaystyle\sim\textrm{Emp}^{n}(\mathbf{y})$
$\displaystyle\mathbf{y}^{\mathrm{new}}_{\rm{boot}}$
$\displaystyle\sim\textrm{Emp}^{n}(\mathbf{y}).$

For large $n$, the amount of overlap between
$\mathbf{y}^{\mathrm{new}}_{\rm{boot}}$ and
$\mathbf{y}^{\textrm{obs}}_{\rm{boot}}$ will be 37%. The new data
$\mathbf{y}^{\mathrm{new}}$ will contain both points that are and are not in
$\mathbf{y}^{\textrm{obs}}$.

The $p$-bootstrap. We can explicitly control the overlap between the two sets
by drawing $\mathbf{y}^{\mathrm{new}}$ from a mixture of points in
$\mathbf{y}^{\textrm{obs}}$ and the remaining points,

$\displaystyle\mathbf{y}^{\textrm{obs}}_{\rm{boot}}$
$\displaystyle\sim\textrm{Emp}^{n}(\mathbf{y})$
$\displaystyle\mathbf{y}^{\mathrm{new}}_{p-\rm{boot}}\sim(p\textrm{Emp}(\mathbf{y}-$
$\displaystyle\mathbf{y}^{\textrm{obs}}_{\rm{boot}})+(1-p)\textrm{Emp}(\mathbf{y}^{\textrm{obs}}_{\rm{boot}}))^{n}.$

This lets us control the relative number of unique points in
$\mathbf{y}^{\mathrm{new}}$ versus those points we also see in
$\mathbf{y}^{\textrm{obs}}$. This is the most general form of bootstrap
sampling distributions for $\mathbf{y}^{\textrm{obs}}$ and
$\mathbf{y}^{\mathrm{new}}$.

Discussion. When the estimation techniques we develop are used with
hierarchical models with local variables Figure 3, we get a form of
hierarchical cross validation and hierarchical bootstrap. There is both cross
validation done across groups for the shared latent variables $\bm{\theta}$
and cross validation within groups for $\mathbf{z}_{j}$, the group specific
latent variables.

Notice Algorithm 1 requires approximating the posterior for each new sample of
$\mathbf{y}^{\textrm{obs}}$. This computational expense can be mitigated using
bootstrap acceleration, such as the bag of little boostraps (Kleiner et al.,,
2012).

Input: The model $\log p(\mathbf{y},\theta)$, the data $\mathbf{y}$, number of
replications $R$

Output: A pop-pc

pop-pc = 0

for  _r < R_ do

Sample $\mathbf{y}^{\textrm{obs}}$ and $\mathbf{y}^{\mathrm{new}}$ from
$\mathbf{y}$ using the bootstrap, cross-validation, or another variant.

Sample $\theta$ from the (approximate) posterior
$p(\theta\,|\,\mathbf{y}^{\textrm{obs}})$.

Sample $\mathbf{y}^{\textrm{rep}}$ from $p(y\,|\,\theta)$. Sample to be the
same size as $\mathbf{y}^{\mathrm{new}}$.

$\lx@glossaries@gls@link{acronym}{PoPC}{\leavevmode\textsc{pop-
pc}}\leftarrow\lx@glossaries@gls@link{acronym}{PoPC}{\leavevmode\textsc{pop-
pc}}+g(d(\mathbf{y}^{\textrm{rep}},\theta),d(\mathbf{y}^{\mathrm{new}},\theta))$

end for

$\lx@glossaries@gls@link{acronym}{PoPC}{\leavevmode\textsc{pop-
pc}}\leftarrow\lx@glossaries@gls@link{acronym}{PoPC}{\leavevmode\textsc{pop-
pc}}/R$

Algorithm 1 Computing Population Predictive Checks

##  5 Empirical Study

We study population predictive checks on a regression model and a hierarchical
model of documents.

###  5.1 Regression

Figure 4: A comparison between the posterior predictive check and the
population predictive check estimated using the bootstrap and the ideal pop-pc
on simulated data. We find that the PPC always prefers flexible models and
cannot measure mismatch in the flexible models. The bootstrap estimated pop-pc
follows the same trends as the idealized check.

Consider Bayesian linear regression, a model for real valued data conditional
on covariates. With a Gaussian prior on the regression weights with variance
$c$, this model is

$\displaystyle\bm{\theta}$ $\displaystyle\sim\textrm{Normal}(0,c)$
$\displaystyle y_{i}$
$\displaystyle\sim\textrm{Normal}(\bm{\theta}^{\top}\mathbf{x}_{i},1).$

Working with Bayesian linear regression requires setting $c$. We compare pop-
pc and ppc as a function of $c$ on a simulation using the means squared error
discrepancy

$\displaystyle
d(\mathbf{y},\bm{\theta})=\frac{1}{N}\sum_{i=1}^{N}(y_{i}-\mathrm{E}[Y_{i}|\bm{\theta}])^{2},$
(16)

with the absolute value distance function $g(a,b)=|a-b|$.

For the simulation, we drew the covariates and observations as

$\displaystyle\mathbf{x}_{i}$ $\displaystyle\sim\textrm{Uniform}(0,1)$
$\displaystyle y_{i}$
$\displaystyle\sim\textrm{Normal}(\bm{\theta}^{\top}\mathbf{x}_{i},1),$

where $\bm{\theta}$ was generated as the product of a Student-$t$ and unit
interval uniform random variables. The simulated data consists of 100
covariates and 50 observations.

We compare the ppc to the ideal and bootstrap estimated pop-pc on this
simulation in Figure 4 (as this is a simulation we can compute the ideal
check). We find that ppcs get smaller as the variance increases (i.e., the
model becomes more flexible), while the pop-pcs detects that as the variance
increases the model overfits. The bootstrap pop-pc slightly overestimates the
ideal pop-pc, but follows the trend of the ideal pop-pc across variances.

We compare different sampling schemes for $\mathbf{y}^{\textrm{obs}}$ and
$\mathbf{y}^{\mathrm{new}}$ detailed in Section 4. For each sampling scheme,
Figure 5 plots the check values versus the log-prior variance. The 632
bootstrap estimator is the $p$-bootstrap estimator with $p=.632$. It has the
closest discrepancy to the ideal pop-pc across variances. All but the double
bootstrap overestimate the the ideal pop-pc discrepancy. Unlike the ppc all
estimators of the pop-pc detect the growing error with respect to the
population.

Figure 5: A comparison of various population predictive check estimation
techniques alongside the ideal check on simulated regression data. We find the
632 Bootstrap estimator to be closes to the ideal check, while all other
methods (except the double bootstrap) overestimate the discrepancy.

###  5.2 Hierarchical Model of Text

We next study pop-pcs on latent Dirichlet allocation (lda) (Blei et al.,,
2003), a hierarchical model of documents. lda models documents by first
drawing a collection of distributions over words, called topics denoted
$\bm{\theta}$, then drawing the words in each document by drawing the topic
proportions for that document and for each word selecting a topic from the
drawn topic proportions and drawing a word from that topic distribution. Its
generative process is

$\displaystyle\bm{\theta}$ $\displaystyle\sim\textrm{Dirichlet}(\eta)$
$\displaystyle\mathbf{z}_{j}$ $\displaystyle\sim\textrm{Dirichlet}(\alpha)$
$\displaystyle w_{k,j}$ $\displaystyle\sim\textrm{Categorical}(z_{j})$
$\displaystyle y_{k,j}$
$\displaystyle\sim\textrm{Categorical}(\bm{\theta}_{w_{k,j}})$

In latent Dirichlet allocation the groups $j$ are the documents and the
parameters shared across groups are the topics $\theta$.

We study the instantaneous mutual information (imi) discrepancy developed to
assess the fitness of each individual topic (Mimno and Blei,, 2011). The imi
discrepancy measures the independence between words and document ids
conditioned on the topic. Let $H$ be the entropy, then mathematically the imi
between a word $w$ and document id $d$ conditional on a topic $k$ is

$\displaystyle\textrm{IMI}(w,d\,|\,k)=H(d\,|\,k)-H(d\,|\,k,w).$

This is a multivariate realized discrepancy with local variables.

lda assumes that conditioned on the topic of a word, the word and document in
which it appears are independent. If the model fits well, the instantaneous
mutual information should be small as $d$ is conditionally independent of $w$
given $k$.

We compare ppcs to pop-pcs on a corpus of 100,000 news documents with a
vocabulary of 1,000 from _the New York Times_. We use out of bag evaluation at
both within and across documents where half of each document forms the
$\mathbf{y}^{\mathrm{new}}_{j}$ for that document and a fixed collection of
1,000 documents form $\mathbf{y}^{\mathrm{new}}$. We set the Dirichlet
hyperparameter on the topics to 0.1 on both the topics and the document
proportions. For inference, we do stochastic variational inference (Hoffman et
al.,, 2013) with minibatches of size 1,000. Following the original use of this
discrepancy, we compute the discrepancy conditioned on a fixed set of topics
shared across groups (Mimno and Blei,, 2011). This focuses our experiment on
comparing the ppc and the pop-pc local latent variables.

Figure 6 plots the deviance ratio between the ppc and the pop-pc as a function
of the number of topics in log-scale. As the number of topics grows large, the
model gets closer to memorizing the data on which its conditioned. We see that
population predictive checks show orders of magnitude more deviance than
posterior predictive checks when approaching this pathology.

Figure 6: The deviance of pop-pcs and ppcs as a function of the number of
topics. This model memorizes its conditioning data as the number of topics
gets large. The y-axis is in log-scale. The pop-pcs better detect this
pathology and outputs the least deviant topic is found when $K=10$.

##  6 Discussion

We developed population predictive checks (pop-pcs), a diagnostic tool that
brings together Bayesian methods for model checking with frequentist
estimation of goodness of fit.

The ideal pop-pc assesses a Bayesian model by comparing samples from the
posterior predictive to samples from the true data generating process, i.e.,
the population distribution. Typically, however new samples from the
population distribution are not available. Thus we developed several
methods—incorporating ideas such as cross validation and the bootstrap—for
using finite data to estimate the ideal pop-pc. We demonstrated the merits of
pop-pcs with Bayesian linear regression models and on probabilistic topic
models of documents.

pop-pcs improve on classical posterior predictive checks (Rubin,, 1984; Gelman
et al.,, 1996) and on prior predictive checks (Box,, 1980). Posterior
predictive checks use the data twice and, if employed without care, naturally
lead to overfitted models. Prior checks do not use the data twice, but focus
on checking the model’s marginal distribution of the data; applied Bayesian
data analysis is more concerned about fidelity of the posterior predictive.

Note that the ideal check of Equation 7 is “Bayesian justifiable,” in Rubin’s
sense, because it always conditions on the observed data. However, the pop-pc
is not—all methods for accounting for the unknown population distribution
require manipulating our sense of what is observed. Indeed, it is the
requirement that we always condition on $\mathbf{y}^{\textrm{obs}}$ that leads
to the double-use of the data in the posterior predictive check (ppc).

There are several avenues for further research in population-based model
checking. First, we can develop methods that handle non-exchangeable data,
such as time series or networks; approximating a population distribution in
those settings is a challenging problem. Second, we can consider alternative
estimators of population-level expectations. Examples include the bag of
little bootstraps for large data sets (Kleiner et al.,, 2012), or the leave-
one-out estimators of Vehtari et al., (2016). Finally, an important goal for
future work is to implement generic software tools for pop-pcs, e.g., in a
probabilistic programming system like Stan (Carpenter et al.,, 2016) or Edward
(Tran et al., 2016b, ). To date, probabilistic programming has focused on
language design and inference compilation; including methods for model
diagnostics is important to building a complete system for applied Bayesian
modeling.

## Acknowledgments

The authors thank Yixin Wang for providing helpful comments on this article
and Scott Linderman for help with Figure 1.

## References

  * Airoldi et al., (2014) Airoldi, E., Blei, D., Erosheva, E., and Fienberg, S., editors (2014).  Handbook of Mixed Membership Models and Their Applications.  CRC Press.
  * Angelino et al., (2016) Angelino, E., Johnson, M., and Adams, R. (2016).  Patterns of scalable Bayesian inference.  Foundations and Trends in Machine Learning, 9(2-3):119–247.
  * Bayarri and Berger, (2000) Bayarri, M. and Berger, J. (2000).  P-values for composite null models.  Journal of the American Statistical Association, 95(452):1127–1142.
  * Bayarri and Castellanos, (2007) Bayarri, M. and Castellanos, M. (2007).  Bayesian checking of the second levels of hierarchical models.  Statistical Science, 22:322–343.
  * Bayarri and Morales, (2003) Bayarri, M. and Morales, J. (2003).  Bayesian measures of surprise for outlier detection.  Journal of Statistical Planning and Inference, 111:3–22.
  * Beaumont et al., (2002) Beaumont, M., Zhang, W., and Balding, D. (2002).  Approximate Bayesian computation in population genetics.  Genetics, 162(4):2025–2035.
  * Belin and Rubin, (1995) Belin, T. and Rubin, D. (1995).  The analysis of repeated-measures data on schizophrenic reaction times using mixture models.  Statistics in Medicine, 14(8):747–768.
  * Blei, (2014) Blei, D. (2014).  Build, compute, critique, repeat: Data analysis with latent variable models.  Annual Review of Statistics and Its Application, 1:203–232.
  * Blei et al., (2017) Blei, D., Kucukelbir, A., and McAuliffe, J. (2017).  Variational inference: A review for statisticians.  Journal of American Statistical Association, 112(518):859–877.
  * Blei et al., (2003) Blei, D., Ng, A., and Jordan, M. (2003).  Latent Dirichlet allocation.  Journal of Machine Learning Research, 3:993–1022.
  * Box, (1980) Box, G. (1980).  Sampling and Bayes’ inference in scientific modeling and robustness.  Journal of the Royal Statistical Society, Series A, 143(4):383–430.
  * Breiman, (1996) Breiman, L. (1996).  Out-of-bag estimation.  Technical report, Statistics Department, U.C. Berkeley.
  * Carpenter et al., (2016) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2016).  Stan: A probabilistic programming language.  Journal of Statistical Software.
  * Chatfield, (1995) Chatfield, C. (1995).  Model uncertainty, data mining and statistical inference.  Journal of the Royal Statistical Society. Series A (Statistics in Society), 158(3):419–466.
  * Dahl et al., (2007) Dahl, F., Gåsemyr, J., and Natvig, B. (2007).  A robust conflict measure of inconsistencies in bayesian hierarchical models.  Scandinavian Journal of Statistics, 34(4):816–828.
  * Dey et al., (1998) Dey, D., Gelfand, A., Swartz, T., and Vlachos, P. (1998).  Simulation based model checking for hierarchical models.  Test.
  * Draper, (1996) Draper, D. (1996).  Comment: Utility,sensitivity analysis, and cross-validation in bayesian model-checking.  Statistica Sinica, 6(760–767.).
  * Efron, (1979) Efron, B. (1979).  Bootstrap methods: Another look at the jackknife.  The Annals of Statistics, pages 1–26.
  * Efron and Tibshirani, (1993) Efron, B. and Tibshirani, R. (1993).  An Introduction to the Bootstrap.  Chapman & Hall.
  * Efron and Tibshirani, (1997) Efron, B. and Tibshirani, R. (1997).  Improvements on cross-validation: The 632+ bootstrap method.  Journal of the American Statistical Association, 92(438):548–560.
  * Ferguson, (1973) Ferguson, T. (1973).  A Bayesian analysis of some nonparametric problems.  The Annals of Statistics, 1:209–230.
  * Geisser, (1975) Geisser, S. (1975).  The predictive sample reuse method with applications.  Journal of the American Statistical Association, 70(350):320–328.
  * Gelfand et al., (1992) Gelfand, A., Dey, D., and Chang, H. (1992).  Model determination using predictive distributions with implementation via sampling-based methods.  Bayesian Statistics, 4.
  * Gelman et al., (1995) Gelman, A., Carlin, J., Stern, H., and Rubin, D. (1995).  Bayesian Data Analysis.  Chapman & Hall, London.
  * Gelman and Hill, (2007) Gelman, A. and Hill, J. (2007).  Data Analysis Using Regression and Multilevel/Hierarchical Models.  Cambridge University Press.
  * Gelman et al., (1996) Gelman, A., Meng, X., and Stern, H. (1996).  Posterior predictive assessment of model fitness via realized discrepancies.  Statistica Sinica, 6:733–807.
  * Goodman and Stuhlmuller, (2015) Goodman, N. and Stuhlmuller, A. (2015).  The design and implementation of probabilistic programming languages.
  * Guttman, (1967) Guttman, I. (1967).  The use of the concept of a future observation in goodness-of-fit problems.  Journal of the Royal Statistical Society. Series B (Methodological), pages 83–100.
  * Hjort et al., (2006) Hjort, N., Dahl, F., and Steinbakk, G. (2006).  Post-processing posterior predictive p values.  Journal of the American Statistical Association, 101(475):1157–1174.
  * Hjort et al., (2010) Hjort, N., Holmes, C., Mueller, P., and Walker, S. (2010).  Bayesian Nonparametrics: Principles and Practice.  Cambridge University Press, Cambridge, UK.
  * Hoffman et al., (2013) Hoffman, M., Blei, D., Wang, C., and Paisley, J. (2013).  Stochastic variational inference.  Journal of Machine Learning Research, 14(1303–1347).
  * Hoffman and Gelman, (2014) Hoffman, M. and Gelman, A. (2014).  The no-U-turn sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(1):1593–1623.
  * Kleiner et al., (2012) Kleiner, A., Talwalkar, A., Sarkar, P., and Jordan, M. (2012).  A scalable bootstrap for massive data.  arXiv":1112.5016.
  * Lewis and Raftery, (1996) Lewis, S. and Raftery, A. (1996).  Comment: Posterior predictive assessment for data subsets in hierarchical models via MCMC.  Statistica Sinica, 6:779–786.
  * Mansinghka et al., (2014) Mansinghka, V., Selsam, D., and Perov, Y. (2014).  Venture: A higher-order probabilistic programming platform with programmable inference.  arXiv:1404.0099.
  * Marshall and Spiegelhalter, (2003) Marshall, E. and Spiegelhalter, D. (2003).  Approximate cross-validatory predictive checks in disease mapping models.  Statistics in Medicine, 22:1649–1660.
  * Meng, (1994) Meng, X.-L. (1994).  Posterior predictive p-values.  The Annals of Statistics, pages 1142–1160.
  * Mimno and Blei, (2011) Mimno, D. and Blei, D. (2011).  Bayesian checking for topic models.  In Empirical Methods in Natural Language Processing.
  * Mimno et al., (2015) Mimno, D., Blei, D., and Engelhardt, B. (2015).  Posterior predictive checks to quantify lack-of-fit in admixture models of latent population structure.  Proceedings of the National Academy of Sciences.
  * O’Hagan, (2003) O’Hagan, A. (2003).  HSSS model criticism.  In Highly Structured Stochastic Systems, volume 27 of Oxford Statist. Sci. Ser., pages 423–453. Oxford Univ. Press, Oxford.
  * Ranganath et al., (2014) Ranganath, R., Gerrish, S., and Blei, D. (2014).  Black box variational inference.  In Artificial Intelligence and Statistics, pages 814–822.
  * Robins et al., (2000) Robins, J., van der Vaart, A., and Ventura, V. (2000).  Asymptotic distribution of $p$-values in composite null models.  Journal of the American Statistical Association, 95(452):1143–1156.
  * Rubin, (1984) Rubin, D. (1984).  Bayesianly justifiable and relevant frequency calculations for the applied statistician.  The Annals of Statistics, 12(4):1151–1172.
  * Salakhutdinov and Mnih, (2008) Salakhutdinov, R. and Mnih, A. (2008).  Probabilistic matrix factorization.  In Neural Information Processing Systems.
  * Srivastava et al., (2018) Srivastava, S., Li, C., and Dunson, D. (2018).  Scalable Bayes via barycenter in Wasserstein space.  The Journal of Machine Learning Research, 19(1):312–346.
  * Stan Development Team, (2015) Stan Development Team (2015).  Stan: A C++ library for probability and sampling, version 2.8.0.
  * Stone, (1974) Stone, M. (1974).  Cross-validatory choice and assessment of statistical predictions.  Journal of the Royal Statistical Society. Series B (Methodological), pages 111–147.
  * (48) Tran, D., Kucukelbir, A., Dieng, A., Rudolph, M., Liang, D., and Blei, D. (2016a).  Edward: A library for probabilistic modeling, inference, and criticism.  arXiv:1610.09787.
  * (49) Tran, D., Kucukelbir, A., Dieng, A. B., Rudolph, M., Liang, D., and Blei, D. M. (2016b).  Edward: A library for probabilistic modeling, inference, and criticism.  arXiv:1610.09787.
  * Vehtari et al., (2016) Vehtari, A., Gelman, A., and Gabry, J. (2016).  Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC.  Statistics and Computing.
  * Wood et al., (2014) Wood, F., van de Meent, J.-W., and Mansinghka, V. (2014).  A new approach to probabilistic programming inference.  In Artificial Intelligence and Statistics, pages 1024–1032.