Направленное тестирование пространственной кластеризации по расстоянию от источника
У меня есть пространственный набор данных местоположений животных, поскольку (x,y) указывает вокруг источника (круговая схема с радиусом 5 км). Мне нужно проверить, сгруппированы ли точки (или отталкиваются) вокруг источника относительно источника дальше, учитывая при этом направленность.
Вещи, которые я пытался:
- проверил ближайшего соседа и K Рипли - не мог понять, как включить расстояние от источника или направленность (плюс, похоже, что Рипли ожидает прямоугольное окно, в то время как мой круговой)
- разделить данные на кардинальные направления (N, E, S, W) и ячейки расстояния и рассчитать плотность точек на ячейку расстояния / направления. Затем снова застрял.
В идеале я бы получил результат, который мог бы сказать мне: "Ваши точки распределены как пончик в направлении X, случайны в направлении Y и подобны горным вершинам в направлении Z". Я чувствую, что этот ответ (пересчет + 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()
1 ответ
Может быть, вы можете получить какое-то уместное понимание, просто изучив анизотропию шаблона (хотя, возможно, он не даст вам все ответы, которые вы ищете).
Инструменты для обнаружения анизотропии в точечных рисунках включают: K-функцию сектора, распределение парной ориентации, функцию корреляции анизотропных пар. Все это описано в разделе 7.9 " Пространственные точечные шаблоны: методология и приложения с R" (заявление об отказе: я соавтор). К счастью, глава 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)
Создано 2018-12-18 пакетом представлением (v0.2.1)