Как построить контурную линию, показывающую, где 95% значений попадают внутрь, в R и в ggplot2

Скажем, что у нас есть:

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

Как использовать ggplot2 для создания графика, содержащего две следующие геометрии:

Двумерное ожидание двух рядов значений
Контурная линия, показывающая, где находятся 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))

И я также знаю о функциях stat_contour() и stat_density2d() в ggplot2.

И я также знаю, что есть опции "бункеров" в stat_contour.

Однако, я думаю, что мне нужно что-то вроде аргумента probs в квантиле, но более двух измерений, а не одного.

Я также видел решение в графическом пакете. Тем не менее, я хотел бы сделать это в ggplot.

Помогите оценить,

Джон

Ответ 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.

Ответ 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)

update: с недавней версией 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")

Выполнение вычисления уровня 95% из сетки kk (для уменьшения числа вычислений ядра до 1) остается в виде упражнения
Я не уверен, почему stat_density2d(geom="tile") и geom_tile дают несколько разные результаты (первый сглажен)
Я не добавил двумерное значение, но что-то вроде annotate("point",x=mean(d$x),y=mean(d$y),colour="red") должно работать.

Ответ 3

Не обращая внимания на ответ Бена Болкера, решение, которое может обрабатывать несколько уровней и работает с 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)

Результат:

Ответ 4

У меня был пример, когда спецификации ширины полосы пропускания MASS::kde2d() были недостаточно гибкими, поэтому я закончил использование пакета ks и функции ks::kde() и, в качестве примера, функцию ks::Hscv() для оценки гибкой полосы пропускания что улучшило гладкость. Это вычисление может быть немного медленным, но в некоторых ситуациях оно имеет гораздо лучшую производительность. Вот пример приведенного выше кода для этого примера:

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)

## using ks::kde
hscv1 <- Hscv(d)
fhat <- ks::kde(d, H=hscv1, compute.cont=TRUE)

dimnames(fhat[['estimate']]) <- list(fhat[["eval.points"]][[1]], 
                                     fhat[["eval.points"]][[2]])
library(reshape2)
aa <- melt(fhat[['estimate']])

ggplot(aa, aes(x=Var1, y=Var2)) +
    geom_tile(aes(fill=value)) +
    geom_contour(aes(z=value), breaks=fhat[["cont"]]["50%"], color="red") +
    geom_contour(aes(z=value), breaks=fhat[["cont"]]["5%"], color="purple")

В этом конкретном примере различия минимальны, но в примере, когда спецификация полосы пропускания требует большей гибкости, эта модификация может быть важной. Обратите внимание, что контур 95% задается с помощью breaks=fhat[["cont"]]["5%"], который я нашел немного интуитивно понятным, потому что здесь он называется "контур 5%".

Ответ 5

Просто tidyverse ответы сверху, разместите их в более tidyverse для tidyverse виде и tidyverse несколько уровней контура. Я использую здесь geom_path(group=probs), добавляя их вручную geom_text. Другой подход заключается в использовании geom_path(colour=probs) который автоматически geom_path(colour=probs) контуры как geom_path(colour=probs) обозначения.

library(ks)
library(tidyverse)

set.seed(1001)

## data
d <- MASS::mvrnorm(1000, c(0, 0.2), matrix(c(1, 0.4, 1, 0.4), ncol=2)) %>% 
  magrittr::set_colnames(c("x", "y")) %>% 
  as_tibble() 

## density function
kd <- ks::kde(d, compute.cont=TRUE, h=0.2)

## extract results
get_contour <- function(kd_out=kd, prob="5%") {
  contour_95 <- with(kd_out, contourLines(x=eval.points[[1]], y=eval.points[[2]],
                                      z=estimate, levels=cont[prob])[[1]])
  as_tibble(contour_95) %>% 
    mutate(prob = prob)
}

dat_out <- map_dfr(c("10%", "20%","80%", "90%"), ~get_contour(kd, .)) %>% 
  group_by(prob) %>% 
  mutate(n_val = 1:n()) %>% 
  ungroup()

## clean kde output
kd_df <- expand_grid(x=kd$eval.points[[1]], y=kd$eval.points[[2]]) %>% 
  mutate(z = c(kd$estimate %>% t))

ggplot(data=kd_df, aes(x, y)) +
  geom_tile(aes(fill=z)) +
  geom_point(data = d, alpha = I(0.4), size = I(0.4), colour = I("yellow")) +
  geom_path(aes(x, y, group = prob), 
            data=filter(dat_out, !n_val %in% 1:3), colour = I("white")) +
  geom_text(aes(label = prob), data = 
              filter(dat_out, (prob%in% c("10%", "20%","80%") & n_val==1) | (prob%in% c("90%") & n_val==20)),
            colour = I("black"), size =I(3))+
  scale_fill_viridis_c()+
  theme_bw() +
  theme(legend.position = "none")

^{Создано в 2019-06-25 пакетом представлением (v0.3.0)}