Определение кластеризации деревьев в лесных промежутках с использованием R?

14

Прилагаемый набор данных показывает приблизительно 6000 саженцев в приблизительно 50 лесных промежутках переменного размера. Мне интересно узнать, как эти саженцы растут в пределах их соответствующих пробелов (то есть кластеризованных, случайных, рассеянных). Как вы знаете, традиционный подход состоял бы в том, чтобы запустить Global Moran's I. Однако скопления деревьев в скоплениях промежутков, по-видимому, являются нецелесообразным использованием Moran I. Я провел несколько тестовых статистик с Moran's I, используя пороговое расстояние 50 метров, который дал бессмысленные результаты (т. е. р-значение = 0,0000000 ...). Взаимодействие между скоплениями пробелов, вероятно, дает эти результаты. Я рассмотрел вопрос создания сценария для прохождения отдельных пробелов купола и определения кластеризации внутри каждого пробела, хотя отображение этих результатов для общественности было бы проблематичным.

Каков наилучший подход для количественной оценки кластеризации внутри кластеров?

введите описание изображения здесь

Аарон
источник
1
Аарон, вы говорите, что пытались запустить Морана I, вам интересно измерить, как атрибут саженца сравнивается с атрибутами соседних саженцев (то есть вы имеете дело с отмеченным точечным образцом)? Название, по-видимому, подразумевает, что вас интересует только расположение саженцев относительно друг друга, а не их атрибуты.
MannyG
@MannyG Да, меня интересует только определение того, сгруппированы ли саженцы относительно расположения других саженцев в пределах какого-либо данного лесного разрыва. Существует только один вид, представляющий интерес, а размер саженцев не представляет интереса.
Аарон

Ответы:

7

У вас нет однородного случайного поля, поэтому попытка проанализировать все ваши данные сразу нарушит допущения любой статистики, которую вы выберете для решения проблемы. Из вашего поста неясно, являются ли ваши данные отмеченным точечным процессом (т. Е. Диаметром или высотой, связанной с каждым местоположением дерева). Если эти данные не представляют отмеченный точечный процесс, я понятия не имею, как вы применили Моран-I. Если данные представляют только пространственные местоположения, я бы рекомендовал использовать Ripley's-K с преобразованием Безаг-L, чтобы стандартизировать нулевое ожидание в нуле. Это позволяет проводить оценку кластеризации в нескольких масштабах. Если ваши данные имеют ассоциированное значение, то ваш лучший вариант - местный Моран-I (LISA). Я бы на самом деле посмотрел на это с обеих статистик. Независимо от вашего выбора, Вам все равно нужно будет пройтись по каждому отдельному сайту, чтобы получить достоверные результаты. Вот пример кода R для симуляции Монте-Карло Ripley's-K / Besag's-L с использованием встроенного набора данных саженцев из красного дерева. Это должно быть довольно просто изменить, чтобы пройтись по вашим сайтам и создать график для каждого из них.

# ADD REQUIRED PACKAGES
require(sp)
require(spatstat)
options(scipen=5)

# USE REDWOOD SAPLING DATASET
spp <- SpatialPoints(coords(redwood))

###################################################
###### START BESAG'S-L MONTE CARLO  ANALYSUS ######
###################################################
# CREATE CONVEX HULL FOR ANALYSIS WINDOW                       
W=ripras(coordinates(spp)) 

# COERCE TO spatstat ppp OBJECT
spp.ppp=as.ppp(coordinates(spp), W)                     
  plot(spp.ppp) 

# ESTIMATE BANDWIDTH
area <- area.owin(W)
lambda <- spp.ppp$n/area
 ripley <- min(diff(W$xrange), diff(W$yrange))/4
   rlarge <- sqrt(1000/(pi * lambda))
     rmax <- min(rlarge, ripley)
bw <- seq(0, rmax, by=rmax/10)  

# CALCULATE PERMUTED CROSS-K AND PLOT RESULTS       
Lenv <- envelope(spp.ppp, fun="Kest", r=bw, i="1", j="2", nsim=99, nrank=5, 
                 transform=expression(sqrt(./pi)-bw), global=TRUE)            
plot(Lenv, main="Besag's-L", xlab="Distance", ylab="L(r)", legend=F, col=c("white","black","grey","grey"), 
    lty=c(1,2,2,2), lwd=c(2,1,1,1) )
     polygon( c(Lenv$r, rev(Lenv$r)), c(Lenv$lo, rev(Lenv$hi)), col="lightgrey", border="grey")
       lines(supsmu(bw, Lenv$obs), lwd=2)
       lines(bw, Lenv$theo, lwd=1, lty=2)
         legend("topleft", c(expression(hat(L)(r)), "Simulation Envelope", "theo"), pch=c(-32,22),
                col=c("black","grey"), lty=c(1,0,2), lwd=c(2,0,2), pt.bg=c("white","grey"))
