# Posterior predictive check looks weird - what can I do?

**URL:** <https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648>\
**Category:** Modeling\
**Tags:** posterior-predictive, brms\
**Created:** [April 1, 2021, 7:33pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648 "2021-04-01T19:33:30Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 1, 2021, 7:33pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/1 "2021-04-01T19:33:31Z")

</div>

Hello together,

I am struggling with a posterior predictive check in brms that looks weird (no matter what I try). But first about my data: I conducted an experiment investigating jumping to conclusions bias. This means collecting few information before deciding. In different trials individuals could draw between 1 and 10 information. After each draw they could decide if they wanted another draw (information) or decide for option A or B. Less draws before deciding means higher jumping to conclusions bias…

Every individual did 24 trials in different scenarios. Furthermore, there were 3 groups.

The distribution of draws to decision (dependent variable) is bimodal with one peak in the middle and a lot of 10´s at the right end.

First I tried a gaussian family. Classic approach. That yielded the following posterior predictive check:

formula:

brm\_DTD ← brms::brm(  
formula = DTD ~ 1 + group\*scen + (1 | chiffre),  
family = gaussian(link = “identity”),  
data = data\_long[!is.na(data\_long$DTD),],  
control = list(adapt\_delta = 0.80, max\_treedepth = 10),  
warmup = 2000, iter = 10000, chains = 4,  
prior = prior, stanvars = stanvars, sample\_prior = TRUE,  
seed = 123, cores = parallel::detectCores())

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/cd70270ad5f2460ed8bdacc63c7d7397bfc82cbd.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/6a9c67b734db09c4e03e0409a1beac91ed9f64bb.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0e2bc85e2db11b8530660b997959e67e99c4b6f8.jpeg)

Then I tried a lot of things:

- student and skew\_normal looked pretty much the same.

- I know I have count data, so I tried poisson family predictive tons of value above 10 (even values of 15), so the gaussian fit was better.

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/4/4f0af62bfc8db09381cfd9cb5a16249c5c90017c.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/ba92cdf74ec0dd31fb08660de0c78f306e195020.png)

- Well, my count data is bounded at 10 - so I tried poisson with truncation between 1 and 10 (and 0 and 9 after subtracting 1 from my dependent variable). That looked pretty good concerning the dens\_overlay predictive check, but the stat\_2d check was not that good. As I am interested in the means I feel not safe interpreting this results (or what do you think?)! Furthermore there were some divergent iterations and I receive a warning message: “11% of all predicted values were invalid. Increasing argument ‘ntrys’ may help.”

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/c5acdeec6754b8a493938728b90cb1327b0b20f6.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/2c9dbc4a9b40212e0893c0f60537d718dc3ffdd0.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/2/20527bedad28fd14e1a90eec68bfef2f544472c9.png)

- I tried some other models like binomial family and some hurdle models (after transforming my data with abs(dependent\_variable -10) so that I had many 0´s instead of 10´s). The dens\_overlay check was between horrible and okay, but the stat\_2d check was never promising…

- I tried a gaussian mixture, but there were thousand of divergent transition and the model took some hours, so that wasn’t very satisfying as well…

What else can I try? I am grateful for every suggestion! Thanks a lot… :-)

Another question would be: If I will not find a solution to improve my posterior predictive check, can I still interpret the model? At least it is not worse than classic frequentist testing with lme4, right?

If you need more information please tell me…

Yours,

Simon

---

<div class="post-metadata">

**Author:** ![JLC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jlc/32/12520_2.png) [@JLC](https://discourse.mc-stan.org/u/JLC)\
**Post date:** [April 1, 2021, 11:11pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/2 "2021-04-01T23:11:31Z")

</div>

How do things look with a negative binomial distribution or a zero inflated negative binomial (if you have some zeros)?

Any chance there’s an additional predictor?

---

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 2, 2021, 7:24am UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/3 "2021-04-02T07:24:25Z")

</div>

Hi JLC,

thank you for the suggestion. I don’t think there is another predictor accounting for the bimodality.

I have tried a negbinomial (with dependent variable -1, so I have some zeros. Naturally it goes from 1 to 10…):

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/30629fca3e9bc4fd1c44c349921ae04218e38de3.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/d1743422dc716b70a04ee8a44665430bbccdcfb2.png)

Same with zero\_inflated negbinomial:

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/09b9f4831b567799b939faffced42f72b90c7bce.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/e/e568fd49c527739284f453f3609b39a9673bd9d6.png)

