3

This question is the follow-up of this previous question: Bayesian inference and testable implications.

For concreteness, consider the following bayesian model. This model is not to be taken literally, it is simply suppose to stand for a model that cannot capture the DGP but we do not know that a priori. However, I very much would like an answer that takes this concrete model to perform an actual posterior predictive check, so we avoid generic answers.

$$ \text{Likelihood:}\\ \\ y \sim \mathcal{N}(\mu_1, \sigma_1)\\ x \sim \mathcal{N}(\mu_2, \sigma_2)\\ \text{Prior:}\\ \\ \mu_1 \sim \mathcal{N}(0, 1000)\\ a \sim \mathcal{U(0,2)}\\ \mu_2 \leftarrow \mu_1 + a\\ \sigma_1 \sim \mathcal{U}(0, 100)\\ \sigma_2 \sim \mathcal{U}(0, 100) $$

Where $\mathcal{N}()$ denotes a gaussian and $\mathcal{U}()$ denotes a uniform distribution. Here is the implementation in rjags:

library(rjags)
  model <- "
model {
  for (i in 1:length(x)){
    x[i] ~ dnorm(mu1, tau1)
  }

  for (i in 1:length(y)){
    y[i] ~ dnorm(mu2, tau2)
  }

  mu1 ~ dnorm(0, .00001)
  a ~ dunif(0, 2)
  mu2 <- mu1 + a

  sigma1 ~ dunif(0,100)
  tau1 <- pow(sigma1, -2)

  sigma2 ~ dunif(0,100)
  tau2 <- pow(sigma2, -2)
}
"

And here is the model fitted to some simulated data that does not conform to the model's assumptions.

n <- 10
dat <- list(x = rnorm(n, mean = 2, sd = 2),
            y = rnorm(n, mean = 10, sd = 10))

jags.model   <- jags.model(textConnection(model), data =dat)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 4
#>    Total graph size: 32
#> 
#> Initializing model
samp <- coda.samples(jags.model, n.iter = 1e4, 
                       variable.names = c("mu1", "mu2", "sigma1", "sigma2"))
post  <- as.data.frame(samp[[1]])
summary(post$mu1)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  -1.732   1.456   1.977   2.004   2.526   6.897
summary(post$mu2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -0.9573  2.4011  3.0740  3.0808  3.7376  8.2234

Now, how do I formally perform a "posterior predictive check" in this model with this data? And how do I formally decide, using the posterior predictive check, that the model misfit is "bad enough" so that I "reject" this model? What "test statistic" would you use? Which "threshold" for decision would you use? And so on. If there are missing details that are required for solving this problem (like, say, you need a cost or loss function) please feel free to add those details in your answer as needed; these details are part of a good answer, since they clarify what we need to know to actually perform the check.

Finally, please try to provide an actual solution to this toy problem. It doesn't need to be code, if you can derive the numerical results by hand that works as well. But the main idea is to have this toy problem actually solved.

  • A simple way to sample from the posterior predictive is to include a missing value in `x` and `y`. JAGS will automatically sample from the PP distribution. Set the monitoring on `x[11]` and `y[11]` (for a sample size of 10) to get the PP distribution for `x` and `y`. Then compare those with the actual values of `x` and `y`. That would be a posterior predictive check in this case, if I'm not mistaken. – COOLSerdash May 24 '20 at 07:45
  • @COOLSerdash thanks, do you want to show how to do this? Also, I'm not asking only about how to sample from the posterior, but also how to actually perform the "check". Which test statistic should you use? What threshold should you use to reject the model? Etc. –  May 24 '20 at 13:29
  • 1
    I think this is a reasonable question and don't quite understand the downvotes. But the request for an implementation is off-topic here, and I'd recommend you remove it. – mkt May 26 '20 at 14:51
  • @mkt-ReinstateMonica why is it off-topic? I have seen several answers where we cannot assess correctness or that are hand-wavy and not very useful. I am asking for implementation just as a sanity check, this model is simple enough that whoever knows what he is doing should be able to carry out. –  May 26 '20 at 14:52
  • 1
    @mkt-ReinstateMonica think of it as just a small cost to avoid those people who are tempted to give generic answers like "there are several ways to do it, you could do it like this, or like that." If the person knows how to do posterior predictive checks, it should be trivial to do it in this example. No? –  May 26 '20 at 14:53
  • Them's the rules, I don't make them. It's a bit of a grey area and so you may be okay, but you also risk getting this closed. Hence my recommendation; it's your call. – mkt May 26 '20 at 14:56
  • @mkt-ReinstateMonica ok thanks for the feedback. Is it the problem that people might think this is about the code? I don't care about the code, I just want to avoid generic answers. –  May 26 '20 at 14:57
  • Yes, it's about the code request. As for generic answers, the best you can do is ask a precise question and then vote/comment to encourage the type of answer you are after. The rest is up to readers and their judgement. I'll add that an unhelpful answer does not hurt you in any way, it will just fail to attract votes and get ignored. – mkt May 26 '20 at 14:59
  • 2
    @mkt-ReinstateMonica thanks I just reworded the question, hope it is a bit better. –  May 26 '20 at 15:03
  • 7
    Too lazy to construct an actual answer, but have you consulted Gelman's Bayesian Data Analysis? The PDF is available for free online, and chapter 6.3 is devoted to posterior predictive checks with some examples partially worked out – klumbard May 26 '20 at 15:43

0 Answers0