Ustalanie, czy drzewa w obrębie luk leśnych są grupowane za pomocą R?

14

Załączony zestaw danych pokazuje około 6000 sadzonek w około 50 lukach leśnych o zmiennej wielkości. Interesuje mnie nauka, w jaki sposób te sadzonki rosną w obrębie swoich odpowiednich luk (tj. Skupione, losowe, rozproszone). Jak wiecie, tradycyjnym podejściem byłoby uruchomienie Globalnego Morana I. Jednak agregacje drzew w obrębie agregatów luk wydają się być niewłaściwym użyciem Morana I. Przeprowadziłem pewne statystyki testowe z Moranem I, stosując odległość progową 50 metrów, co dało bezsensowne wyniki (tj. wartość p = 0,0000000 ...). Interakcja między agregacjami luk prawdopodobnie przyniesie te wyniki. Zastanawiałem się nad stworzeniem skryptu do przechodzenia przez poszczególne luki w czaszy i określania klastrowania w obrębie każdej luki, chociaż wyświetlanie tych wyników opinii publicznej byłoby problematyczne.

Jakie jest najlepsze podejście do kwantyfikacji klastrów w klastrach?

wprowadź opis zdjęcia tutaj

Aaron
źródło
1
Aaron, mówisz, że próbowałeś uruchomić Morana I, czy chcesz zmierzyć, jak atrybut drzewka porównuje się z atrybutami sąsiednich drzewek (tj. Czy masz do czynienia z zaznaczonym wzorem punktowym)? Tytuł wydaje się sugerować, że interesuje Cię tylko lokalizacja sadzonek względem siebie, a nie ich atrybuty.
MannyG
@MannyG Tak, interesuje mnie tylko ustalenie, czy drzewka są skupione w stosunku do lokalizacji innych drzewek w obrębie danej luki leśnej. Istnieje tylko jeden interesujący gatunek, a wielkość sadzonek nie jest interesująca.
Aaron

Odpowiedzi:

7

Nie masz jednolitego losowego pola, więc próba analizy wszystkich twoich danych naraz naruszy założenia każdej statystyki, którą wybierzesz na problem. Z twojego postu nie jest jasne, czy dane są procesem oznaczonym punktem (tj. Średnicą lub wysokością związaną z każdą lokalizacją drzewa). Jeśli te dane nie reprezentują procesu oznaczonego punktu, nie mam pojęcia, w jaki sposób zastosowałeś Morana. Jeśli dane reprezentują tylko lokalizacje przestrzenne, zaleciłbym użycie Ripley's-K z transformacją Besag-L, aby ustandaryzować zerowe oczekiwanie na zero. Pozwala to na wieloskalową ocenę klastrowania. Jeśli dane mają powiązaną wartość, najlepszą opcją jest lokalny Moran's-I (LISA). Spojrzałbym na to z obu statystyk. Bez względu na twój wybór nadal będziesz musiał przeglądać poszczególne witryny, aby uzyskać prawidłowe wyniki. Oto przykładowy kod R dla symulacji Monte Carlo dla Ripley's-K / Besag's-L przy użyciu wbudowanego zestawu danych drzewka sekwoi. Modyfikacja tego powinna być dość prosta, aby przeglądać witryny i tworzyć wykresy dla każdej z nich.

