使用 R 将四次核热图转换为大多边形

Posted 2023-02-14

【中文标题】使用 R 将四次核热图转换为大多边形

【英文标题】：Converting quartic kernel heatmap into large polygon with R

【发布时间】：2021-08-28 18:18:17

【问题描述】：

我有欧胡岛海岸附近的点数据。其他人使用这些相同的数据创建了一个大的polygon。我相信他首先使用quartic (biweight) kernel 创建了一个heatmap，每个点的半径为1 公里，像素大小可能为1 平方公里。他引用了 Silverman（1986 年，第 76 页，方程 4.5，我相信它指的是《统计和数据分析的密度估计》一书）。我相信他将他的heatmap 转换为他的polygon。我正在尝试使用R 和Windows 10 用虚假数据来近似他的polygon。我可以使用ks 包中的kde 函数来接近（见下图）。但该软件包仅包含Gaussian kernels。是否可以使用quartic kernel 创建类似的polygon？

另一个分析者实际上创建了polygon 的两个版本。一个的边界被标记为“> 1 per km density”；另一个边界被标记为“> 0.5 每公里密度”。我不知道他是用R、QGIS、ArcGIS还是别的什么。我无法在 QGIS 中创建单个大 polygon 并且没有 ArcGIS。

感谢您就如何创建与所示类似但使用quartic kernel 而不是Gaussian kernel 的polygon 提出任何建议。如果我可以提供更多信息，请告诉我。

这是我的CSV 和QGIS 格式的假数据的链接：enter link description here（编辑：希望现在任何人都可以访问假数据。我以前可以，但我猜其他人不能。）

1. fake_points_oahu.csv

     a. raw data

2. fake_points_oahu_utm (.shp, .dbf, .prj, .shx) 

     a. vector point layer 

3. fake_points_oahu_June11_2021.png

     a. the figure shown above

这是我的R 代码：

setwd('C:/Users/mark_/Documents/ctmm/density_in_R/density_files_for_***/')

library(sf) # to read shapefile
library(ks) # to use kde function

my.data <- read.csv("fake_points_oahu.csv", header = TRUE, stringsAsFactors = FALSE, na.strings = "NA")
head(my.data)

# Import shapefile
st_layers("fake_points_oahu_utm.shp")

points_utm <- st_read(dsn = "fake_points_oahu_utm.shp", layer = 'fake_points_oahu_utm')
st_crs(points_utm)
plot(points_utm)

my.matrix <- as.matrix(my.data[,2:3])
head(my.matrix)

# This uses the Guassian kernel
my_gps_hpi <- Hpi(x = my.matrix, pilot = "samse", pre = "scale")

my.fhat <- kde(x = my.matrix, compute.cont = TRUE, h = my_gps_hpi,
               xmin = c(min(my.data$longitude), min(my.data$latitude)),
               xmax = c(max(my.data$longitude), max(my.data$latitude)),
               bgridsize = c(500, 500))

my.contours <- c(96.5)

contourLevels(my.fhat, cont = my.contours)
contourSizes(my.fhat, cont = my.contours, approx = TRUE)

plot(my.data$longitude, my.data$latitude)
plot(my.fhat, lwd = 3, display = "filled.contour", cont = my.contours, add = TRUE)

png(file="fake_points_oahu_June11_2021.png")

     plot(my.data$longitude, my.data$latitude)
     plot(my.fhat, lwd = 3, display = "filled.contour", cont = my.contours, add = TRUE)

dev.off()

【问题讨论】：

您的样本数据不可访问。

@kwes 我现在使用了 Google Drive 上的复制链接功能，并选择与拥有该链接的任何人共享。然后我在这里复制了那个链接。希望您现在可以访问虚假数据。

【参考方案1】：

您可以通过稍微修改 MASS 包中的 kde2d 函数来执行您的估计。据我所知，R 中目前没有包可以使用四次（双权重）内核实现双变量 KDE 估计。

单变量双权核可以通过多种方式扩展到多变量核，最简单的方法是使用乘积核，您可以对每个维度使用单变量核，然后将结果相乘。您可以找到二重乘积内核here 的数学表达式。 当您将此内核合并到来自MASS 包的kde2d 密度估计器时，它看起来如下所示

kde_biweight_kernel <- function(x,y, bw_x, bw_y, xrange, yrange)
  # This function is based on the kde2d function from 
  # the MASS package. The only difference is that the Gaussian
  # kernel is substituted with a biweight product kernel
  
  # product kernel:
  biweight_kernel <- function(u)
    mask = abs(u) > 1
    kernel_val = (15/16)*((1-u^2)^2)
    kernel_val[mask] = 0
    return(kernel_val)
  
  
  lims = c(xrange, yrange)
  n = 500
  nx <- length(x)
  n <- rep(n, length.out = 2L)
  # get grid on which we want to estimate the density
  gx <- seq.int(lims[1L], lims[2L], length.out = n[1L])
  gy <- seq.int(lims[3L], lims[4L], length.out = n[2L])
  
  # inputs to kernel
  ax <- outer(gx, x, "-" )/bw_x
  ay <- outer(gy, y, "-" )/bw_y
  
  # evaluate and multiply kernel results along both axes
  res = tcrossprod(biweight_kernel(ax), biweight_kernel(ay))/(nx * bw_x * bw_y)
  return(list(x = gx, y = gy, z = res))

