How to calculate and label peak value of distribution by multiple conditions/facets in R ggplot?

Question

While the question appears similar to others, there's a key difference in my mind.

1. I want to be able to calculate and/or print (graphing it would be the ultimate goal, but calculating it in the data frame the primary goal) the peak value of a density curve of EACH SUB-CONDITION BY FACET The density graph looks like this:

So, ideally, I would be able to know the intensity (x-axis value) corresponding to the highest peak of the density curves for each condition.

Here's some dummy data:

set.seed(1234)

library(tidyverse)
library(fs)
n = 100000
silence = factor(c("sil1", "sil2", "sil3", "sil4", "sil5"))
treat = factor(c("con", "uos", "uos+wnt5a", "wnt5a"))
silence = rep(silence, n)
treat = rep(treat, n)
intensity = sample(4000:10000, n)

df <- cbind(silence, treat, intensity)
df$silence <- silence
df$treat <- treat

2. What I've tried:

• Subsetting the primary DF and going through and calculating the density of each condition, but this could take days
• Something close to this answer: Calculating peaks in histograms or density functions but not quite. I think the data look better as a histogram personally, but that constructs an arbitrary number of bins for intensity data (a continuous measure). The histogram looks like this:

Again, it would be sufficient to get the peak values for each of these groups (i.e., treatments by silencing subdistributions) just in the console, but adding them as a vertical line in the graphs would be a sweet cherry on top (it could also make it hella busy, so I will see about that piece later)

Thank you!!

Answer 1

Depending on the way you're producing the density plots, there may be a more direct way to recreate the density calculation before it goes into ggplot. That'll be the easiest way to get the peak values and keep them in the format of your data.

Without that, here's a hack that should work in general, but requires some kludging to fit the extracted points back into the form of your original data.

Here's a plot like yours:

mtcars %>% 
  mutate(gear = as.character(gear)) %>%
  ggplot(aes(wt, fill = gear, group = gear)) +
  geom_density(alpha = 0.2) +
  facet_wrap(~am) ->my_plot

Here are the components that make up that plot:

ggplot_build(my_plot) -> my_plot_innards

With some ugly hacking we can extract the points that make up the curves and make them look kind of like our original data. Some info is destroyed, e.g. the gear values 3/4/5 become group 1/2/3. There might be a cool way to convert back, but I don't know it yet.

extracted_points <- tibble(
  wt = my_plot_innards[["data"]][[1]][["x"]],
  y = my_plot_innards[["data"]][[1]][["y"]],
  gear = (my_plot_innards[["data"]][[1]][["group"]] + 2) %>% as.character, # HACK
  am = (my_plot_innards[["data"]][[1]][["PANEL"]] %>% as.numeric) - 1 # HACK
)

ggplot(extracted_points, aes(wt, y, fill = gear)) +
  geom_point(size = 0.3) +
  facet_wrap(~am)

extracted_points_notes <- extracted_points %>%
  group_by(gear, am) %>%
  slice_max(y)


my_plot +
  geom_point(data = extracted_points_notes,
             aes(y = y), color = "red", size = 3, show.legend = FALSE) +
  geom_text(data = extracted_points_notes, hjust = -0.5,
             aes(y = y, label = scales::comma(y)), color = "red", size = 3, show.legend = FALSE)