I repeated this with abs(dependent\_varable-10), so 1 becomes 9 and 10 becomes 0. I hoped his would improve results, as I have far more 10´s than 1`s. Well, not bad, but not much better than the gaussian and much harder to interpret!

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/6bac7bc1669381b40a7c470fb9db97d90cf083f9.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/3259ba119856f7671234fa956453bb1c42dcf738.png)

I also tried a zero inflated poisson on abs(dependent\_varable-10)… same thing: Not satisfying and very hard to interpret…

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/f/f37ada0bdd7c8ec0880c7b0a1d029c12de8477b3.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/00dcd5a969fd2b110b1cdf84acc4db5e09c6ddbd.png)

What do you think? Any other distribution that makes sense? And in doubt, can I interpret my gaussian model?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 2, 2021, 12:45pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/4 "2021-04-02T12:45:13Z")

</div>

Great that you are posting those detailed ploty. If your data are integers and bounded on both sides, then I think ordinal models would be a good choice. If I understand your experiment well, the continuation ratio (`cratio` in brms) family might be particularly suitable. Also check out the ordinal model tutorial by Paul Burkner for more details (currently on phone, so don’t have the link ready, sorry)

Best of luck with the model!

---

<div class="post-metadata">

**Author:** ![JLC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jlc/32/12520_2.png) [@JLC](https://discourse.mc-stan.org/u/JLC)\
**Post date:** [April 2, 2021, 1:02pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/5 "2021-04-02T13:02:43Z")

</div>

Could you post some pp\_checks with `type = "bars"` to `type = "bars_grouped`?

Bars can be helpful for investigating discrete outcomes.

---

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 2, 2021, 2:19pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/6 "2021-04-02T14:19:03Z")

</div>

Of course,

here is my groped bars pp\_check for the truncated poisson model (which is the best but still unsatisfying and harder to interpret alternative to the gaussian so far):

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/3c95ba30163d8b168a1d44ecc4e2409eacf24daf.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/5/518af99b90385c63361e91b0539735f4c014585c.png)

Its not that bad. The shape of the distribution is captured by the truncated poisson quite good, but the mean value and standard deviation is better reproduced by the gaussian model (which fails with respect to the dens\_overlay check…). I am not 100% happy with both so far…

With regard to the faceting I sent you am PM. And thank you @martinmodrak , I will try the cratio as soon as I can!

---

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 2, 2021, 2:26pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/7 "2021-04-02T14:26:14Z")

</div>

Ok I just naively changed the brmsfamily argument to “cratio”. The posterior predictive checks really look too good to be true!

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/aab291d1b5419963f5d19809c1ce1589c6671174.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/c/cdbc4c349a6be13fec34db0bac9bbcb0b8e73cc0.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/d73544e3c2667a4a001ea9ce15ffafc6bdb0074b.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/aaa8c0e5245383c0ca0689748c4a77ff109bb33e.png)

![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/9/929042e72be64ae840082fd77daf04b86a1bea3d.png)

I have to admit that I have no clue how to interpret the results. I have to read the vignette you mentioned and figure out if I can use the cratio model to answer my hypotheses as good as with the gaussian model but if I can, it is very promising!

Apart from that: Do you think with the given pp\_checks is my gaussian standard model a valid choice at all? On the one hand the dens\_overlay is really not fine. On the other hand its basically not worse than a typical lme4 frequentist model 99% of the scientists would choose. I am new to bayesian statistics and have no idea how serious this is!

---

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 2, 2021, 3:53pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/8 "2021-04-02T15:53:11Z")

</div>

@martinmodrak, do you mean this vignette?

[https://journals.sagepub.com/doi/10.1177/2515245918823199](https://journals.sagepub.com/doi/10.1177/2515245918823199)

I can’t find a description of the cratio family there. Generally, I do not find information about how the cratio family works and how to interprete the results. Maybe you do have a source or a link for me? Thank you very much for your support!

Simon

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 2, 2021, 4:30pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/9 "2021-04-02T16:30:32Z")

</div>

Yes, that’s the one. It seems it is called “Sequential models” there (I forgot about the naming conflict)

