Extracting-peak-values-from-a-density-ggplot.R

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

plot <- data.frame(Values=c(1,1,1,2,2,3,3,5,4,5,5,5,5,6,6,6,6,5,5,5,
                            4,4,4,7,8,8,8,8,8,8,8,4,3,2,2,3,5,6,7,4,
                            2,2,3,9,9,9,9,9,9,9,9,9,9,5,1,3,3,1)) %>% 
  ggplot(aes(x=Values))+
  geom_density()+
  theme_classic() 
plot

p <- ggplot_build(plot)

plot +
  geom_vline(xintercept = p$data[[1]]$x[which.max(p$data[[1]]$y)]) +
  geom_vline(xintercept = p$data[[1]]$x[p$data[[1]]$x>7][which.max(p$data[[1]]$y[p$data[[1]]$x>7])]) 

######################################method_2
d <- density(mtcars$disp)
d

## 
## Call:
##  density.default(x = mtcars$disp)
## 
## Data: mtcars$disp (32 obs.); Bandwidth 'bw' = 55.77
## 
##        x                y            
##  Min.   :-96.22   Min.   :4.209e-06  
##  1st Qu.: 87.67   1st Qu.:2.932e-04  
##  Median :271.55   Median :1.417e-03  
##  Mean   :271.55   Mean   :1.358e-03  
##  3rd Qu.:455.43   3rd Qu.:2.159e-03  
##  Max.   :639.32   Max.   :3.100e-03

modes <- function(d){
  i <- which(diff(sign(diff(d$y))) < 0) + 1
  data.frame(x = d$x[i], y = d$y[i])
}

modes(d)

##          x           y
## 1 128.3295 0.003100294
## 2 305.3759 0.002204658

d$x[which.max(d$y)]    # double-check

## [1] 128.3295

ggplot(data.frame(x = d$x, y = d$y), aes(x, y)) + 
  geom_density(stat = "identity", fill = 'mistyrose', alpha = 0.3) + 
  geom_vline(xintercept = modes(d)$x)

d <- density(diamonds$carat, n = 2048)
modes(d)

##            x            y
## 1  0.3292376 1.7747630053
## 2  0.5160818 0.9802698796
## 3  0.7153823 0.8664724666
## 4  1.0193155 1.0640973891
## 5  1.1937035 0.4611193525
## 6  1.5175667 0.4129392437
## 7  2.0232917 0.1923433463
## 8  2.4966304 0.0135706931
## 9  3.0098291 0.0038717759
## 10 3.2265684 0.0003081410
## 11 3.5006065 0.0003225040
## 12 3.6600469 0.0003016810
## 13 4.0088228 0.0004634136
## 14 4.4996002 0.0001531126
## 15 5.0103077 0.0001532038

ggplot(data.frame(x = d$x, y = d$y * d$n), aes(x, y)) + 
  geom_density(stat = "identity", fill = 'papayawhip', alpha = 0.3) + 
  geom_point(data = modes(d), aes(y = y * d$n))

#ref https://stackoverflow.com/questions/65886842/extracting-peak-values-from-a-density-ggplot
#ref https://stackoverflow.com/questions/58785930/r-find-maximum-of-density-plot