# 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"))
Jeffrey Evans
źródło
1
Ale nie możesz po prostu użyć wypukłego kadłuba jako okna dla swojego wzoru punktowego! Pamiętaj, że okno to obszar, w którym działa wzór tworzący punkty. Wiesz z góry, że drzewa rosną tylko w tych określonych regionach, i musisz ustawić okno, aby to odzwierciedlić. Możesz to złagodzić, ustawiając zakres K (r) na coś bardzo małego, rzędu wielkości 0,3x twoich prześwitów, ale uzyskasz stronnicze wyniki z powodu braku korekt efektu krawędzi. Jeffrey używa wielkości całego obszaru badań do zdefiniowania swojego rmax.
Spacedman
1
W moim przykładzie: Tak, używam całego regionu. Właśnie dlatego zaleciłem zapętlanie każdej przykładowej strony (luki). Za każdym razem, gdy wybierzesz określony obszar próbki, ponownie uruchomisz analizę. Nie możesz traktować całego obszaru badań jako pola losowego, ponieważ nie masz ciągłego próbkowania. Mając tylko próbkowane luki, w rzeczywistości masz niezależne wykresy. Wywołana przeze mnie funkcja Kest domyślnie wykorzystuje korekcję krawędzi „granicy”. Dostępne są inne opcje korekcji krawędzi. Twierdziłbym, że twoją jednostką eksperymentalną jest luka baldachimowa i powinna być jako taka analizowana.
Jeffrey Evans
1
Myśląc o tym trochę więcej. Naprawdę powinieneś używać wielokątów, które reprezentują każdą przerwę jako okno. Jeśli poddziałsz swój problem, aby odzwierciedlić jednostkę eksperymentalną, wówczas CSR i K będą tendencyjne, ponieważ obszar nie odzwierciedla faktycznie rozmiaru szczeliny czaszy. Jest to problem zarówno w moich zaleceniach, jak i @ Spacedman.
Jeffrey Evans
2
Zauważ, że mój rozszerzony przykład użyłem tylko grubej siatki, ponieważ był to dość prosty sposób na stworzenie czegoś z mniej więcej odpowiednią strukturą. Twoja maska ​​powinna wyglądać jak mapa otwartych obszarów leśnych. Z technicznego punktu widzenia próba zdefiniowania maski na podstawie danych jest błędna!
Spacedman
1
@Spacedman. Podoba mi się twoje podejście i na pewno jest skuteczne. Moją szczególną troską jest to, że luki czaszy są jednostką eksperymentalną. W twoim podejściu, jeśli dwie przerwy są bliższe, szerokość pasma może prawdopodobnie obejmować obserwacje z różnych jednostek próbkowania. Ponadto wynikowa statystyka nie powinna odzwierciedlać „puli” jednostek eksperymentalnych, ale powinna być reprezentatywna dla każdej jednostki i wnioskować na temat procesu przestrzennego czerpanego ze wspólnych wzorców w jednostkach eksperymentalnych. Jeśli traktowane globalnie, reprezentuje niestacjonarny proces intensywności, który narusza założenia statystyczne.
Jeffrey Evans
4

To, co masz, to wzór punktowy z oknem, które jest liczbą małych rozłączonych regionów wielokąta.

Powinieneś być w stanie skorzystać z dowolnego z testów package:spatstatCSR, o ile podasz je odpowiednim oknem. Może to być liczba zestawów par (x, y) definiujących każde czyszczenie lub macierz binarna wartości (0,1) w przestrzeni.

Najpierw zdefiniujmy coś, co wygląda trochę jak twoje dane:

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)

i udawajmy, że nasze polany to kwadratowe komórki, które akurat są takie:

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)

Możemy więc wykreślić funkcję K tych punktów w tym oknie. Oczekujemy, że nie będzie to CSR, ponieważ punkty wydają się skupione w komórkach. Zauważ, że muszę zmienić zakres odległości, aby był mały - rzędu wielkości komórki - w przeciwnym razie funkcja K będzie oceniana na odległościach wielkości całego wzoru.

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

Jeśli wygenerujemy punkty CSR w tych samych komórkach, możemy porównać wykresy funkcji K. Ten powinien być bardziej podobny do CSR:

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

dwa wzorce punktowe w niejednolitych oknach

Tak naprawdę nie widać punktów skupionych w komórkach w pierwszym wzorze, ale jeśli narysujesz go samodzielnie w oknie graficznym, będzie to jasne. Punkty w drugim wzorze są jednolite w komórkach (i nie istnieją w czarnym obszarze), a funkcja K jest wyraźnie różna od Kpois(r)funkcji K CSR dla danych klastrowych i podobna dla danych jednolitych.

Spacedman
źródło
2

Oprócz postu Andy'ego:

To, co chcesz obliczyć, to miara jednorodności przestrzennej (ergo hipoteza: „Czy twoje punkty są skupione?”), Takie jak funkcja L i K Ripleya .

Ten post na blogu całkiem dobrze wyjaśnia, jak to zrobić w R. W oparciu o opisany kod najpierw oznaczyłbym każdy klaster w zbiorze danych, a następnie obliczyłem w każdej pętli krytyczną obwiednię za pomocą K Ripleya

Kulik
źródło
Obecnie usunąłem swoją odpowiedź. Pewna krótka analiza sugeruje, że oportunistyczna identyfikacja wykresów w oparciu o K oznacza, że ​​lokalne statystyki są bardziej skupione niż sugerowałoby to przypadek. Ta odpowiedź wciąż obowiązuje, ale +1 (samo tworzenie okien z danych jest bardziej problematyczne niż sugerowałaby moja pierwotna odpowiedź).
Andy W