【问题标题】：Identify spikes/peaks in density plot by group按组识别密度图中的尖峰/峰值
【发布时间】：2016-02-01 15:37:42
【问题描述】：

我使用ggplot2 包为R 创建了一个密度图。我想确定图中出现在 0.01 和 0.02 之间的峰值/峰值。传说太多无法挑出，所以我删除了所有传说。我试图过滤掉我的数据，以找到一个组在 0.01 和 0.02 之间的大多数行数。然后我过滤掉选定的组以查看尖峰/峰值是否消失但没有，它仍然在那里绘制。您能否提出一种方法来识别这些图中的这些尖峰/峰值？

这里是一些代码：

ggplot(NumofHitsnormalized, aes(NumofHits_norm, fill = name)) + geom_density(alpha=0.2) + theme(legend.position="none") + xlim(0.0 , 0.15) 

## To filter out the data that is in the range of first spike
test <- NumofHitsnormalized[which(NumofHitsnormalized$NumofHits_norm > 0.01 & NumofHitsnormalized$NumofHits_norm <0.02),] 

## To figure it out which group (name column) has the most number of rows ##thus I thought maybe I could get the data that lead to spike
testMatrix <- matrix(ncol=2, nrow= length(unique(test$name))) 
for (i in 1:length(unique(test$name))){ 
testMatrix[i,1] <- unique(test$name)[i] 
testMatrix[i,2] <- nrow(unique(test$name)[i])}

康拉德，

这是我用极值包过滤掉我的数据后制作的新图。有新的峰值，它们位于不同的间隔，它还说 96% 的初始组在新图中有数据（尽管过滤数据中的行数减少到初始数据集的 0.023%）所以我无法确定哪个峰属于哪些组。

【问题讨论】：

我质疑“异常值”这个词。我似乎隐约记得一位著名的统计学家说过这样的话：“每当我看到一个非常出乎意料的结果时，我都不确定是否应该将其丢弃或申请专利。”。
当然可以显示一些代码？？？
看看extremevalues (PDF) 包，它应该能满足您的需求。
添加到您的帖子中，而不是在评论部分中
你也可以看看 pastecs 包中的 turnpoints()

标签： r ggplot2 density-plot

【解决方案1】：

我遇到了类似的问题。

我的做法是用 3 个窗口创建 y 值的滚动平均值和 sd。

计算你的基线数据的平均标准差（你知道的数据不会有峰值）

设置阈值

如果高于阈值，则为 1，否则为 0。

d5$roll_mean = runMean(d5$`Current (pA)`,n=3)
d5$roll_sd = runSD(x = d5$`Current (pA)`,n = 3)
d5$delta = ifelse(d5$roll_sd>1,1,0)
currents = subset(d5,d5$delta==1,na.rm=TRUE) # Finds all peaks

我的阈值是 sd > 1。根据您的数据，您可能希望使用均值或 sd。对于缓慢上升的峰值，意味着比 sd 更好。

【讨论】：

这是caTools包runmean函数吗？
是的，zoo 包中还有 rollmean，或者您可以很容易地自己编写代码。任你选！

【解决方案2】：

没有看代码，我草拟了这个简单的函数，将TRUE/FALSE 标志添加到表示异常值的变量中：

GenerateOutlierFlag <- function(x) {
  # Load required packages
  Vectorize(require)(package = c("extremevalues"), char = TRUE)
  # Run check for ouliers
  out_flg <- ifelse(1:length(x) %in% getOutliers(x, method = "I")$iLeft,
                    TRUE,FALSE)
  out_flg <- ifelse(1:length(x) %in% getOutliers(x, method = "I")$iRight,
                    TRUE,out_flg)
  return(out_flg)
}

如果您愿意阅读 extremevalues 软件包，您会发现它在识别异常值方面提供了一些灵活性，但从广义上讲，它是查找各种 峰值的好工具em> 或 spikes 在数据中。

侧点

您实际上可以通过创建一个与getOutliers(x, method = "I") 对应的对象来显着优化它，而不是调用该方法两次。

更合理的语法

GenerateOutlierFlag <- function(x) {
  # Load required packages
  require("extremevalues")
  # Outliers object
  outObj <- getOutliers(x, method = "I")
  # Run check for ouliers
  out_flg <- ifelse(1:length(x) %in% outObj$iLeft,
                    TRUE,FALSE)
  out_flg <- ifelse(1:length(x) %in% outObj$iRight,
                    TRUE,out_flg)
  return(out_flg)
}

结果

x <- c(1:10, 1000000, -99099999)
table(GenerateOutlierFlag(x))
FALSE  TRUE 
   10     2

【讨论】：