Attaching labels estimating proportion of posterior density within a fixed interval in the ggdist package

General
2

I’m not sure how to do this within ggdist (or if it can be done), but @mjskay probably knows. For now, here’s a way to construct the plot by directly calculating the desired density values from the posterior draws and then creating a custom ggplot.

The overlapping shading of the ROPE and HDI might be confusing, but hopefully you can play around with the code to get something that meets your needs. (I’ve also added this answer to the Stackoverflow question.)

# Calculate HDIs
hdi = compCat %>% 
  group_by(Species, SepalLength_cat) %>% 
  summarise(as_tibble(hdi(.epred))) %>%
  rename(hdi.lwr=V1, hdi.upr=V2)

# Calculate densities
adj = 0.5
dens.hist = compCat %>% 
  group_by(Species, SepalLength_cat) %>% 
  group_split() %>% 
  map_df(~{
    dens = density_bounded(.x$.epred, adjust=adj)
    tibble(Species=.x$Species[1], SepalLength_cat=.x$SepalLength_cat[1],
           x=dens$x, .epred=dens$y)
  }) %>% 
  left_join(hdi)

# Calculate ROPE density
rope.rng = c(-0.25, 0.25)

rope.dens = dens.hist %>% 
  filter(x >= min(rope.rng), x <= max(rope.rng)) %>% 
  group_by(Species, SepalLength_cat) %>% 
  summarise(rope.dens = sum(mean(diff(x)) * .epred))

rope.color = hcl(0, 80, 30, alpha=c(1, 0.7))
hdi.color = hcl(240, 100, 30, alpha=c(1,0.5))
lwd = 0.3
theme_set(theme_minimal(base_size=14))

dens.hist %>% 
  ggplot(aes(x, .epred)) +
    geom_hline(yintercept=0, linewidth=0.1, colour="black") +
    # Overall densities
    geom_line(linewidth=lwd) +
    # Shade ROPE
    geom_area(data=. %>% 
                filter(x >= min(rope.rng), x <= max(rope.rng)),
              aes(fill="ROPE")) +
    # Shade HDI
    geom_area(data=. %>% 
                filter(x >= hdi.lwr, x <= hdi.upr),
              aes(fill="HDI")) +
    # HDI segment
    geom_segment(data=hdi, 
                 aes(x=hdi.lwr, xend=hdi.upr, y=-0.1, yend=-0.1,
                     colour="HDI")) +
    # ROPE segment
    geom_line(data=. %>% 
                group_by(Species, SepalLength_cat) %>% 
                reframe(x=rope.rng, .epred=-0.2),
              aes(colour="ROPE")) +
    # ROPE density value
    geom_text(data=rope.dens, 
              size=3.5, vjust=0.5,
              aes(label=paste0(round(rope.dens*100, 1),"%"), 
                  x = mean(rope.rng), 
                  y = -0.4,
                  colour="ROPE")) +
    facet_grid(cols=vars(Species), row=vars(fct_rev(SepalLength_cat)), switch="y") +
    theme(
      strip.text.y.left=element_text(angle=0, vjust=0.02),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank()
    ) +
    scale_y_continuous(expand=c(0.03, 0.03)) +
    scale_fill_manual(values=c(ROPE=rope.color[2], HDI=hdi.color[2])) +
    scale_colour_manual(values=c(ROPE=rope.color[1], HDI=hdi.color[1])) +
    guides(colour="none") +
    labs(y="Sepal Length comparison", x="Sepal Width", fill=NULL)

Rplot05
Rplot05919×611 51.3 KB

2 Likes