Search code examples
rggplot2intervalsfacetdensity-plot

Draw interval on geom_density


How do I draw a horizontal line indicating the Highest (Posterior) Density interval for faceted density plots in ggplot2? This is what I have tried:

# Functions to calculate lower and upper part of HPD.
hpd_lower = function(x) coda::HPDinterval(as.mcmc(x))[1]
hpd_upper = function(x) coda::HPDinterval(as.mcmc(x))[2]

# Data: two groups with different means
df = data.frame(value=c(rnorm(500), rnorm(500, mean=5)), group=rep(c('A', 'B'), each=500))

# Plot it
ggplot(df, aes(x=value)) + 
  geom_density() +
  facet_wrap(~group) + 
  geom_segment(aes(x=hpd_lower(value), xend=hpd_upper(value), y=0, yend=0), size=3)

enter image description here

As you can see, geom_segment computes on all data for both facets whereas I would like it to respect the faceting. I would also like a solution where HPDinterval is only run once per facet.


Solution

  • Pre-calculate the hpd intervals. ggplot evaluates the calculations in the aes() function in the entire data frame, even when data are grouped.

    # Plot it
    library(dplyr)
    df_hpd <- group_by(df, group) %>% summarize(x=hpd_lower(value), xend=hpd_upper(value))
    
    ggplot(df, aes(x=value)) + 
      geom_density() +
      facet_wrap(~group) + 
      geom_segment(data = df_hpd, aes(x=x, xend=xend, y=0, yend=0), size=3)
    

    enter image description here