> [@Simon\_Wilms](#):
>
> Apart from that: Do you think with the given pp\_checks is my gaussian standard model a valid choice at all? On the one hand the dens\_overlay is really not fine. On the other hand its basically not worse than a typical lme4 frequentist model 99% of the scientists would choose. I am new to bayesian statistics and have no idea how serious this is!

Tough question. I’ll go with Box: “Since all models are wrong the scientist must be alert to what is importantly wrong. It is inappropriate to be concerned about mice when there are tigers abroad.” In other words, depends on you problem and the inferences you need to make. Since it appears the Gaussian matches your mean and standard deviation OK, then if it also does this well in subgroups you care about and you only make inferences about mean, it might not really matter that much.

With that said I personally think the `cratio` is overall a safer choice and IMHO very closely matches your experiment and is thus highly interpretable (but make sure it matches, I might have misunderstood things).

> [@Simon\_Wilms](#):
>
> The posterior predictive checks really look too good to be true!

That’s because the ordinal families are very flexible and generally have N-1 parameters for N outcomes - so it is not suprpring they fit the overall distribution. The fits in subgroups are of more interest.

---

<div class="post-metadata">

**Author:** ![JLC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jlc/32/12520_2.png) [@JLC](https://discourse.mc-stan.org/u/JLC)\
**Post date:** [April 2, 2021, 5:04pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/10 "2021-04-02T17:04:43Z")

</div>

Would `loo_compare` provide useful comparative diagnostics between the models in this case?

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [April 2, 2021, 8:23pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/11 "2021-04-02T20:23:17Z")

</div>

> [@JLC](#):
>
> Would `loo_compare` provide useful comparative diagnostics between the models in this case?

Good question! First, there is a technical problem for this specific case: `loo` can’t directly compare fits with continuous families and discrete families. (see the relevant entry in Aki’s [Cross-validation FAQ](https://avehtari.github.io/modelselection/CV-FAQ.html#12_Can_cross-validation_be_used_to_compare_different_observation_models __response_distributions__ likelihoods) )

Second, there is the general consideration that `loo` while useful should (in my personal opinion) be of second-order interest and the most important input should be formed by the question being asked. `loo` will try to estimate which model is more likely to be better at predicting out-of-sample data. For many questions, we may not want to optimize prediction capacity, but rather understanding, interpretability or “correctness” in some sense. For example, in causal inference, a model that does not adjust for a confounder or that conditions on a collider (or both) might easily be preferred by `loo` despite not providing a good answer to our question. Similarly when choosing a response distribution, it might sometimes be sensible to choose say a more interpretable option over the one favored by `loo`.

---

<div class="post-metadata">

**Author:** ![JLC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jlc/32/12520_2.png) [@JLC](https://discourse.mc-stan.org/u/JLC)\
**Post date:** [April 2, 2021, 9:00pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/12 "2021-04-02T21:00:57Z")

</div>

Thank you very much for the great explanation!

---

<div class="post-metadata">

**Author:** ![Simon\_Wilms](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/simon_wilms/32/11162_2.png) [@Simon\_Wilms](https://discourse.mc-stan.org/u/Simon_Wilms)\
**Post date:** [April 2, 2021, 9:09pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/13 "2021-04-02T21:09:28Z")

</div>

Thank you, I like your point of view! Tomorrow I will read the vignette you recommended… I am excited if it solves my problem!

One question that comes to my mind: If I simply have more parameters for the cratio model and as a result my pp\_check always looks fantastic - am I not overfitting my model? However, maybe this will become clearer after reading the text tomorrow… I will come back here after that!

For today, big thanks for your support @JLC and @martinmodrak. Your answers helped me a lot!

---

<div class="post-metadata">

**Author:** ![JLC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jlc/32/12520_2.png) [@JLC](https://discourse.mc-stan.org/u/JLC)\
**Post date:** [April 2, 2021, 10:18pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/14 "2021-04-02T22:18:34Z")

</div>

In case it’s at all helpful, here is how I have munged some ppc’s for discrete outcomes with more than 1 grouping factor:

```
# Get posterior predictions
preds <- data %>%
  add_predicted_draws(model,
                      re_formula = NULL,
                      prediction = "Pred",
                      n = 100)

# Calculate median counts per outcome and CI's
pred_counts <- preds %>%
  ungroup() %>%
  select(-.row, -.chain, -.iteration) %>%
  group_by(Cov1,
           Cov2,
           Prediction, 
           .draw) %>%
  summarize(Count = n()) %>%
  ungroup() %>%
  group_by(Cov1,
           Cov2,
           Pred) %>%
  point_interval(
    Count,
    .width = 0.89,
    .point = median,
    .interval = hdci,
    .simple_names = TRUE,
    na.rm = TRUE
  ) 

# Plot it out
ggplot() +
  geom_bar(data = data,
           aes(x = Y, fill = "Observed"), colour = "black") +
  facet_grid(Cov1 ~ Cov2) +
  geom_point(data = pred_counts,
             aes(x = Pred, y = Count, shape = "Predicted"),
             size = 2) +
  geom_errorbar(data = pred_counts,
                aes(x = Pred,
                    ymin = .lower,
                    ymax = .upper),
                width = 0.2)

```

---

<div class="post-metadata">

**Author:** ![jerlich](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jerlich/32/2131_2.png) [@jerlich](https://discourse.mc-stan.org/u/jerlich)\
**Post date:** [April 6, 2021, 9:25am UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/15 "2021-04-06T09:25:21Z")

</div>

To avoid overfitting, do cross-validation. You can use the `kfold` function. But depending on the size of your data, you may want to generate the folds yourself (e.g. 10% from each subject for each fold).

---

<div class="post-metadata">

**Author:** ![nguyenthiphuong](https://avatars.discourse-cdn.com/v4/letter/n/59ef9b/32.png) [@nguyenthiphuong](https://discourse.mc-stan.org/u/nguyenthiphuong)\
**Post date:** [February 2, 2024, 1:38am UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/16 "2024-02-02T01:38:23Z")

</div>

Hi Simon @Simon_Wilms, Can you share your codes to make the graph for pp\_check() with stat = “mean” that can show values on y-axis! I used the default function and cannot show y-axis values and scales! Thanks

---

<div class="post-metadata">

**Author:** ![aaditya44](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/aaditya44/32/18367_2.png) [@aaditya44](https://discourse.mc-stan.org/u/aaditya44)\
**Post date:** [April 24, 2024, 1:33pm UTC](https://discourse.mc-stan.org/t/posterior-predictive-check-looks-weird-what-can-i-do/21648/17 "2024-04-24T13:33:40Z")

</div>

The default theme used within the pp\_check() function masks the y-axis for some reason. just add + theme\_classic() for example, and you have the y-axis again.