13

我最初的目标是绘制单个点的种群，然后绘制一个凸包，包围 80% 的人口，以人口的质量为中心。

在尝试了许多想法之后，我想出的最佳解决方案是使用ggplot's stat_density2d. 虽然这对于定性分析很有用，但我仍然需要指出 80% 的边界。我开始寻找一种方法来勾勒 80% 的人口边界，但我可以以 80% 的概率工作密度边界来代替。

这里是我寻求帮助的地方。（用于）的bin参数没有明确记录。如果我在下面的示例中设置 = 4，我是否正确地将中央（绿色）区域解释为包含 25% 的概率质量，而黄色、红色和绿色的组合区域代表 75% 的概率质量？如果是这样，通过将 bin 更改为 = 5，内接区域是否等于 80% 的概率质量？kde2dstat_density2dbin

set.seed(1)
n=100

df <- data.frame(x=rnorm(n, 0, 1), y=rnorm(n, 0, 1))

TestData <- ggplot (data = df) +
  stat_density2d(aes(x = x, y = y, fill = as.factor(..level..)), 
  bins=4, geom = "polygon", ) +
  geom_point(aes(x = x, y = y)) +
  scale_fill_manual(values = c("yellow","red","green","royalblue", "black"))

TestData

在此处输入图像描述

我重复了一些测试用例并手动计算了排除的点[希望找到一种方法来根据它们所包含的 ..level.. 来计算它们] 但是考虑到数据的随机性（我的真实数据和测试数据）该区域外的点数stat_density2d变化足以需要寻求帮助。

总而言之，是否有一种实用的方法可以在数据框中的中心 80% 的点人口周围绘制一个多边形？或者，除此之外，我是否可以安全使用stat_density2d并将 bin 设置为 5 以产生 80% 的概率质量？

Bryan Hanson 的出色回答消除了我可以bin在stat_density2d. 结果看起来接近bin4 到 6 左右的值，但正如他所说，实际功能是未知的，因此无法使用。

我使用 DWin 接受的答案中提供的 HDRegionplot 来解决我的问题。为此，我从包中添加了重心 ( COGravity) 和多边形中的点 ( )以完成分析。pnt.in.polySDMTools

library(MASS)
library(coda)
library(SDMTools)
library(emdbook)
library(ggplot2)


theme_set(theme_bw(16))
set.seed(1)
n=100

df <- data.frame(x=rnorm(n, 0, 1), y=rnorm(n, 0, 1))

HPDregionplot(mcmc(data.matrix(df)), prob=0.8)
with(df, points(x,y))
ContourLines <- as.data.frame(HPDregionplot(mcmc(data.matrix(df)), prob=0.8))
df$inpoly <- pnt.in.poly(df, ContourLines[, c("x", "y")])$pip

dp <- df[df$inpoly == 1,]
COG100 <- as.data.frame(t(COGravity(df$x, df$y)))
COG80 <- as.data.frame(t(COGravity(dp$x, dp$y)))

TestData <- ggplot (data = df) +
  stat_density2d(aes(x = x, y = y, fill = as.factor(..level..)), 
  bins=5, geom = "polygon", ) +
  geom_point(aes(x = x, y = y, colour = as.factor(inpoly)), alpha = 1) +
  geom_point(data=COG100, aes(COGx, COGy),colour="white",size=2, shape = 4) +
  geom_point(data=COG80, aes(COGx, COGy),colour="green",size=4, shape = 3) +
  geom_polygon(data = ContourLines, aes(x = x, y = y), color = "blue", fill = NA) +
  scale_fill_manual(values = c("yellow","red","green","royalblue", "brown", "black", "white", "black", "white","black")) +
  scale_colour_manual(values = c("red", "black"))
TestData 
nrow(dp)/nrow(df) # actual number of population members inscribed within the 80% probability polgyon

在此处输入图像描述

4

3 回答 3

4

好吧，首先让我说我不完全确定这个答案，这只是部分答案！没有bin参数MASS::kde2d是 所使用的函数stat_density2d。查看帮助页面kde2d和它的代码（只需在控制台中输入函数名称即可看到），我认为bin参数是h（但是这些函数如何知道传递bin给h尚不清楚）。在帮助页面之后，我们看到如果h没有提供，它是由 计算的MASS:bandwidth.nrd。该功能的帮助页面是这样说的：

# The function is currently defined as
function(x)
{
    r <- quantile(x, c(0.25, 0.75))
    h <- (r[2] - r[1])/1.34
    4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
}

基于此，我认为您最后一个问题（“我安全吗……”）的答案绝对是否定的。 r在上面的函数中是你的假设是安全的，但它显然被修改了，所以你不安全。HTH。

额外的想法：您是否有任何证据表明您的代码正在使用您的bins论点？我想知道它是否被忽略了。如果是这样，请尝试通过h，bins看看它是否在听。

于 2013-10-12T13:08:36.510 回答
2

package:emdbook 中的 HPDregionplot 应该这样做。它确实使用了 MASS::kde2d 但它使结果标准化。我认为它的缺点是它需要一个 mcmc 对象。

library(MASS)
library(coda)
HPDregionplot(mcmc(data.matrix(df)), prob=0.8)
with(df, points(x,y))

在此处输入图像描述

于 2013-10-12T13:04:30.127 回答
1

在 42 的答案的基础上，我进行了简化HPDregionplot()以减少依赖关系并消除使用mcmc-objects 的要求。该函数适用于两列data.frame并且不创建中间图。但是请注意，一旦grDevices::contourLines()返回多个轮廓，这种方法就会中断。

hpd_contour <- function (x, n = 50, prob = 0.95, ...) {
  post1 <- MASS::kde2d(x[[1]], x[[2]], n = n, ...)

  dx <- diff(post1$x[1:2])
  dy <- diff(post1$y[1:2])
  sz <- sort(post1$z)
  c1 <- cumsum(sz) * dx * dy

  levels <- sapply(prob, function(x) {
    approx(c1, sz, xout = 1 - x)$y
  })

  as.data.frame(grDevices::contourLines(post1$x, post1$y, post1$z, levels = levels))
}

theme_set(theme_bw(16))
set.seed(1)
n=100

df <- data.frame(x=rnorm(n, 0, 1), y=rnorm(n, 0, 1))
ContourLines <- hpd_contour(df, prob=0.8)

ggplot(df, aes(x = x, y = y)) +
  stat_density2d(aes(fill = as.factor(..level..)), bins=5, geom = "polygon") +
  geom_point() +
  geom_polygon(data = ContourLines, color = "blue", fill = NA) +
  scale_fill_manual(values = c("yellow","red","green","royalblue", "brown", "black", "white", "black", "white","black")) +
  scale_colour_manual(values = c("red", "black"))

在此处输入图像描述

此外，工作流现在很容易扩展到分组数据。

ContourLines <- iris[, c("Species", "Sepal.Length", "Sepal.Width")] %>% 
  group_by(Species) %>% 
  do(hpd_contour(.[, c("Sepal.Length", "Sepal.Width")], prob=0.8))

ggplot(data = iris, aes(x = Sepal.Length, y = Sepal.Width, color = Species)) +
  geom_point(size = 3, alpha = 0.6) +
  geom_polygon(data = ContourLines, fill = NA) +
  guides(color = FALSE) +
  theme(plot.margin = margin())

在此处输入图像描述

于 2019-12-04T09:58:48.017 回答