Dear @ASKurz, Reading chapter 20, section 20.2.3, I found out that the HDIs and estimates that you calculated do not match the ones in the book. In that section, Krushcke ran a model with Salary as DV and two IDs (Pos and Org). Then, he presented the main effect comparisons as the figure below shows. For example, the left panel shows that the mode is 13400 with 9720 to 16700 HDIs. Screen Shot 2021-08-26 at 2 02 48 pm

In the bookdown, to plot this figure, the `fitted` function is used: ```r # define the new data nd <- tibble(Pos = c("Assis", "Assoc"), Org = "mu") # feed the new data into `fitted()` f <- fitted(fit20.1, newdata = nd, summary = F, allow_new_levels = T) %>% as_tibble() %>% set_names("Assis", "Assoc") %>% mutate(`Assoc vs Assis` = Assoc - Assis) # plot f %>% ggplot(aes(x = `Assoc vs Assis`, y = 0)) + stat_beedrill() + scale_y_continuous(NULL, breaks = NULL) + labs(title = "Assoc vs Assis", x = "Difference") ``` And, this creates the following figure: Screen Shot 2021-08-26 at 1 58 33 pm

As you can see, the HDIs do not match with the ones in the book. Just to make sure, let's check them more carefully: ```r f %>% mode_hdi(`Assoc vs Assis`) %>% select(`Assoc vs Assis`:.upper) %>% mutate_if(is.double, round, digits = 0) ## # A tibble: 1 x 3 ## `Assoc vs Assis` .lower .upper ## ## 1 11260 -13545 40783 ``` I am not familiar with the `fitted` function, but it seems that HDIs of main effect comparisons are not correct. Interestingly, using this method for plotting simple effect and interaction comparisons produces correct HDIs. I guess the solution would be easy, but since I am not familiar at all with `fitted`, I used another method that you have used in previous chapters. This method gave the correct HDIs for main effect comparisons: ```r post_pos <- coef(fit20.1, summary = F)$Pos %>% as_tibble() %>% select(Assis = Assis.Intercept, Assoc = Assoc.Intercept) %>% mutate(`Assoc vs Assis` = Assoc - Assis) post_pos %>% ggplot(aes(x = `Assoc vs Assis`, y = 0)) + stat_beedrill() + scale_y_continuous(NULL, breaks = NULL) + labs(title = "Assoc vs Assis", x = "Difference") ``` Screen Shot 2021-08-26 at 1 59 15 pm

```r post_pos %>% mode_hdi(`Assoc vs Assis`) # `Assoc vs Assis` .lower .upper .width .point .interval # #1 13681. 8475. 18027. 0.95 mode hdi ``` However, the problem with this method is that while it gives the correct HDIs and estimates for main and interaction comparisons, it does not work well with simple effects: ```r post_pos_org <- coef(fit20.1, summary = F)$`Pos:Org` %>% as_tibble() %>% select(Assis_PSY.Intercept, Full_PSY.Intercept, Assis_CHEM.Intercept,Full_CHEM.Intercept) %>% transmute (`Full - Assis @ PSY` = Full_PSY.Intercept - Assis_PSY.Intercept, `Full - Assis @ CHEM` = Full_CHEM.Intercept - Assis_CHEM.Intercept) %>% mutate(`Full.v.Assis (x) CHEM.v.PSY` = `Full - Assis @ CHEM` - `Full - Assis @ PSY`) %>% gather() post_pos_org %>% group_by(key) %>% mode_hdi(value) # key value .lower .upper .width .point .interval # # 1 Full - Assis @ CHEM 25846. 10907. 40284. 0.95 mode hdi # 2 Full - Assis @ PSY -11968. -26198. 2658. 0.95 mode hdi # 3 Full.v.Assis (x) CHEM.v.PSY 38480. 17523. 57330. 0.95 mode hdi ``` If you compare the output to the following figure, you can see that while they match the values in the interaction plot (the right one), they are extremely inconsistent with the left and middle plots. Screen Shot 2021-08-26 at 2 19 54 pm

While it would be great to know the problem with the first method (i.e., `fitted`), I would be grateful if you could tell me what is wrong with the second method (i.e., `coef(fit20.1, summary = F)$Pos:Org`). Thanks a lot.

GitHub Explorer

Doing-Bayesian-Data-Analysis-in-brms-and-the-tidyverse

Score Breakdown

Issues Activity: Last 6 months

Top Labels

Hottest Issues