Проверка гипотез на основе симуляции на пространственных точечных гиперкартах с использованием функции "envelope" в spatstat
Я хочу использовать реплицированные пространственные точечные шаблоны для проверки гипотез в spatstat
, spatstat
имеет прекрасную документацию, и вы можете найти подробную информацию об анализе точечных образцов здесь: https://cran.r-project.org/web/packages/spatstat/vignettes/replicated.pdf
Несколько картированных лесных участков будут представлять собой точечный процесс, работающий в разных местах. Каждый точечный образец отмечен, где каждый знак представляет разновидность дерева. Кроме того, каждый точечный шаблон связан с ковариатным растровым изображением. Во-первых, я хочу создать нулевую модель, которая предполагает, что каждая метка имеет разные отношения с ковариатой. Затем я хочу использовать эту нулевую модель, чтобы проверить, связаны ли определенные виды (или избегать) друг с другом, моделируя гипотезу случайной маркировки и нанося на карту полученные конверты моделирования. (Гипотеза о случайной маркировке гласит, что метки случайным образом назначаются точкам в схеме точек.)
Сначала я покажу, как я обычно делаю этот анализ, используя единый точечный паттерн. Затем я объясню две проблемы, которые у меня возникают с использованием гиперкадров для выполнения одного и того же анализа.
Допустим, у вас есть отмеченный точечный рисунок, где каждая отметка является видом дерева. Я буду использовать набор данных "lansing", доступный в spatstat:
library(spatstat)
data(lansing)
par(mar=rep(0.5,4))
plot(split(lansing),main="")
Теперь предположим, что вы хотите взглянуть на некоторый пространственный ковариат (например, питательные вещества или влажность почвы), поэтому вы создаете растровое изображение сглаженной плотности измерения:
sim1 <- rpoispp(function(x,y) {500 * exp(-3*x)}, win=owin(c(0,1),c(0,1)))
sim1 <- density(sim1)
Сначала создайте нулевую модель:
single.mod <- ppm(lansing ~ marks*sim1)
Где функция "ppm" распознает "метки" как столбец с названиями видов, а "sim1" как ковариату.
Затем выполните проверку гипотезы на основе моделирования, где вас интересует, находятся ли черный дуб и клен в одних и тех же местах.
single.E<-envelope.ppm(single.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(lansing)))
par(mar=rep(4,4))
plot(E,legend=FALSE,ylab="L-function",xlab="Spatial scale (m)",main="Testing random label hypothesis \nfor single point pattern")
Это отлично работает. Теперь, если мы выйдем и попробуем еще пару графиков, чтобы сделать наш анализ более надежным, мы можем включить дополнительные графики в гиперкадр, где каждый график имеет свой точечный шаблон и считается "экспериментальной" копией. Каждый участок также получит свою пространственную ковариату:
sim2 <- rpoispp(function(x,y) {500 * exp(-2*y)}, win=owin(c(0,1),c(0,1)))
sim2 <- density(sim2)
sim3 <- rpoispp(500, win=owin(c(0,1),c(0,1)))
sim3 <- density(sim3)
hyper <- hyperframe(pp=list(lansing,lansing,lansing),sims=list(sim1,sim2,sim3))
Sim2 и sim3 являются пространственными ковариатами для двух дополнительных собранных нами графиков, а функция "гиперкадра" объединяет три точечных шаблона с соответствующими пространственными ковариатами в один гиперкадр.
Я хотел бы построить модель, используя "mppm" (используется для создания моделей для нескольких точечных моделей), где каждый точечный шаблон объясняется их пространственными ковариатами, "sims":
hyper.mod <- mppm(pp ~ sims, data = hyper)
Первая проблема возникает, когда я пытаюсь разрешить каждому знаку иметь разные отношения с ковариатой:
int.mod <- mppm(pp ~ marks*sims, data=hyper)
Выдается следующее сообщение об ошибке:
Ошибка в checkvars(формула, data.sumry$col.names, extra = c("x", "y",: переменная "marks" в формуле не является одним из имен в данных "
Я получаю ту же ошибку, используя:
int.mod <- mppm(pp ~ pp$marks*sims, data=hyper)
Вторая проблема - выяснить, как запустить тестирование гипотез на основе симуляции на гиперкадре. Давайте используем модель гиперкадра, которая работала (hyper.mod), чтобы попробовать это:
E.hyper <- envelope(hyper.mod,Lcross,i="blackoak",j="maple",nsim=39, nrank=1,global=TRUE,correction="best",simulate=expression(rlabel(pp)))
Вы получаете сообщение об ошибке:
Ошибка в UseMethod("конверт"): нет применимого метода для "конверта", примененного к объекту класса "c (" mppm "," list ")"
Подразумевается, что "конверт" не работает на объектах mppm (только ppp или ppm). Я подозреваю, что есть способ обойти это ограничение, но я еще не нашел его. Любые предложения или рекомендации будут очень полезны!
2 ответа
Здесь нет envelope.mppm
в spatstat
потому что некоторые статистические вопросы (связанные с проверкой нескольких гипотез) еще не решены.
Самое быстрое решение, вероятно, использовать cdf.test.mppm
который выполнит тест и даст некоторый графический вывод.
Было бы разумно объединить огибающие (скажем, из K- функций) из разных точечных рисунков, чтобы получить единую огибающую, при условии, что подобранная модель подразумевает, что разные рисунки должны быть статистически эквивалентными. Это не будет действительным, если, например, модель включает в себя ковариаты, которые различны для разных шаблонов.
Лучшей стратегией, вероятно, является построение глобальных огибающих для каждого шаблона и использование правила продукта для многократного тестирования. Предположим, что в данных есть M точечных паттернов, и вам нужен тест (подобранной модели) уровня значимости альфа (обычно альфа = 0,05). Затем вы собираетесь построить M конвертов, по одному для каждого шаблона, каждый с уровнем значимости gamma = 1 - (1-alpha)^(1/M)
, Каждый конверт будет создан envelope.ppp с global=TRUE
а также nsim = 1/gamma - 1
примерно. Пример: если М = 10 и альфа = 0,05, то gamma=1 - 0.95^(1/10) = 0.0051
так nsim=1/gamma-1 = 194.45
, назови это nsim=195
, Поскольку вы собираетесь делать глобальные конверты, вам может потребоваться вдвое больше симуляций, как объясняется в справке для envelope
, Так что сделайте следующее, где fit
это подогнанная модель и H
исходный гиперкадр данных:
sims <- simulate(fit, nsim=2*195)
SIMS <- list()
for(i in 1:nrow(sims)) SIMS[[i]] <- as.solist(sims[i,,drop=TRUE])
Hplus <- cbind(H, hyperframe(Sims=SIMS))
Это увеличивает исходный гиперкадр, добавляя столбец имитируемых шаблонов (каждая запись в столбце является "солистом", содержащим 2*195 шаблонов). Тогда делай (где X
это столбец H
содержащий исходные наборы данных точечных рисунков)
EE <- with(Hplus, envelope(X, Lest, global=TRUE, nsim=195, simulate=Sims))
plot(EE)
Это дает сюжет с множеством панелей конвертов. Их интерпретация заключается в том, что если какая-либо из наблюдаемых L-функций выходит за пределы соответствующих огибающих, результат является значительным - тест отвергает нулевую гипотезу о том, что подобранная модель верна.
Первая ошибка - это ошибка, исправленная несколько дней назад в версии для разработчиков. spatstat
, Если у вас есть devtools
пакет, который вы можете получить последний (и лучший) spatstat
без труда:
devtools::install_github('spatstat/spatstat')
Дайте мне знать, если это не вариант для вас (а также на какой платформе вы находитесь).
Вторая ошибка действительно потому, что envelope
не реализовано для класса mppm
Таким образом, вы должны разработать обходной путь на данный момент, как вы говорите.
Я думаю, что с тем, что вы уже сделали, есть несколько проблем: модель неоднородна, но вы используете Lcross
скорее, чем Lcross.inhom
и я думаю, что ваш single.E
эквивалентно (где вы не используете подобранную модель вообще):
single.E <- envelope(lansing, Lcross, i="blackoak", j="maple", nsim=39, rank=1, global=TRUE, correction="best", simulate=expression(rlabel(lansing)))
Дайте мне знать, как вы прогрессируете с этим. Я мог бы найти время, чтобы дать более подробную информацию об обходе пропавших без вести envelope.mppm
(объединение итоговых функций для каждого шаблона).