Джеффри Эванс
источник
1
Но вы не можете просто использовать выпуклый корпус в качестве окна для вашего точечного рисунка! Помните, что окно - это область, в которой работает шаблон, который создает точки. Вы априори знаете, что деревья растут только в этих заданных областях, и вам нужно настроить окно, чтобы это отражать. Вы можете смягчить это, установив диапазон K (r) на очень маленький размер порядка 0,3x ваших клирингов, но вы получите предвзятые результаты из-за отсутствия коррекции краевого эффекта. Джеффри использует размер всей области исследования, чтобы определить его rmax.
Spacedman
1
В моем примере, да, я использую весь регион. Именно поэтому я рекомендовал проходить через каждый образец сайта (пробел). Каждый раз, когда вы устанавливаете под определенную область выборки, вы запускаете анализ заново. Вы не можете рассматривать всю область исследования как свое случайное поле, потому что у вас нет непрерывной выборки. Имея только пробные пробелы, вы фактически получаете независимые графики. Функция Kest, которую я вызываю, по умолчанию использует коррекцию края "border". Доступны и другие варианты коррекции края. Я бы сказал, что ваша экспериментальная установка - это пробел в куполе, и ее следует проанализировать как таковую.
Джеффри Эванс
1
Думая об этом немного больше. Вы действительно должны использовать полигоны, которые представляют каждый пробел в качестве вашего окна. Если вы зададите свою задачу для отражения экспериментальной единицы, тогда CSR и K будут смещены, потому что площадь не отражает фактический размер зазора купола. Это проблема в рекомендациях my и @ Spacedman.
Джеффри Эванс
2
Обратите внимание, что в моем расширенном примере использовалась только грубая сетка, потому что это был довольно простой способ создания чего-либо с примерно правильной структурой. Ваша маска должна выглядеть как карта ваших открытых лесных массивов. Технически неправильно пытаться определить маску из данных!
Spacedman
1
@Spacedman. Мне нравится ваш подход, и он, безусловно, эффективен. Моя конкретная проблема заключается в том, что разрывы навеса являются экспериментальной единицей. В вашем подходе, если два промежутка проксимальны, полоса пропускания может правдоподобно включать наблюдения из разных единиц выборки. Кроме того, результирующая статистика не должна отражать «пул» экспериментальных единиц, но должна быть репрезентативной для каждой единицы и выводить о пространственном процессе, извлеченном из общих моделей по экспериментальным единицам. Если рассматривать глобально, это представляет нестационарный процесс интенсивности, который нарушает статистические предположения.
Джеффри Эванс
4

То, что у вас есть, это точечный шаблон с окном, которое представляет собой множество небольших разрозненных многоугольных областей.

Вы должны быть в состоянии использовать любой из тестов package:spatstatдля CSR, если вы вводите его с правильным окном. Это может быть либо количество наборов (x, y) пар, определяющих каждую очистку, либо двоичная матрица значений (0,1) в пространстве.

Сначала давайте определим что-то похожее на ваши данные:

set.seed(310366)
nclust <- function(x0, y0, radius, n) {
               return(runifdisc(n, radius, centre=c(x0, y0)))
             }
c = rPoissonCluster(15, 0.04, nclust, radius=0.02, n=5)
plot(c)

и давайте притворимся, что наши расчеты - это квадратные ячейки, которые просто таковы:

m = matrix(0,20,20)
m[1+20*cbind(c$x,c$y)]=1
imask = owin(c(0,1),c(0,1),mask = t(m)==1 )
pp1 = ppp(x=c$x,y=c$y,window=imask)
plot(pp1)

Таким образом, мы можем построить K-функцию этих точек в этом окне. Мы ожидаем, что это будет не CSR, потому что точки кажутся сгруппированными внутри ячеек. Обратите внимание, что я должен изменить диапазон расстояний, чтобы он был небольшим - порядка размера ячейки - в противном случае K-функция оценивается на расстоянии, равном размеру всего шаблона.

plot(Kest(pp1,r=seq(0,.02,len=20)))

Если мы сгенерируем несколько точек CSR в тех же ячейках, мы сможем сравнить графики K-функций. Этот должен быть больше похож на CSR:

ppSim = rpoispp(73/(24/400),win=imask)
plot(ppSim)
plot(Kest(ppSim,r=seq(0,.02,len=20)))

двухточечные узоры в пятнистых окнах

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

Spacedman
источник
2

В дополнение к посту Энди:

То , что вы хотите вычислить является мерой пространственной однородности (эрго гипотезы: «группируются ваши очки?») , Такие как L Рипли и функции K .

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

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