r - 使用 R 显示跨面密度的主要峰值

Question

我正在尝试使用方面内的 ggplot 绘制数据的分布/密度。这是我现在所拥有的，其中红线显示平均值，每个方面显示平均值。现在在这里，平均值没有意义，我希望有类似的绘图，其中密度中的峰值用 xintercept 和文本显示。
我用于手段的代码是这样的:

data <- read.table("sample.csv", header=F, sep=',')
colnames(data) <- c("frame", "val")
attach(data)
library(ggplot2)
library(grid)

library(plyr)
xdat <- ddply(data,"frame", transform, val_mean = signif(mean(val),3), med.x = signif(mean(val),3), med.y=signif(mean(density(val)$y),3))

ppi <- 500
png("sample.png", width=4*ppi, height=4*ppi, res=ppi)

hp <-ggplot(data=data, aes(x=val))+
geom_density() +
geom_vline(aes(xintercept=val_mean),xdat, color="red",linetype="dashed",size=1) +
theme_bw()

hp<-hp + facet_wrap (~ frame, ncol=2, scales="free_y") +
geom_text(data = xdat, aes(x=med.x,y=med.y,label=val_mean))

print(hp)
dev.off()

用于绘制此图的数据是:

data <- data.frame(
    "frame"=c(rep("A",9), rep("B", 13), rep("C", 7)), 
    "val"=c(1, rep(2,4), 4, 5, 6, rep(1,6), 2, rep(3,7), 1, rep(4,6))
    )

我知道有一些帖子使用 R 来查找值的峰值。但是我希望绘制密度的峰值并且我无法找到任何解决方案(或者我可能错过了它)。是否可以在 R 中即时计算峰值并在不同方面进行绘图？
非常感谢您的时间和帮助!!

Answer 1

我假设您想确定每个方面的单个最大峰值 - 这将是分布的模式。如果您的分布是多峰的，我的回答只会识别最大的峰值。 This answer to another question 解释说 geom_density() 使用带有默认参数的 density() 函数。

话虽如此，以下代码应该适合您:

library(ggplot2)
library(grid)
library(plyr)

data <- data.frame("frame"=c(rep("A",9), rep("B", 13), rep("C", 7)), "val"=c(1,rep(2,4),4,5,6,rep(1,6),2,rep(3,7),1,rep(4,6)))
attach(data)

densMode <- function(x){
    td <- density(x)
    maxDens <- which.max(td$y)
    list(x=td$x[maxDens], y=td$y[maxDens])
}
xdat <- ddply(data,"frame", transform, val_mean = signif(densMode(val)$x,3), med.x = signif(densMode(val)$x,3), med.y=signif(densMode(val)$y,3))

hp <- ggplot(data=data, aes(x=val)) + 
    geom_density() + 
    geom_vline(aes(xintercept=val_mean),xdat, color="red",linetype="dashed",size=1) + 
    theme_bw()

hp<- hp + 
    facet_wrap (~ frame, ncol=2, scales="free_y") + 
    geom_text(data = xdat, aes(x=med.x,y=med.y,label=val_mean))

hp

我改变的唯一几行是那些决定图形是如何创建的(我没有使用 png() )，插入 densMode() 函数，并在 densMode() 的定义中使用 xdat 。我还根据您的示例数据创建了一个 data.frame(我已将其作为对您问题的编辑提交，以方便其他可能想要回答的人)。

代码产生下图: