Printing the posterior predictive p-value with ppc_stat

[dkaplan]

Hi all,

Does anyone know a convenient way to obtain the posterior predictive p-value and have it printed along with the ppc_stat plot?

Thanks in advance,

David

[jonah]

Hey David,

Here’s an example using stat=median but you can do the same with any stat:

library(bayesplot)
library(ggplot2)


# just using example objects that come with bayesplot
y <- example_y_data()
yrep <- example_yrep_draws()

# using stat=median but doesn't matter which stat
plot <- ppc_stat(y, yrep, stat = "median")

# calculate proportion of stat(yrep) > stat(y)
p <- mean(apply(yrep, 1, median) > median(y))
plot + 
  yaxis_text() + # just so I can see y-axis values for specifying them in annotate() below, but can remove this if you don't want the useless y-axis values displayed 
  annotate("text", x = 89, y = 40, label = paste("p =", p))

Here’s the plot this creates:

[dkaplan]

Perfect. Thanks!