# Discrepancy in Manually- vs Automatically-Generated pp\_check() Figures

**URL:** <https://discourse.mc-stan.org/t/discrepancy-in-manually-vs-automatically-generated-pp-check-figures/17545>\
**Category:** brms\
**Tags:** techniques, specification\
**Created:** [August 19, 2020, 1:41pm UTC](https://discourse.mc-stan.org/t/discrepancy-in-manually-vs-automatically-generated-pp-check-figures/17545 "2020-08-19T13:41:01Z")\
**Posts on this page:** 3\
**Page:** 1

<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:** [August 19, 2020, 1:41pm UTC](https://discourse.mc-stan.org/t/discrepancy-in-manually-vs-automatically-generated-pp-check-figures/17545/1 "2020-08-19T13:41:01Z")

</div>

I am working with a multivariate ordered probit model fit with `brms` with the following form:

```
Modality_ord_mv_complete_fmla <- bf(mvbind(Valence, Arousal) | resp_thres(9, gr = Category) ~ Category +
                                   Condition +
                                   Category:Condition +
                                   Intensity.c +
                                   (1|o|Token) +
                                   (1|p|Condition) +
                                   (1|q|Condition:Subject)) +
                lf(disc ~ 0 + Category, resp = "Valence", cmc = FALSE) +
                lf(disc ~ 0 + Category, resp = "Arousal", cmc = FALSE)

```

Using the built in pp\_check functions, I get reasonable looking posterior predictions:

```
pp_check(Modality_Ord_Complete_MV_Mod_Const_Thresh, type = "bars_grouped", group = "Token", resp = "Valence", nsamples = 10)

```

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

But, when I try to manually created this plot, or one similar I get very different results, so I’m wondering where I’m going wrong:

```
# Create a data frame of original data and add predictions
Modality_Preds <- cbind(Modality_Ord_Full %>%
        select(Subject, Condition, Category, Token, Pitch.c, Intensity.c, Valence) %>%
        na.omit(),
        as.data.frame(apply(posterior_predict(Modality_Ord_Complete_MV_Mod_Const_Thresh, newdata = Modality_Ord_Full %>%
                                                      select(Subject, Condition, Category, Token, Pitch.c, Intensity.c, Valence) %>%
                                                      na.omit(),
                                              re_formula = NULL,
                                              nsamples = 10),
                            c(2,3),
                            mean)) %>% # Take the mean across the sample dimension
                rename("Pred_Valence" = Valence,
                       "Pred_Arousal" = Arousal) %>%
                mutate(Pred_Valence = as.integer(floor(Pred_Valence)),
                       Pred_Arousal = as.integer(floor(Pred_Arousal)))) %>%
        select(-Pred_Arousal) %>%
        mutate(Pred_Valence = factor(Pred_Valence, levels = c(1:9), ordered = TRUE))

## Plot combined original data with predicted responses
Modality_Preds %>%
        select(Condition, Category, Token, Valence, Pred_Valence) %>% 
        ggplot() +
        facet_wrap( ~ Token, ncol = 8) +
        geom_bar(aes(x = Valence)) +
        geom_point(aes(x = Pred_Valence), stat = "count") +
        theme_bw() +
        scale_x_discrete("Rating (Valence)") +
        scale_y_continuous("Count") +
        theme( axis.text.x = element_text(angle = 0, hjust = 0.5),
                plot.title = element_text(hjust = 0.5),
                panel.border = element_rect(fill = NA, colour = "black"),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                axis.line = element_line(colour = "black"),
                strip.text.x = element_text(size = 12),
                strip.background = element_blank(),
                legend.position = c(1, 0), legend.justification = c(1, 0),
                legend.background = element_rect(fill = "transparent"),
                legend.title = element_blank()
        )

```

 ![Token_Rating_Count_Pred](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/5/53857e56f41d048dd5846950962c14a1c4b6efeb.png)

Could anyone help me understand how to reproduce the pp\_check() plots and/or where I’m going wrong?

---

<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:** [August 19, 2020, 3:15pm UTC](https://discourse.mc-stan.org/t/discrepancy-in-manually-vs-automatically-generated-pp-check-figures/17545/2 "2020-08-19T15:15:55Z")

</div>

I think I’ve found at least part of the difference, if I manually work through creating counts with the expected predicted probabilities, I get much closer to the automatic plot:

```
med_epred <- as.data.frame(apply(posterior_epred(Modality_Ord_Complete_MV_Mod_Const_Thresh, newdata = Modality_Ord_Full %>%
                                                         select(Subject, Condition, Category, Token, Pitch.c, Intensity.c, Valence) %>%
                                                         na.omit(),
                                                 re_formula = NULL,
                                                 resp = "Valence",
                                                nsamples = 10),
                              c(2,3),
                              median)) %>%
        mutate(Subject = Modality_Ord_Full %>%
                       pull(Subject),
               Token = Modality_Ord_Full %>%
                       pull(Token)) %>%
        pivot_longer(cols = 1:9,names_to = "Rating", values_to = "Pred_Prob")

Rating_Counts <- Modality_Ord_Full %>%
        #filter(Condition == "AO") %>%
        group_by(Token) %>%
        summarize(Count = n())
        

Epred_Counts <- left_join(med_epred,
                          Rating_Counts,
                          by = "Token") %>%
        mutate(Pred_Count = Pred_Prob*Count) %>%
        group_by(Token, Rating) %>%
        summarize(Pred_Count = median(Pred_Count))

```

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

---

<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:** [August 25, 2020, 1:18pm UTC](https://discourse.mc-stan.org/t/discrepancy-in-manually-vs-automatically-generated-pp-check-figures/17545/3 "2020-08-25T13:18:00Z")

</div>

Sorry that we didn’t react sooner - your new plot looks to me identical to the “official” pp\_check, so if I understand it correctly you did resolve the issue. Or is there something you still need help with?