r - 通过距源的距离进行空间聚类的定向测试

Question

我有一个动物位置的空间数据集，作为源周围的 (x,y) 点(半径为 5 公里的圆形图案)。我需要测试这些点是否相对于远离源的位置聚集(或排斥)在源周围，同时考虑方向性。
我尝试过的事情:

1. 检查了最近邻居和 Ripley 的 K - 无法弄清楚如何合并与源的距离或方向性(另外，看起来 Ripley 的窗口是矩形的，而我的是圆形的)
2. 将数据划分为基本方向(N、E、S、W)和距离箱，并计算每个距离/方向箱的点密度。然后又卡住了。

理想情况下，我会得到一个结果，可以告诉我“你的点在 X 方向上像 donut 一样分布，在 Y 方向上是随机的，在 Z 方向上是山峰状”。我感觉像this answer (重采样 + mad.test)对我来说可能是正确的方向，但不能完全适应我的场景......

这是围绕点源的圆形分布的假数据集:

library(ggplot2)
library(spatstat)
library(dplyr)

set.seed(310366)
nclust <- function(x0, y0, radius, n) {
           return(runifdisc(n, radius, centre=c(x0, y0)))
         }

c <- rPoissonCluster(1000, 0.1, nclust, radius=0.02, n=2)

df <- data.frame(x = c$x - 0.5, y = c$y - 0.5) %>%
    mutate(Distance = sqrt(x^2 + y^2)) %>%
    filter(Distance < 0.5)

ggplot(df) + 
    geom_point(aes(x = x, y = y)) +
    # source point
    geom_point(aes(x = 0, y = 0, colour = "Source"), size = 2) +
    coord_fixed()

Answer 1

也许您可以通过研究模式的各向异性来获得一些相关的见解(即使它可能不会为您提供您正在寻找的所有答案)。

检测点图案各向异性的工具包括:扇区 K 函数、对方向分布、各向异性对相关函数。这些都在 Spatial Point Patterns: Methodology and Applications with R 的第 7.9 节中进行了描述。 (免责声明:我是合著者)。幸运的是，第 7 章是免费示例章节之一，因此您可以在这里下载:http://book.spatstat.org/sample-chapters.html .

它不会以特殊方式处理您的源位置，因此它不能解决整个问题，但当您考虑要做什么时，它可能会为您提供灵感。

您可以制作一个泊松模型，其强度取决于距源的距离和方向，看看这是否能给您带来任何见解。

下面是一些轻微注释的代码片段，因为我没有时间详细说明(请记住这些只是粗略的想法 - 可能有更好的替代方案)。请随意改进。

单位圆盘中的均匀点:

library(spatstat)
set.seed(42)
X <- runifdisc(2000)
plot(X)

W <- Window(X)

极坐标作为协变量:

rad <- as.im(function(x,y){sqrt(x^2+y^2)}, W)
ang <- as.im(atan2, W)
plot(solist(rad, ang), main = "")

north <- ang < 45/180*pi & ang > -45/180*pi
east <- ang > 45/180*pi & ang < 135/180*pi
west <- ang < -45/180*pi & ang > -135/180*pi
south <- ang< -135/180*pi | ang > 135/180*pi
plot(solist(north, east, west, south), main = "")

plot(solist(rad*north, rad*east, rad*west, rad*south), main = "")

拟合一个简单的对数线性模型(可以拟合更复杂的关系与ippm():

mod <- ppm(X ~ rad*west + rad*south +rad*east)
mod
#> Nonstationary Poisson process
#> 
#> Log intensity:  ~rad * west + rad * south + rad * east
#> 
#> Fitted trend coefficients:
#>   (Intercept)           rad      westTRUE     southTRUE      eastTRUE 
#>    6.37408999    0.09752045   -0.23197347    0.18205119    0.03103026 
#>  rad:westTRUE rad:southTRUE  rad:eastTRUE 
#>    0.32480273   -0.29191172    0.09064405 
#> 
#>                  Estimate      S.E.    CI95.lo   CI95.hi Ztest       Zval
#> (Intercept)    6.37408999 0.1285505  6.1221355 6.6260444   *** 49.5843075
#> rad            0.09752045 0.1824012 -0.2599794 0.4550203        0.5346480
#> westTRUE      -0.23197347 0.1955670 -0.6152777 0.1513307       -1.1861588
#> southTRUE      0.18205119 0.1870798 -0.1846184 0.5487208        0.9731206
#> eastTRUE       0.03103026 0.1868560 -0.3352008 0.3972613        0.1660651
#> rad:westTRUE   0.32480273 0.2724648 -0.2092185 0.8588240        1.1920904
#> rad:southTRUE -0.29191172 0.2664309 -0.8141066 0.2302832       -1.0956377
#> rad:eastTRUE   0.09064405 0.2626135 -0.4240690 0.6053571        0.3451614
plot(predict(mod))

非统一模型:

lam <- 2000*exp(-2*rad - rad*north - 3*rad*west)
plot(lam)

set.seed(4242)
X2 <- rpoispp(lam)[W]
plot(X2)

适合:

mod2 <- ppm(X2 ~ rad*west + rad*south +rad*east)
plot(predict(mod2))
plot(X2, add = TRUE, col = rgb(.9,.9,.9,.5))

在中心添加点并查看限制为该点的 Ksector() 作为引用点(对于这个例子来说信息不是很多，但可能会有所帮助在其他情况下？？):

X0 <- ppp(0, 0, window = W)
plot(X2[disc(.1)], main = "Zoom-in of center disc(0.1) of X2")
plot(X0, add = TRUE, col = "red")
dom <- disc(.01)
plot(dom, add = TRUE, border = "blue")

X3 <- superimpose(X2, X0)

估计的北区 K 函数高于西区(绘制差异):

Knorth <- Ksector(X3, 45, 135, domain = dom)
Kwest <- Ksector(X3, 135, 225, domain = dom)
plot(eval.fv(Knorth-Kwest), iso~r)

^{由reprex package于2018年12月18日创建(v0.2.1)}