stat_密度2d - 图例是什么意思？

Question

我用 R 制作了一张地图stat_density2d。这是代码：

ggplot(data, aes(x=Lon, y=Lat)) + 
  stat_density2d(aes(fill = ..level..), alpha=0.5, geom="polygon",show.legend=FALSE)+
  geom_point(colour="red")+
  geom_path(data=map.df,aes(x=long, y=lat, group=group), colour="grey50")+
  scale_fill_gradientn(colours=rev(brewer.pal(7,"Spectral")))+
  xlim(-10,+2.5) +
  ylim(+47,+60) +
  coord_fixed(1.7) +
  theme_void()

它产生这个：

伟大的。有用。但我不知道这个传说意味着什么。我确实找到了这个维基百科页面：

https://en.wikipedia.org/wiki/Multivariate_kernel_密度_估计

他们使用的示例（包含红色、橙色和黄色）指出：

彩色轮廓对应于包含相应概率质量的最小区域：红色 = 25%，橙色 + 红色 = 50%，黄色 + 橙色 + 红色 = 75%

然而，使用 stat_密度2d，我的地图中有 11 个等高线。有谁知道 stat_密度2d 是如何工作的以及图例的含义是什么？理想情况下，我希望能够说明诸如红色轮廓包含 25% 的图之类的内容。

我读过这篇文章： https: //ggplot2.tidyverse.org/reference/geom_密度_2d.html，但我仍然一无所知。

Answer 1

让我们以faithfulggplot2 为例：

ggplot(faithful, aes(x = eruptions, y = waiting)) +
  stat_density_2d(aes(fill = factor(stat(level))), geom = "polygon") +
  geom_point() +
  xlim(0.5, 6) +
  ylim(40, 110)

（提前道歉，没有让这个变得更漂亮）

水平面是 3D“山”被切片的高度。我不知道有什么方法（其他人可能）将其转换为百分比，但我确实知道如何得到你所说的百分比。

如果我们查看该图表，水平0.002包含绝大多数点（除了 2 个点之外的所有点）。关卡0.004实际上是 2 个多边形，它们包含除了大约十几个点之外的所有点。如果我明白了您所问的要点，那就是您想知道的，除了不是计数而是给定级别的多边形所包含的点的百分比。使用涉及的各种 ggplot2“统计数据”的方法可以直接计算该值。

请注意，当我们导入tidyverse和sp包时，我们将使用一些其他完全限定的函数。现在，让我们faithful稍微重塑一下数据：

library(tidyverse)
library(sp)

xdf <- select(faithful, x = eruptions, y = waiting)

（更容易输入x和y）

现在，我们将按照 ggplot2 的方式计算二维核密度估计：

h <- c(MASS::bandwidth.nrd(xdf$x), MASS::bandwidth.nrd(xdf$y))

dens <- MASS::kde2d(
  xdf$x, xdf$y, h = h, n = 100,
  lims = c(0.5, 6, 40, 110)
)

breaks <- pretty(range(zdf$z), 10)

zdf <- data.frame(expand.grid(x = dens$x, y = dens$y), z = as.vector(dens$z))

z <- tapply(zdf$z, zdf[c("x", "y")], identity)

cl <- grDevices::contourLines(
  x = sort(unique(dens$x)), y = sort(unique(dens$y)), z = dens$z,
  levels = breaks
)

我不会用输出来混淆答案str()，但看看那里发生的事情很有趣。

我们可以使用空间操作来计算有多少个点落在给定的多边形内，然后我们可以将多边形分组在同一级别以提供每个级别的计数和百分比：

SpatialPolygons(
  lapply(1:length(cl), function(idx) {
    Polygons(
      srl = list(Polygon(
        matrix(c(cl[[idx]]$x, cl[[idx]]$y), nrow=length(cl[[idx]]$x), byrow=FALSE)
      )),
      ID = idx
    )
  })
) -> cont

coordinates(xdf) <- ~x+y

data_frame(
  ct = sapply(over(cont, geometry(xdf), returnList = TRUE), length),
  id = 1:length(ct),
  lvl = sapply(cl, function(x) x$level)
) %>% 
  count(lvl, wt=ct) %>% 
  mutate(
    pct = n/length(xdf),
    pct_lab = sprintf("%s of the points fall within this level", scales::percent(pct))
  )
## # A tibble: 12 x 4
##      lvl     n    pct pct_lab                              
##    <dbl> <int>  <dbl> <chr>                                
##  1 0.002   270 0.993  99.3% of the points fall within this level
##  2 0.004   259 0.952  95.2% of the points fall within this level
##  3 0.006   249 0.915  91.5% of the points fall within this level
##  4 0.008   232 0.853  85.3% of the points fall within this level
##  5 0.01    206 0.757  75.7% of the points fall within this level
##  6 0.012   175 0.643  64.3% of the points fall within this level
##  7 0.014   145 0.533  53.3% of the points fall within this level
##  8 0.016    94 0.346  34.6% of the points fall within this level
##  9 0.018    81 0.298  29.8% of the points fall within this level
## 10 0.02     60 0.221  22.1% of the points fall within this level
## 11 0.022    43 0.158  15.8% of the points fall within this level
## 12 0.024    13 0.0478  4.8% of the points fall within this level 

ggalt::geom_bkde2d()我只是将其拼写出来，以避免更多废话，但百分比会根据您如何修改密度计算的各种参数而变化（对于使用不同估计器的我来说也是如此）。

如果有一种方法可以在不重新执行计算的情况下梳理出百分比，那么没有比让其他 SO R 人员展示他们比写这个答案的人聪明得多（希望以更外交的方式）来指出这一点更好的方法了。方式比似乎是最近的模式）。