How to best summarize this evidence

[Eric_Green]

I have a conceptual question about how to present and summarize results. I’d be happy to delete if this is not the best forum.

I fit 3 models in brms and used a combination of tidybayes, bayesplot, and bayestestR functions to create and annotate my plot.

As a budding Bayesian, I’d like to know how others might summarize this evidence. Here’s what I wrote:

Figure 1 displays the results of three separate difference-in-differences regressions that estimate the effect of disclosure on child, caregiver, and family well-being. Each plot shows Markov chain Monte Carlo draws from the posterior distribution of the coefficient of the interaction between time (12-month follow-up = 1) and disclosure status (disclosed = 1). These coefficients represent the standardized effect of disclosing on the well-being outcomes.

Overall, the evidence is inconclusive. It appears in Model 3, where the dependent variable is family relationship quality, that relationship quality may have decreased in a meaningful way following disclosure. The median of the posterior of the coefficient is -0.35 standard deviations, and 83% of the posterior distribution is consistent with a negative effect (i.e., less than -0.1 standard deviations, the lower bound to our defined region of practical equivalence, ROPE). However, this is only weak evidence, as a non-trivial proportion of the posterior distribution falls within or to the right of the ROPE.

Model 1 suggests the possibility of a slight decrease in child problems following disclosure, but the uncertainty around the estimate is too large to state this with confidence. Likewise, Model 2 suggests that there was no practical effect of disclosure on caregiver depression severity, but we cannot rule out a modest increase or decrease.

[bbbales2]

This has been sitting here a while. I don’t know the specifics of the model, or application, or goals, but for lack of better answer I’ll make some suggestions :D.

Since you have a pretty graph to offer up, you can offer it up to the blog, though if it gets on the blog roll it might take a while to get around to it and you may or may not get good feedback.

A couple generic comments though – what do the dashed verticle lines (the ROPE things tell me)?

Also the probabilities printed on the plot aren’t the things shaded in yellow. Why is the yellow shaded? Why is the median of the distribution plotted (seems easy enough to eyeball)? Those probabilities too – is 57.3% meaningful vs. 57% or \approx 57% or something?

These are interaction effects too – could you display more information about the model (the non-interaction effects, or other interations) by going to different summaries? Somehow these things look pretty well behaved – no big skewness and whatnot. Presumably there are other parameters in the model and maybe instead of a super-deep dive on these a reader would be better served knowing them all? You can get much more mileage with dot/line plots, see here

[avehtari]

I think your plots are good and your summarizing interpretations are valid.

[Eric_Green]

Thanks for sharing some ideas, @bbbales2. The region between the ROPE lines represents values that are small enough to be practically null. (Or maybe your question is meant to suggest I should clarify in the plot?).

The yellow region covers the 89% credible interval. The text brings attention to quantities of interest in each plot. In Model 2, the posterior is centered on zero, but only 31% falls within the ROPE, so IMO there is not a strong conclusion to draw about the size or direction of the effect. The other two model results fall outside the ROPE so the text calls attention to the directional questions we want to consider.

[Eric_Green]

I appreciate your taking the time to weigh in, @avehtari. Thanks.

ps. Really enjoying Regression and Other Stories right now.

[JimBob]

Those are really nice and very informative graphs - I’m surprised you only call yourself a ‘budding’ Bayesian. If you were willing, I would really be interested to see the code (if you don’t wish to share the data itself, even just seeing the plot codes would be really interesting).

The only two things I think might be useful would be to make sure, if not in the graph or figure caption itself, at least in the paper, what exactly the ROPE means like bbbales noted.

The second thing which I think might go in the plot itself is to make sure you really mean ‘credible interval’. Is it a Highest Density Interval (the 89% most likely values in the posterior) or is it an Equal Tailed Interval (the 89% middle values). People who say Credible Interval are often referring to an ETI. However, it is a bit weird to say that simply because something is a middle value that it is a credible value. Long story short, I would always says directly what the interval is, rather than giving it an interpretive name, even though this is what a lot of people do.

[Eric_Green]

Thanks, @JimBob. I should be able to post all before too long.

I was trying to get clarity on credible interval vs HDI, so I’m glad you raised it. The core plotting is done with bayesplot::mcmc_areas. The width is controlled by the prob parameter:

