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

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?

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"))

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)))

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.

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