如何绘制轮廓线,显示95%的值落在R和ggplot2中的位置

Question

说我们有:

x <- rnorm(1000)
y <- rnorm(1000)

如何使用ggplot2生成包含以下两个geom的图:

1. 两个系列值的双变量期望
2. 一条轮廓线显示95%的估计值在哪里？

我知道如何做第一部分:

 df <- data.frame(x=x, y=y)
 p <- ggplot(df, aes(x=x, y=y))
 p <- p + xlim(-10, 10) + ylim(-10, 10) # say
 p <- p + geom_point(x=mean(x), y=mean(y))

我也知道ggplot2中的stat_contour()和stat_density2d()函数.

我也知道stat_contour中有'bins'选项.

但是,我想我需要的是像分位数中的probs参数,但是超过两个维度而不是一个维度.

我也在图形包中看到了一个解决方案.但是,我想在ggplot中这样做.

非常感谢,

乔恩

Answer 1

不幸的是,接受的答案目前失败Error: Unknown parameters: breaks的ggplot2 2.1.0.我根据本答案中的代码拼凑了一种替代方法,该方法使用该ks包来计算内核密度估计:

library(ggplot2)

set.seed(1001)
d <- data.frame(x=rnorm(1000),y=rnorm(1000))

kd <- ks::kde(d, compute.cont=TRUE)
contour_95 <- with(kd, contourLines(x=eval.points[[1]], y=eval.points[[2]],
                                    z=estimate, levels=cont["5%"])[[1]])
contour_95 <- data.frame(contour_95)

ggplot(data=d, aes(x, y)) +
  geom_point() +
  geom_path(aes(x, y), data=contour_95) +
  theme_bw()

这是结果:

提示:ks包取决于rgl包,这可能是手动编译的痛苦.即使您使用的是Linux,也可以更轻松地获得预编译版本,例如sudo apt install r-cran-rgl在Ubuntu上设置相应的CRAN存储库.

Answer 2

这可行,但效率很低,因为您实际上必须计算三次内核密度估计值.

set.seed(1001)
d <- data.frame(x=rnorm(1000),y=rnorm(1000))
getLevel <- function(x,y,prob=0.95) {
    kk <- MASS::kde2d(x,y)
    dx <- diff(kk$x[1:2])
    dy <- diff(kk$y[1:2])
    sz <- sort(kk$z)
    c1 <- cumsum(sz) * dx * dy
    approx(c1, sz, xout = 1 - prob)$y
}
L95 <- getLevel(d$x,d$y)
library(ggplot2); theme_set(theme_bw())
ggplot(d,aes(x,y)) +
   stat_density2d(geom="tile", aes(fill = ..density..),
                  contour = FALSE)+
   stat_density2d(colour="red",breaks=L95)

(在http://comments.gmane.org/gmane.comp.lang.r.ggplot2/303的帮助下)

更新:使用最新版本的ggplot2(2.1.0)它似乎无法传递breaks给stat_density2d(或至少我不知道如何),但下面的方法geom_contour仍然似乎工作...

通过计算一次核密度估计并绘制同一网格中的切片和轮廓,可以使事情更有效:

kk <- with(dd,MASS::kde2d(x,y))
library(reshape2)
dimnames(kk$z) <- list(kk$x,kk$y)
dc <- melt(kk$z)
ggplot(dc,aes(x=Var1,y=Var2))+
   geom_tile(aes(fill=value))+
   geom_contour(aes(z=value),breaks=L95,colour="red")

• 从kk网格执行95%级别的计算(将内核计算的数量减少到1)留作练习
• 我不知道为什么stat_density2d(geom="tile"),并geom_tile给出的结果稍有不同(前者是平滑)
• 我没有添加双变量的意思,但是annotate("point",x=mean(d$x),y=mean(d$y),colour="red")应该有用.

Answer 3

重复一下Ben Bolker的答案,一个可以处理多个级别并与ggplot 2.2.1一起工作的解决方案:

library(ggplot2)
library(MASS)
library(reshape2)
# create data:
set.seed(8675309)
Sigma <- matrix(c(0.1,0.3,0.3,4),2,2)
mv <- data.frame(mvrnorm(4000,c(1.5,16),Sigma))

# get the kde2d information: 
mv.kde <- kde2d(mv[,1], mv[,2], n = 400)
dx <- diff(mv.kde$x[1:2])  # lifted from emdbook::HPDregionplot()
dy <- diff(mv.kde$y[1:2])
sz <- sort(mv.kde$z)
c1 <- cumsum(sz) * dx * dy

# specify desired contour levels:
prob <- c(0.95,0.90,0.5)

# plot:
dimnames(mv.kde$z) <- list(mv.kde$x,mv.kde$y)
dc <- melt(mv.kde$z)
dc$prob <- approx(sz,1-c1,dc$value)$y
p <- ggplot(dc,aes(x=Var1,y=Var2))+
  geom_contour(aes(z=prob,color=..level..),breaks=prob)+
  geom_point(aes(x=X1,y=X2),data=mv,alpha=0.1,size=1)
print(p)

结果: