在散点图中绘制95%置信区间限制。

Question

我需要绘制几个数据点，它们被定义为

c(x,y, stdev_x, stdev_y)

作为一个散点图，显示它们的95%置信区间的表示，例如显示该点和周围的一个轮廓。理想情况下，我想在点周围绘制一个椭圆形，但不知道如何实现。我考虑构建样本并绘制它们，添加stat_density2d()，但需要将轮廓数量限制为1，并且无法弄清楚如何实现。

require(ggplot2)
n=10000
d <- data.frame(id=rep("A", n),
                se=rnorm(n, 0.18,0.02), 
                sp=rnorm(n, 0.79,0.06) )
g <- ggplot (d, aes(se,sp)) +
  scale_x_continuous(limits=c(0,1))+
  scale_y_continuous(limits=c(0,1)) +
  theme(aspect.ratio=0.6)
g + geom_point(alpha=I(1/50)) +
  stat_density2d()

- K Owen - Reinstate Monica

Answer 1

首先，将所有绘图保存为对象（更改限制）。

g <- ggplot (d, aes(se,sp, group=id)) +
  scale_x_continuous(limits=c(0,0.5))+
  scale_y_continuous(limits=c(0.5,1)) +
  theme(aspect.ratio=0.6) + 
  geom_point(alpha=I(1/50)) +
  stat_density2d()

使用ggplot_build()函数保存绘图所用的所有信息。等高线存储在对象data[[2]]中。

gg<-ggplot_build(g)
str(gg$data)
head(gg$data[[2]])
  level         x         y piece group PANEL
1    10 0.1363636 0.7390318     1   1-1     1
2    10 0.1355521 0.7424242     1   1-1     1
3    10 0.1347814 0.7474747     1   1-1     1
4    10 0.1343692 0.7525253     1   1-1     1
5    10 0.1340186 0.7575758     1   1-1     1
6    10 0.1336037 0.7626263     1   1-1     1

总共有12条轮廓线，但如果只保留外部轮廓线，您应该只对group == "1-1" 进行子集剪裁，并替换原始信息。

gg$data[[2]]<-subset(gg$data[[2]],group=="1-1")

然后使用ggplot_gtable()和grid.draw()来获取您的图形。

p1<-ggplot_gtable(gg)
grid.draw(p1)

- Didzis Elferts

---

谢谢你的答复。我们怎么知道外轮廓是95%，而不是97%或99%呢？这可能很明显，但我在文档中（包括kde2d的文档）没有找到相关说明。 - K Owen - Reinstate Monica

现在这个解决方案只展示了如何保留一个轮廓线。关于置信度还需要进一步探讨。 - Didzis Elferts

Answer 2

latticeExtra 提供了 panel.ellipse 函数，它可以从二元数据中计算并绘制置信度椭球，可能还会按第三个变量进行分组。

在这里，我使用您的数据绘制水平为 0.65 和 0.95 的椭圆。

library(latticeExtra)
xyplot(sp~se,data=d,groups=id,
       par.settings = list(plot.symbol = list(cex = 1.1, pch=16)),
       panel = function(x,y,...){
         panel.xyplot(x, y,alpha=0.2)
         panel.ellipse(x, y, lwd = 2, col="green", robust=FALSE,  level=0.65,...)
         panel.ellipse(x, y, lwd = 2, col="red", robust=TRUE,  level=0.95,...)

       })

- agstudy

Answer 3

看起来你找到的 stat_ellipse 函数确实是一个很好的解决方案，但这里还有另一个（非 ggplot）的解决方案，只是为了记录，使用 car 包中的 dataEllipse。

# some sample data
n=10000
g=4
d <- data.frame(ID = unlist(lapply(letters[1:g], function(x) rep(x, n/g))))
d$x <- unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))) 
d$y <- unlist(lapply(1:g, function(i) rnorm(n/g, runif(1)*i^2))) 

# plot points with 95% normal-probability contour
# default settings...
library(car)
with(d, dataEllipse(x, y, ID, level=0.95, fill=TRUE, fill.alpha=0.1))

# with a little more effort...
# random colours with alpha-blending
d$col <- unlist(lapply(1:g, function (x) rep(rgb(runif(1), runif(1), runif(1), runif(1)),n/g)))
# plot points first
with(d, plot(x,y, col=col, pch="."))
# then ellipses over the top
with(d, dataEllipse(x, y, ID, level=0.95, fill=TRUE, fill.alpha=0.1, plot.points=FALSE, add=TRUE,  col=unique(col), ellipse.label=FALSE, center.pch="+"))

- Ben

Answer 4

刚刚发现了这个函数stat_ellipse()在这里(还有这里)，它可以很好地处理这个问题。

g + geom_point(alpha=I(1/10))  +
  stat_ellipse(aes(group=id), color="black")

当然，数据集不同：

- K Owen - Reinstate Monica

Answer 5

我不了解ggplot2库，但是您可以使用plotrix绘制椭圆。这个图形是否与您要求的类似?

library(plotrix)
n=10
d <- data.frame(x=runif(n,0,2),y=runif(n,0,2),seX=runif(n,0,0.1),seY=runif(n,0,0.1))
plot(d$x,d$y,pch=16,ylim=c(0,2),xlim=c(0,2))
draw.ellipse(d$x,d$y,d$seX,d$seY)

- slammaster

---

哇，这是一个优雅的解决方案，非常感谢！不过我仍然希望能够在ggplot2中实现这个。 - K Owen - Reinstate Monica