stat_density2d - Что означает легенда?

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_density_estimation

И пример, который они использовали (который содержит красный, оранжевый и желтый), гласил:

Цветные контуры соответствуют наименьшей области, которая содержит соответствующая масса вероятности: красный = 25%, оранжевый + красный = 50%, желтый + оранжевый + красный = 75%

Однако, используя stat_density2d, у меня 11 контуров на карте. Кто-нибудь знает, как работает stat_density2d и что означает легенда? В идеале я хотел бы иметь возможность утверждать, что красный контур содержит 25% графиков и т. Д.

Я прочитал это: https://ggplot2.tidyverse.org/reference/geom_density_2d.html и я все еще не мудрый.

Answer 1

Давайте возьмем пример faithful из ggplot2:

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)

(заранее извиняюсь за то, что не сделал это красивее)

Уровень - это высота, на которой были вырезаны трехмерные "горы". Я не знаю, как (другие могут) перевести это в процент, но я знаю, как получить эти проценты.

Если мы посмотрим на этот график, уровень 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 показать, насколько они умнее, чем человек, пишущий этот ответ (надеюсь, более дипломатичным образом, чем кажется в последнее время).