The probability mass to include in the inner interval (for mcmc_intervals() ) or in the shaded region (for mcmc_areas() ). The default is 0.5 (50% interval) and 1 for mcmc_areas_ridges() .

Someone commenting on this SO question said that mcmc_areas uses ETI. Once I confirm that’s the case I’ll update the plot to clarify.

As for removing the “credible interval” label, I wonder if you can point me to more on this topic. This explanation from the bayestestR docs frames HDI and ETI as two different methods for computing a credible interval.

[JimBob]

I think that explanation from the bayestestr package is a pretty good summary. I don’t have much more information on it. Personally, I prefer the HDI because then you can really say ‘the 89% most likely values were between bla bla’ - which is tempting to do with the ETI but not quite accurate. Usually the ETI and HDI will correspond or approximately correspond, but I’d rather just be able to say what exactly I mean. One thing that is sometimes considered a ‘benefit’ of the ETI is that if you transform the estimates, it would remain on the same numbers (though transformed), whereas with the HDI the transformed values can lead to a different set of estimates and perhaps affect assessments when using a ROPE, for example.

I have a bit of a sense that some people use the term credible interval because you can shorten it to CI and then certain reviewers won’t want extra explanation because it seems more like a ‘confidence interval’ (I’ve even seen this humorously and explicitly mentioned as a benefit of saying CI in a talk). It’s just my personal preference to simply say either HDI or ETI directly to avoid confusion, and also I don’t see why there is anything strictly credible about the 89% middle values in an ETI if there could be other values that are more likely. Certainly not everyone feels the same way though!

[Eric_Green]

Thanks for the additional info, @JimBob.

Adding this from the {bayesplot} docs for my future self and for anyone who lands on this thread:

For models fit using MCMC we can compute posterior uncertainty intervals (sometimes called “credible intervals”) in various ways. bayesplot currently provides plots of central intervals based on quantiles, although additional options may be provided in future releases (e.g., HDIs, which can be useful in particular cases).

[JimBob]

Great - and please do update when you might be willing to share your plot code. The annotations and titling are very sleek (are the plots actually patched together in R, or did you make separate plots and then combine them in photoshop or something like that?).

If you do wish to get a HDI, then I believe there is a HDI function in the tidybayes package, and also the R package HDInterval is specifically for that purpose :

HDInterval: Highest (Posterior) Density Intervals

[Eric_Green]

Thanks for the additional resource, @JimBob.

Here’s a version of the code applied to a generic brms example. I have no idea if the second fit makes sense—I just wanted to create a second model to show how to piece everything together. It’s a 100% R plot. No post-processing.

library(tidyverse)
library(brms)
library(bayesplot)
library(tidybayes)
library(patchwork)
library(ggtext)

fit1 <- brm(count ~ zAge + zBase * Trt + (1|patient), 
            data = epilepsy, family = poisson())

fit2 <- brm(count ~ zAge + zBase * Trt + (1|patient), 
            data = epilepsy, family = zero_inflated_poisson())

fit1.p <- as_tibble(as.matrix(fit1)) %>%
  rename(`Treatment` = `b_Trt1`) %>%
  select(`Treatment`)

fit2.p <- as_tibble(as.matrix(fit2)) %>%
  rename(`Treatment` = `b_Trt1`) %>%
  select(`Treatment`)

fit1H <- as.numeric((hypothesis(fit1, 
                                    "Trt1 > 0.1")$hypothesis["Post.Prob"])*100)

fit2H <- as.numeric((hypothesis(fit2, 
                                "Trt1 > 0.1")$hypothesis["Post.Prob"])*100)

fit1.median <- median_hdi(fit1.p) %>%
  pull(`Treatment`)

fit2.median <- median_hdi(fit2.p) %>%
  pull(`Treatment`)


color_scheme_set(c("#FDE725", "#440154", "#440154",
                   "#440154", "#440154", "#440154"))
bayesplot_theme_set(theme_bw())


p1 <- fit1.p %>%
  mcmc_areas(prob = 0.89, point_est = "median") +
  xlab("Title") +
  theme(#plot.title.position = "plot",
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    plot.title = element_markdown(face="bold", size=16)) +
  xlim(-1,1) +
  geom_vline(xintercept = 0, linetype="solid") +
  geom_vline(xintercept = -0.1, linetype="dashed") +
  geom_vline(xintercept = 0.1, linetype="dashed") +
  labs(title="Model 1: Some title",
       subtitle = "Something here") +
  theme(axis.text = element_text(size=14),
        axis.title = element_text(size=15),
        plot.subtitle = element_text(size=14)) +
  annotate("text", x = 0.3, y = 1.8, hjust = 0,
           label = "Some label",
           size = 4, color="#818283") + 
  geom_curve(aes(x = 0, y = 1.8, xend = 0.25, yend = 1.8), 
             color = "#818283", 
             arrow = arrow(type = "open",
                           length = unit(0.1, "inches")), 
             curvature = 0, angle = 100, ncp =15) +
  annotate("text", x = -Inf, y = 1.7, hjust = -0.08,
           label = paste0("P(effect > 0.1): ",
                          round(fit1H, 1), "%"),
           size = 4.3, color="black") + 
  annotate("text", x = -0.7, y = 1.1, hjust = 0,
           label = paste0("Median: ",
                          round(fit1.median, 2)),
           size = 4.3, color="#440154") 

p2 <- fit2.p %>%
  mcmc_areas(prob = 0.89, point_est = "median") +
  xlab("Title") +
  theme(#plot.title.position = "plot",
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    plot.title = element_markdown(face="bold", size=16)) +
  xlim(-1,1) +
  geom_vline(xintercept = 0, linetype="solid") +
  geom_vline(xintercept = -0.1, linetype="dashed") +
  geom_vline(xintercept = 0.1, linetype="dashed") +
  labs(title="Model 2: Some title",
       subtitle = "Something here") +
  theme(axis.text = element_text(size=14),
        axis.title = element_text(size=15),
        plot.subtitle = element_text(size=14)) +
  annotate("text", x = 0.3, y = 1.7, hjust = 0,
           label = "Some label",
           size = 4, color="#818283") + 
  geom_curve(aes(x = 0, y = 1.8, xend = 0.25, yend = 1.8), 
             color = "#818283", 
             arrow = arrow(type = "open",
                           length = unit(0.1, "inches")), 
             curvature = 0, angle = 100, ncp =15) +
  annotate("text", x = -Inf, y = 1.8, hjust = -0.08,
           label = paste0("P(effect > 0.1): ",
                          round(fit2H, 1), "%"),
           size = 4.3, color="black") + 
  annotate("text", x = -0.7, y = 1.1, hjust = 0,
           label = paste0("Median: ",
                          round(fit2.median, 2)),
           size = 4.3, color="#440154") 

((p1 / p2)) + plot_annotation(
  title = 'Overall title here',
  subtitle = 'Overall subtitle here',
  caption = paste0("Data from ", nrow(distinct(epilepsy)), " participants, etc."),
  theme = theme(plot.title = element_text(face="bold", size=16))
) & theme(plot.caption = element_markdown(size=12, color="#686869")) 

[JimBob]

Thanks for sharing - that’s awesome!

[jerlich]

These plots are nice! and pretty clear.

But I’d like to 2nd some of the comments from @bbbales2.

You are plotting the interaction effect of followup and disclosure. But we don’t know the main effect of time or disclosure, so i’m not sure these plots tell us enough about the effects of disclosure. Like, if there was a massive effect of disclosure that didn’t change across time, then this interaction is going to be small, but your plot title de-emphasizes the time-interaction element of the plots.

I also think you make the readers work harder by having some of the axes have a more is worse direction and some have a more is better direction.

why is yellow 89%? not 90% or 95%? is there something about 89%

I’m also not a huge fan of the bright yellow, but that’s personal taste :)

[Eric_Green]

Thanks for the feedback, @jerlich. In this diff-in-diff model we’re most interested in the interaction, but we’ll include a table of full model results in the appendix.

I think you’re right about the direction. We need to decide on the tradeoff between presenting the results on the original scale or reversing one so they all go in the same direction. My sense is that our audience might be more offended by mucking with the scale of the original measure, but I could be wrong!

I’m following McElreath’s convention to be unconventional: https://twitter.com/rlmcelreath/status/846674692815704066?s=20

The {bayestestR} package does the same: Credible Intervals (CI) • bayestestR