Graphical posterior predictive checks (PPCs)

Thebayesplot package provides various plottingfunctions forgraphical posterior predictive checking, that is,creating graphical displays comparing observed data to simulated datafrom the posterior predictive distribution (Gabryet al, 2019).

The idea behind posterior predictive checking is simple: if a modelis a good fit then we should be able to use it to generate data thatlooks a lot like the data we observed. To generate the data used forposterior predictive checks (PPCs) we simulate from theposteriorpredictive distribution. This is the distribution of the outcomevariable implied by a model after using the observed data\(y\) (a vector of\(N\) outcome values) to update our beliefsabout unknown model parameters\(\theta\). The posterior predictivedistribution for observation\(\widetilde{y}\) can be written as\[p(\widetilde{y} \,|\, y) = \intp(\widetilde{y} \,|\, \theta) \, p(\theta \,|\, y) \, d\theta.\]Typically we will also condition on\(X\) (a matrix of predictor variables).

For each draw (simulation)\(s = 1, \ldots,S\) of the parameters from the posterior distribution,\(\theta^{(s)} \sim p(\theta \,|\, y)\), wedraw an entire vector of\(N\) outcomes\(\widetilde{y}^{(s)}\) from theposterior predictive distribution by simulating from the data modelconditional on parameters\(\theta^{(s)}\). The result is an\(S \times N\) matrix of draws\(\widetilde{y}\).

When simulating from the posterior predictive distribution we can useeither the same values of the predictors\(X\) that we used when fitting the model ornew observations of those predictors. When we use the same values of\(X\) we denote the resultingsimulations by\(y^{rep}\), as they canbe thought of as replications of the outcome\(y\) rather than predictions for futureobservations (\(\widetilde{y}\) usingpredictors\(\widetilde{X}\)). Thiscorresponds to the notation from Gelman et al. (2013) and is thenotation used throughout the package documentation.

Using the replicated datasets drawn from the posterior predictivedistribution, the functions in thebayesplot packagecreate various graphical displays comparing the observed data\(y\) to the replications. The names of thebayesplot plotting functions for posterior predictivechecking all have the prefixppc_.

Setup

In addition tobayesplot we’ll load the followingpackages:

• ggplot2, in case we want to customize the ggplotobjects created bybayesplot
• rstanarm, for fitting the example models usedthroughout the vignette

library("bayesplot")library("ggplot2")library("rstanarm")

Example models

To demonstrate some of the various PPCs that can be created with thebayesplot package we’ll use an example of comparingPoisson and Negative binomial regression models from one of therstanarmpackagevignettes (Gabry and Goodrich, 2017).

We want to make inferences about the efficacy of a certain pestmanagement system at reducing the number of roaches in urban apartments.[…] The regression predictors for the model are the pre-treatment numberof roachesroach1, the treatment indicatortreatment, and a variablesenior indicatingwhether the apartment is in a building restricted to elderly residents.Because the number of days for which the roach traps were used is notthe same for all apartments in the sample, we include it as an exposure[…].

First we fit a Poisson regression model with outcome variabley representing the roach count in each apartment at the endof the experiment.

head(roaches)# see help("rstanarm-datasets")

    y roach1 treatment senior exposure21 153 308.00         1      0  0.8000002 127 331.25         1      0  0.6000003   7   1.67         1      0  1.0000004   7   3.00         1      0  1.0000005   0   2.00         1      0  1.1428576   0   0.00         1      0  1.000000

roaches$roach100<- roaches$roach1/100# pre-treatment number of roaches (in 100s)

# using rstanarm's default priors. For details see the section on default# weakly informative priors at https://mc-stan.org/rstanarm/articles/priors.htmlfit_poisson<-stan_glm(  y~ roach100+ treatment+ senior,offset =log(exposure2),family =poisson(link ="log"),data = roaches,seed =1111,refresh =0# suppresses all output as of v2.18.1 of rstan)

print(fit_poisson)

stan_glm family:       poisson [log] formula:      y ~ roach100 + treatment + senior observations: 262 predictors:   4------            Median MAD_SD(Intercept)  3.1    0.0  roach100     0.7    0.0  treatment   -0.5    0.0  senior      -0.4    0.0  ------* For help interpreting the printed output see ?print.stanreg* For info on the priors used see ?prior_summary.stanreg

We’ll also fit the negative binomial model that we’ll compare to thePoisson:

fit_nb<-update(fit_poisson,family ="neg_binomial_2")

print(fit_nb)

stan_glm family:       neg_binomial_2 [log] formula:      y ~ roach100 + treatment + senior observations: 262 predictors:   4------            Median MAD_SD(Intercept)  2.8    0.2  roach100     1.3    0.3  treatment   -0.8    0.2  senior      -0.3    0.3  Auxiliary parameter(s):                      Median MAD_SDreciprocal_dispersion 0.3    0.0   ------* For help interpreting the printed output see ?print.stanreg* For info on the priors used see ?prior_summary.stanreg

Definingy andyrep

In order to use the PPC functions from thebayesplotpackage we need a vectory of outcome values,

y<- roaches$y

and a matrixyrep of draws from the posterior predictivedistribution,

yrep_poisson<-posterior_predict(fit_poisson,draws =500)yrep_nb<-posterior_predict(fit_nb,draws =500)dim(yrep_poisson)

[1] 500 262

dim(yrep_nb)

[1] 500 262

Each row of the matrix is a draw from the posterior predictivedistribution, i.e. a vector with one element for each of the data pointsiny.

Since we fit the models usingrstanarm we used itsspecialposterior_predict function, but if we were using amodel fit with therstan package we could createyrep in thegenerated quantities block of theStan program or by doing simulations in R after fitting the model. Drawsfrom the posterior predictive distribution can be used withbayesplot regardless of whether or not the model wasfit using an interface to Stan.bayesplot just requiresayrep matrix that hasnumber_of_draws rowsandnumber_of_observations columns.

ppc_dens_overlay

The first PPC we’ll look at is a comparison of the distribution ofy and the distributions of some of the simulated datasets(rows) in theyrep matrix.

color_scheme_set("brightblue")ppc_dens_overlay(y, yrep_poisson[1:50, ])

In the plot above, the dark line is the distribution of the observedoutcomesy and each of the 50 lighter lines is the kerneldensity estimate of one of the replications ofy from theposterior predictive distribution (i.e., one of the rows inyrep). This plot makes it easy to see that this model failsto account for the large proportion of zeros iny. That is,the model predicts fewer zeros than were actually observed.

To see the discrepancy at the lower values of more clearly we can usethexlim function fromggplot2 to restrictthe range of the x-axis:

ppc_dens_overlay(y, yrep_poisson[1:50, ])+xlim(0,150)

See Figure 6 inGabry et al. (2019) foranother example of usingppc_dens_overlay.

ppc_hist

We could see the same thing from a different perspective by lookingat separate histograms ofy and some of theyrep datasets using theppc_hist function:

ppc_hist(y, yrep_poisson[1:5, ])

The same plot for the negative binomial model looks muchdifferent:

ppc_hist(y, yrep_nb[1:5, ])

The negative binomial model does better handling the number of zerosin the data, but it occasionally predicts values that are way too large,which is why the x-axes extend to such high values in the plot and makeit difficult to read. To see the predictions for the smaller values moreclearly we can zoom in:

ppc_hist(y, yrep_nb[1:5, ],binwidth =20)+coord_cartesian(xlim =c(-1,300))

ppc_stat

First we define a function that takes a vector as input and returnsthe proportion of zeros:

prop_zero<-function(x)mean(x==0)prop_zero(y)# check proportion of zeros in y

[1] 0.3587786

Thestat argument toppc_stat accepts afunction or the name of a function for computing a test statistic from avector of data. In our case we can specifystat = "prop_zero" since we’ve already defined theprop_zero function, but we also could have usedstat = function(x) mean(x == 0).

ppc_stat(y, yrep_poisson,stat ="prop_zero",binwidth =0.005)

The dark line is at the value\(T(y)\), i.e. the value of the teststatistic computed from the observed\(y\), in this caseprop_zero(y). The lighter area on the left is actually ahistogram of the proportion of zeros in in theyrepsimulations, but it can be hard to see because almost none of thesimulated datasets inyrep have any zeros.

Here’s the same plot for the negative binomial model:

ppc_stat(y, yrep_nb,stat ="prop_zero")

Again we see that the negative binomial model does a much better jobpredicting the proportion of observed zeros than the Poisson.

However, if we look instead at the distribution of the maximum valuein the replications, we can see that the Poisson model makes morerealistic predictions than the negative binomial:

ppc_stat(y, yrep_poisson,stat ="max")

ppc_stat(y, yrep_nb,stat ="max")

ppc_stat(y, yrep_nb,stat ="max",binwidth =100)+coord_cartesian(xlim =c(-1,5000))

See Figure 7 inGabry et al. (2019) foranother example of usingppc_stat.

available_ppc()

available_ppc(pattern ="_grouped")

ppc_stat_grouped

For example,ppc_stat_grouped is the same asppc_stat except that the test statistic is computed withinlevels of the grouping variable and a separate plot is made for eachlevel:

ppc_stat_grouped(y, yrep_nb,group = roaches$treatment,stat ="prop_zero")

See Figure 8 inGabry et al. (2019) foranother example of usingppc_stat_grouped.

Defining app_check method

Here is an example for how to define a simplepp_checkmethod in a package that creates fitted model objects of class"foo". We will define a methodpp_check.foothat extracts the datay and the draws from the posteriorpredictive distributionyrep from an object of class"foo" and then calls one of the plotting functions frombayesplot.

Suppose that objects of class"foo" are lists with namedcomponents, two of which arey andyrep.Here’s a simple methodpp_check.foo that offers the userthe option of two different plots:

# @param object An object of class "foo".# @param type The type of plot.# @param ... Optional arguments passed on to the bayesplot plotting function.pp_check.foo<-function(object,type =c("multiple","overlaid"), ...) {  type<-match.arg(type)  y<- object[["y"]]  yrep<- object[["yrep"]]stopifnot(nrow(yrep)>=50)  samp<-sample(nrow(yrep),size =ifelse(type=="overlaid",50,5))  yrep<- yrep[samp, ]if (type=="overlaid") {ppc_dens_overlay(y, yrep, ...)  }else {ppc_hist(y, yrep, ...)  }}

To try outpp_check.foo we can just make a list withy andyrep components and give it classfoo:

x<-list(y =rnorm(200),yrep =matrix(rnorm(1e5),nrow =500,ncol =200))class(x)<-"foo"

color_scheme_set("purple")pp_check(x,type ="multiple",binwidth =0.3)

color_scheme_set("darkgray")pp_check(x,type ="overlaid")

Examples ofpp_check methods in other packages

Several packages currently use this approach to provide an interfacetobayesplot’s graphical posterior predictive checks.See, for example, thepp_check methods in therstanarmandbrmspackages.