使用kde_biweight_kernel函数，您可以计算所需的密度，如下所示

library(MASS)
library(birk)
library(kedd)
library(sf)
library(ks)


# load data
my.data <- read.csv("fake_points_oahu.csv", header = TRUE, stringsAsFactors = FALSE, na.strings = "NA")
# Import shapefile
st_layers("fake_points_oahu_utm.shp")
points_utm <- st_read(dsn = "fake_points_oahu_utm.shp", layer = 'fake_points_oahu_utm')

x = my.data$longitude
y = my.data$latitude

# determine bandwidth for biweight kernel along both axes
bw_x = h.amise(x, deriv.order = 0, kernel = "biweight")$h
bw_y = h.amise(y, deriv.order = 0, kernel = "biweight")$h

# get ranges in which you want to estimate density
xrange = c(min(my.data$longitude), max(my.data$longitude))
yrange = c(min(my.data$latitude),  max(my.data$latitude))

# get 2d density estimate with quartic (biweight) kernel
result = kde_biweight_kernel(x,y, bw_x, bw_y, xrange, yrange)

请注意，带宽是专门针对双权核情况计算的。 生成的密度对象与ks::kde 对象有点不同。例如，它还没有等高线级别。我们可以通过使用rmngbpackage 中的kde2dQuantile 函数的略微修改版本计算分位数来获得等高线水平

# get quantiles of interest:
kde2dQuantile <- function(d, X, Y, probs = .05) 
  xInd <- sapply(X, function(x) which.closest(d$x, x))
  yInd <- sapply(Y, function(x) which.closest(d$y, x))
  zValues <- d$z[cbind(xInd, yInd)]
  quantile(zValues, probs=probs)

# get quantiles
quantiles = kde2dQuantile(result, x, y, seq(0,1,by=0.001))

根据您的问题，我不确定您对哪个分位数感兴趣，所以我只选择了 1% 的分位数。 为了能够以与问题相同的方式绘制数据，我们必须以与 kde 类中的对象相同的方式格式化密度结果：

# to make the kde estimate compatible with the other density estimates
# from the ks package, the result can be converted to a named list.
# -> create ks::KDE object:
axes = matrix(c(result$x,result$y), ncol = 2)
colnames(axes) = c('longitude', 'latitude')

my.fhat_biweight = list('x' = axes,
                        'eval.points' = list(result$x, result$y),
                        'estimate' = result['z']$z,
                        'gridtype' = 'linear', 'gridded' = TRUE,
                        'binned' = TRUE, 'names' = c("longitude","latitude" ))

# add quantile to ks::KDE object
my.fhat_biweight$cont = quantiles

# change class (make sure ks package is loaded for this)
class(my.fhat_biweight) <- "kde"

最后绘制数据上的二权核密度

plot(my.data$longitude, my.data$latitude)
plot(my.fhat_biweight, lwd = 3, display = "filled.contour", cont = cont=c(96.5), add = TRUE)

这个输出：

【讨论】：

我有点担心它在欧胡岛的土地上延伸了多少。高斯核轮廓似乎非常接近海岸线。我不确定我看到你在哪里使用了 1% 分位数，但我需要更多地研究这个。我不确定，但我认为另一位分析师正在根据密度的绝对值（0.5 或 1）绘制他的轮廓，但我可能错了。无论如何，感谢您为此所做的所有工作。

是的，你完全正确，这是错误的链接，我更新了它。乘积核的表达式在 6.3.2 Multivariate Density Estimation 下。

我也不太确定置信水平，并使用了与您在问题中提供的相同的置信水平。还有更复杂的变体可以将单变量双权核扩展到双变量情况 - 乘积核确实是最简单的方法。您可以通过调用my.fhat_biweight$cont 在分位数向量中找到 1% 的分位数值。但是从我对 plot + "display.quantile" 函数的理解来看，cont 参数是一个分位数水平的向量，所以你实际上通过了 96.5% 的分位数？如果我错了，请纠正我！

谢谢你，yuki。您的代码中的哪一行选择了 1% 分位数？我认为这条线会生成所有分位数：quantiles = kde2dQuantile(result, x, y, seq(0,1,by=0.001))。我可能会误会。我使用了 96.5% 的等值线，但怀疑最初的分析师没有。我怀疑他使用了基于密度值> 0.5（或> 1）的网格单元的轮廓。我可能错了。是否可以根据网格单元内的密度值生成多边形？您的代码中的第二个 plot 语句是否应该是：plot(my.fhat_biweight, lwd = 3, display = "filled.contour", cont = c(96.5), add = TRUE)？

我实际上没有选择 1% 的分位数。您可以使用 quantiles["1.0%"] 从分位数列表中获取 1% 的分位数，然后将其传递给图中的 contargument。但是你确定你想要 1% 的分位数吗？也许您可以发布该问题的文献链接，以便“> 1 每公里密度”的描述可能会变得更清楚一点？