如何从 ks R 包中估计 kde 对象 95% 轮廓的面积

Question

我正在尝试从 R 中的 ks 包中估计一个 kde 对象的 95% 轮廓的面积。

如果我使用 ks 包中的示例数据集，我将按如下方式创建内核对象：

library(ks)
data(unicef)
H.scv <- Hscv(x=unicef)
fhat <- kde(x=unicef, H=H.scv)

我可以使用绘图功能轻松绘制 25、50、75% 的等高线：

plot(fhat)

但我想估计轮廓内的面积。

我看到一个类似的问题here，但是提出的答案并没有解决问题。

在我的实际应用程序中，我的数据集是动物坐标的时间序列，我想使用二元正态核来测量该动物的家庭范围大小。我正在使用 ks 包，因为它允许使用插件和平滑交叉验证等方法估计内核分布的带宽。

任何帮助将不胜感激！

Answer 1

这里有两种方法。它们在概念上都相当复杂，但实际上代码非常简单。

fhat <- kde(x=unicef, H=H.scv,compute.cont=TRUE)
contour.95 <- with(fhat,contourLines(x=eval.points[[1]],y=eval.points[[2]],
                                     z=estimate,levels=cont["95%"])[[1]])
library(pracma)
with(contour.95,polyarea(x,y))
# [1] -113.677

library(sp)
library(rgeos)
poly <- with(contour.95,data.frame(x,y))
poly <- rbind(poly,poly[1,])    # polygon needs to be closed...
spPoly <- SpatialPolygons(list(Polygons(list(Polygon(poly)),ID=1)))
gArea(spPoly)
# [1] 113.677

说明

首先，kde(...) 函数返回一个kde 对象，它是一个包含 9 个元素的列表。您可以在文档中阅读相关内容，或者您​​可以在命令行中键入 str(fhat)，或者，如果您使用的是 RStudio（强烈推荐），您可以通过在 Environment 选项卡中展开 fhat 对象来查看此内容。

其中一个元素是$eval.points，即评估核密度估计值的点。默认值是在 151 个等间距点进行评估。 $eval.points 本身就是一个列表，在您的情况下是 2 个向量。因此，fhat$eval.points[[1]] 代表“Under-5”沿线的点，fhat$eval.points[[2]] 代表“Ave life exp”沿线的点。

另一个元素是$estimate，它具有核密度的 z​​ 值，在 x 和 y 的每个组合处进行评估。所以$estimate 是一个 151 X 151 矩阵。

如果您用compute.cont=TRUE 调用kde(...)，您会在结果中得到一个额外的元素：$cont，它包含$estimate 中对应于从1% 到99% 的每个百分位数的z 值。

因此，您需要提取对应于 95% 轮廓的 x 和 y 值，并使用它来计算面积。你可以这样做：

fhat <- kde(x=unicef, H=H.scv,compute.cont=TRUE)    
contour.95 <- with(fhat,contourLines(x=eval.points[[1]],y=eval.points[[2]],
                                     z=estimate,levels=cont["95%"])[[1]])

现在，contour.95 的 x 和 y 值对应于fhat 的 95% 轮廓。有（至少）两种方法可以获取该区域。一个使用pracma 包并计算 直接用。

library(pracma)
with(contour.95,polyarea(x,y))
# [1] -113.677

负值的原因与 x 和 y 的顺序有关：polyarea(...) 将多边形解释为“洞”，因此它的面积为负。

另一种方法是使用rgeos（GIS 包）中的面积计算例程。不幸的是，这需要您首先将您的坐标转换为“SpatialPolygon”对象，这有点像熊。不过，它也很简单。

library(sp)
library(rgeos)
poly <- with(contour.95,data.frame(x,y))
poly <- rbind(poly,poly[1,])    # polygon needs to be closed...
spPoly <- SpatialPolygons(list(Polygons(list(Polygon(poly)),ID=1)))
gArea(spPoly)
# [1] 113.677

Answer 2

另一种方法是使用 kde 包中的contourSizes() 函数。我也有兴趣使用这个包来比较生态学中的 2D 和 3D 空间使用，但我不确定如何提取 2D 密度估计。我通过估计一个“动物”的面积来测试这种方法，该面积仅限于一个已知半径的圆的面积。下面是代码：

set.seed(123)
require(GEOmap)
require(kde)
# need this library for the inpoly function

# Create a data frame centered at  coordinates 0,0
data = data.frame(x=0,y=0)

# Create a vector of radians from 0 to 2*pi for making a circle to
# test the area
circle = seq(0,2*pi,length=100) 

# Select a radius for your circle
radius = 10
# Create a buffer for when you simulate points (this will be more clear below)
buffer = radius+2

# Simulate x and y coordinates from uniform distribution and combine
# values into a dataframe

createPointsX = runif(1000,min = data$x-buffer, max = data$x+buffer)
createPointsY = runif(1000,min = data$y-buffer, max = data$y+buffer)
data1 = data.frame(x=createPointsX,y=createPointsY)

# Plot the raw data
plot(data1$x,data1$y)

# Calculate the coordinates used to create a cirle with center 0,0 and
# with radius specified above
coords = as.data.frame(t(rbind(data$x+sin(circle)*radius,
                           data$y+cos(circle)*radius)))
names(coords) = c("x","y")

# Add circle to plot with red line
lines(coords$x,coords$y,col=2,lwd=2)

# Use the inpoly function to calculate whether points lie within
# the circle or not.
inp = inpoly(data1$x, data1$y, coords)
data1 = data1[inp == 1,]

# Finally add points that lie with the circle as blue filled dots
points(data1$x,data1$y,pch=19,col="blue")

# Radius of the circle (known area)
pi * radius^2
#[1] 314.1593


# Sub in your own data here to calculate 95% homerange or 50% core area usage
H.pi = Hpi(data1,binned=T)
fhat = kde(data1,H=H.pi)
ct1 = contourSizes(fhat, cont = 95, approx=TRUE)

# Compare the known area of the circle to the 95% contour size
ct1
#     5% 
# 291.466 

我还尝试创建 2 个未连接的圆圈并测试 contourSizes() 函数，它似乎在脱节的分布上运行良好。