Проверка гипотез на основе моделирования на гиперкадрах пространственных точечных шаблонов с использованием функции огибающей в 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» распознает «mark» как столбец с названиями видов, а «sim1» как ковариату.

Затем выполните проверку гипотез на основе моделирования, где вас интересует, находятся ли Black Oak и Maple в одних и тех же местах.

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)

Выдается следующее сообщение об ошибке:

Ошибка в контрольных переменных (формула, data.sumry $ col.names, extra = c ("x", "y",: переменная "mark" в формуле не является одним из имен в данных "

Я получаю ту же ошибку, используя:

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 ("envelope"): нет применимого метода для "envelope", примененного к объекту класса "c ('mppm', 'list')"

Подразумевая, что "конверт" не работает с объектами mppm (только ppp или ppm). Я подозреваю, что есть способ обойти это ограничение с умом, но я его еще не нашел. Любые предложения или рекомендации были бы очень полезны!


person A. Brown    schedule 02.07.2016    source источник


Ответы (2)


В spatstat нет envelope.mppm, потому что некоторые статистические проблемы (связанные с проверкой множественных гипотез) еще не решены.

Самым быстрым решением, вероятно, является использование cdf.test.mppm, который выполнит тест и выдаст некоторый графический результат.

Было бы разумно объединить огибающие (например, функций K) из разных точечных шаблонов для получения единого конверта, при условии, что подобранная модель подразумевает, что разные шаблоны должны быть статистически эквивалентным. Это было бы неверно, если, например, модель включает ковариаты, которые различны для разных паттернов.

Лучшей стратегией, вероятно, является построение глобальных огибающих для каждого шаблона и использование правила продукта для множественного тестирования. Предположим, что в данных есть M точечных паттернов, и вам нужен тест (подобранной модели) уровня значимости альфа (обычно альфа = 0,05). Затем вы собираетесь построить M конвертов, по одному для каждого шаблона, каждый с уровнем значимости gamma = 1 - (1-alpha)^(1/M). Каждый конверт будет сгенерирован envelope.ppp примерно с global=TRUE и nsim = 1/gamma - 1. Пример: если M = 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-функций выходит за пределы соответствующих огибающих, результат является значительным - тест отклоняет нулевую гипотезу о том, что подобранная модель верна.

person Adrian Baddeley    schedule 04.07.2016
comment
Ваш код потрясающе работает, @Adrian Baddeley, спасибо! Если бы я хотел проверить гипотезу случайной метки в гиперкадре, я предполагаю, что мог бы использовать аналогичный подход: EE <- with(H, envelope.ppp(spp, Lest, global=TRUE, nsim=195, simulate=expression(rlabel(spp))) Где я использую шаблон необработанных точек (spp) для имитации случайных меток на шаблоне и правильное количество имитаций для приходиться на многократное тестирование? - person A. Brown; 08.07.2016

Первая ошибка - это ошибка, которая была исправлена ​​несколько дней назад в разрабатываемой версии 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 (объединение сводных функций для каждого шаблона).

person Ege Rubak    schedule 02.07.2016
comment
Спасибо, доктор Рубак! Установка новой версии spatstat устранила первую проблему. Я продолжу мозговой штурм по второй проблеме; хотя моя новизна в программировании означает, что это может быть очень далеко. - person A. Brown; 03.07.2016
comment
Очень простое решение проблемы № 2 (выполнение проверки гипотез с помощью конвертов моделирования для нескольких точечных шаблонов), по-видимому, состоит в том, чтобы взять среднее (или взвешенное среднее, взвешенное по размеру выборки) obs, mmean, hi и lo (каждый из которых является выходом из функционального конверта) и создайте график вручную, используя эти числа. Это слишком упрощенное решение? - person A. Brown; 03.07.2016