虽然之前已经有过关于这个问题的讨论，但我想我会在这里发布一个答案，说明如何确保每个轮廓线中包含特定比例的点。
为此，我们可以使用MASS::kde2d获取2d密度，然后使用terra转换为栅格。然后，我们可以根据相关2d密度格网中的密度对点进行排序，并使用approx找到通过分位数的密度

density_quantiles <- function(x, y, quantiles) {
  dens <- MASS::kde2d(x, y, n = 500)
  df   <- cbind(expand.grid(x = dens$x, y = dens$y), z = c(dens$z))
  r    <- terra::rast(df)
  ind  <- sapply(seq_along(x), function(i) cellFromXY(r, cbind(x[i], y[i])))
  ind  <- ind[order(-r[ind][[1]])]
  vals <- r[ind][[1]]
  ret  <- approx(seq_along(ind)/length(ind), vals, xout = quantiles)$y
  replace(ret, is.na(ret), max(r[]))
}

这意味着，如果我们可以指定一个分位数向量：

quantiles <- c(0, 0.2, 0.4, 0.6, 0.8)

我们可以得到一个易于解释的图，显示包含20%、40%、60%和80%的点的区域，如下所示：

ggplot(dat, aes(X, Y)) +
  geom_density2d_filled(
    aes(fill = after_stat(level)),
    contour_var = "density",
    breaks = density_quantiles(dat$X, dat$Y, quantiles)) +
  geom_point(alpha = 0.1) +
  coord_equal() +
  scale_fill_viridis_d('Quantiles', l
                       abels = scales::percent(quantiles[-1]),
                       direction = -1)

四分位数是这样的：

quartiles <- c(0, 0.25, 0.5, 0.75)

ggplot(dat, aes(X, Y)) +
  geom_density2d_filled(
    aes(fill = after_stat(level)),
    contour_var = "density",
    breaks = density_quantiles(dat$X, dat$Y, quartiles)) +
  geom_point(alpha = 0.1) +
  coord_equal() +
  scale_fill_viridis_d('Quartiles', labels = scales::percent(quartiles[-1]),
                       direction = -1)

请注意，这与ndensity图中的水平非常不同，后者是 * 最大密度的比例 *，而不是包含固定比例点的区域。换句话说，如果您从侧面观察ndensity图的3-D表示，条带将全部具有相同的高度，如以下动画所示（请参阅脚注了解生成此动画的代码）：

当存在高密度囊袋时，代码似乎给予了“奇怪”的结果，例如：

set.seed(123)
dat <-
  data.frame(
    X = c(rnorm(300, 3, 2.5), rnorm(150, 7, 2), rnorm(450, 4, 0.5)),
    Y = c(rnorm(300, 6, 2.5), rnorm(150, 2, 2), rnorm(450, 5, 0.5)))

ggplot(dat, aes(X, Y)) +
  geom_density2d_filled(
    aes(fill = after_stat(level)),
    contour_var = "density",
    breaks = density_quantiles(dat$X, dat$Y, quantiles)) +
  geom_point(alpha = 0.1) +
  coord_equal() +
  scale_fill_viridis_d('Quartiles', labels = scales::percent(quantiles[-1]),
                       direction = -1)

我们可以看到第80百分位数包含“岛屿”或非连续区域。然而，这仅仅是因为我们仍在绘制 * 密度 *，这些区域高于包含正确点数的阈值密度值。无论在何处设置阈值，都不能保证密度带是单个连续区域。
我们可以使用density_quantiles的quantiles版本的3D图清楚地看到这一点，其中低密度的小“块”在这里和那里突破了我们的阈值。

dens <- MASS::kde2d(dat$X, dat$Y, n = 1000)

levels <- as.character(cut(dens$z[-1, -1], 
                           breaks = density_quantiles(dat$X, dat$Y, quantiles),
                           labels = c(scales::viridis_pal()(4))))

persp(dens, col = levels, phi = 20, theta = -20, axes = FALSE, border = NA)

相反，如果您希望区域是连续的，那么您的问题将变得难以定义。例如，如果您希望区域是连续的，则不清楚下面的密度图应该是什么样子：

set.seed(123)
dat <-
  data.frame(
    X = c(rnorm(150, 3, 2), rnorm(150, 10, 2)),
    Y = c(rnorm(150, 3, 2), rnorm(150, 10, 2)))

ggplot(dat, aes(X, Y)) +
  geom_density2d_filled(
    aes(fill = after_stat(level)),
    contour_var = "density",
    breaks = density_quantiles(dat$X, dat$Y, quartiles)) +
  geom_point(alpha = 0.1) +
  coord_equal() +
  scale_fill_viridis_d('Quartiles', labels = scales::percent(quartiles[-1]),
                       direction = -1)

在这个例子中，为了得到连续的区域，我们必须选择其中一个聚类作为中心点，从中心点开始计算点的数量，这将给予非常人为和误导性的结果，这些结果根本不能反映密度，而必须是距离最高密度点的距离函数，因此总是一组嵌套的圆，就像靶心一样。

动画代码

library(magick)

p1 <- ggplot(dat, aes(X, Y)) +
  geom_density2d_filled(
    aes(fill = after_stat(level)),
    contour_var = "ndensity",
    breaks = seq(0.25, 1, length.out = 4),
    show.legend = FALSE
  ) +
  geom_point(alpha = 0.1) +
  coord_fixed(0.95, expand = FALSE)  +
  theme(plot.margin = margin(75, 50, 60, 50))

ggsave('gg.png', p1, width = 480, height = 480, units = 'px', dpi = 72)

dens <- MASS::kde2d(dat$X, dat$Y, n = 1000)
df2 <- data.frame(x = dens$x, y = apply(dens$z, 1, max)/max(dens$z))

p2 <- ggplot(df2, aes(x, y)) +
  geom_area(fill = scales::viridis_pal()(3)[3]) +
  geom_area(fill = scales::viridis_pal()(3)[2],
            aes(y = ifelse(y > 0.75, 0.75, y))) +
  geom_area(fill = scales::viridis_pal()(3)[1],
            aes(y = ifelse(y > 0.5, 0.5, y))) +
  geom_area(fill = 'gray92',
            aes(y = ifelse(y > 0.25, 0.25, y))) +
  coord_fixed(diff(range(dens$x))) +
  geom_hline(yintercept = c(0, 0.25, 0.5, 0.75, 1)) +
  theme_classic() +
  theme(plot.margin = margin(58, 50, 50, 45))

ggsave('gg2.png', p2, width = 480, height = 480, units = 'px', dpi = 72)

levels <- as.character(cut(dens$z[-1, -1], 
                           breaks = c(0, 0.25, 0.5, 0.75, 1) * max(dens$z),
                           labels = c('gray92', scales::viridis_pal()(3))))

for(i in seq(0, 90, 3)) {
  ragg::agg_png(paste0("persp", sprintf("%02d", i), ".png"))
  persp(dens, col = levels, phi = i, d = 1000, axes = FALSE, 
        box = FALSE, border = NA)
  dev.off()
}

f <- list.files(pattern = 'persp\\d+\\.png', full.names = TRUE) 

c(c(rep(f[31], 10), rev(f), rep(f[1], 10),
    rep('gg2.png', 10), rep(f[1], 10), f, rep(f[31], 5), 
    rep('gg.png', 10))
  ) %>%
  rev(.) %>%
  image_read() %>% 
  image_join() %>% 
  image_animate(fps = 10) %>% 
  image_write("D:\\persp.gif")