Menentukan apakah pohon dalam celah hutan dikelompokkan menggunakan R?

14

Dataset terlampir menunjukkan sekitar 6000 pohon muda di sekitar 50 celah hutan berukuran variabel. Saya tertarik mempelajari bagaimana anakan ini tumbuh di dalam celahnya masing-masing (mis. Bergerombol, acak, tersebar). Seperti yang Anda ketahui, pendekatan tradisional adalah menjalankan Global Moran I. Namun, agregasi pohon dalam agregasi kesenjangan tampaknya merupakan penggunaan yang tidak tepat dari Moran I. Saya menjalankan beberapa statistik uji dengan Moran I menggunakan jarak ambang 50 meter, yang menghasilkan hasil yang tidak masuk akal (yaitu p-value = 0,0000000 ...). Interaksi di antara agregasi gap kemungkinan menghasilkan hasil ini. Saya telah mempertimbangkan membuat skrip untuk loop melalui celah kanopi individu dan menentukan pengelompokan dalam setiap celah, meskipun menampilkan hasil ini kepada publik akan bermasalah.

Apa pendekatan terbaik untuk mengukur pengelompokan dalam kelompok?

masukkan deskripsi gambar di sini

Harun
sumber
1
Aaron, Anda mengatakan bahwa Anda mencoba menjalankan Moran's I, apakah Anda tertarik untuk mengukur bagaimana atribut anak pohon dibandingkan dengan atribut anak pohon tetangga (yaitu, apakah Anda berurusan dengan pola titik yang ditandai )? Judul tersebut sepertinya menyiratkan bahwa Anda hanya tertarik pada lokasi anakan relatif satu sama lain dan bukan atributnya.
MannyG
@MannyG Ya, saya hanya tertarik untuk menentukan apakah anakan dikelompokkan relatif terhadap lokasi anakan lainnya di dalam celah hutan tertentu. Hanya ada satu spesies yang menarik dan ukuran anakan tidak menarik.
Aaron

Jawaban:

7

Anda tidak memiliki bidang acak seragam, jadi mencoba menganalisis semua data Anda sekaligus akan melanggar asumsi statistik apa pun yang Anda pilih untuk dilemparkan pada masalah. Tidak jelas dari pos Anda apakah data Anda merupakan titik proses yang ditandai (yaitu, diameter atau tinggi yang terkait dengan setiap lokasi pohon). Jika data ini tidak mewakili proses titik yang ditandai, saya tidak tahu bagaimana Anda menerapkan Moran's-I. Jika data hanya mewakili lokasi spasial, saya akan merekomendasikan menggunakan Ripley's-K dengan transformasi Besag-L untuk membakukan harapan nol di nol. Hal ini memungkinkan untuk penilaian multiskala pengelompokan. Jika data Anda memiliki nilai terkait, maka opsi terbaik Anda adalah Moran's-I (LISA) setempat. Saya benar-benar akan melihatnya dengan kedua statistik. Apa pun pilihan Anda, Anda masih perlu mengulang-ulang setiap situs untuk menghasilkan hasil yang valid. Berikut adalah beberapa contoh kode R untuk simulasi Monte Carlo dari Ripley's-K / Besag's-L menggunakan dataset sapling redwood bawaan. Seharusnya cukup mudah untuk memodifikasi ini untuk mengulang melalui situs Anda dan menghasilkan grafik untuk masing-masing situs.

# 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
sumber
1
Tetapi Anda tidak bisa begitu saja menggunakan cembung cembung sebagai jendela untuk pola titik Anda! Ingat, jendela adalah area di mana pola yang menghasilkan titik beroperasi. Anda tahu a-apriori bahwa pohon-pohon hanya tumbuh di daerah-daerah tertentu ini, dan Anda harus mengatur jendela Anda untuk mencerminkannya. Anda mungkin mengurangi ini dengan mengatur rentang K (r) untuk sesuatu yang sangat kecil, dari urutan ukuran 0,3x dari pembukaan Anda, tetapi Anda akan mendapatkan hasil yang bias karena kurangnya koreksi efek tepi. Jeffrey menggunakan ukuran seluruh area studi untuk mendefinisikan rmax-nya.
Spacedman
1
Dalam contoh saya, Ya saya menggunakan seluruh wilayah. Dengan ini persis mengapa saya merekomendasikan perulangan melalui setiap situs sampel (celah). Setiap kali Anda subset ke area sampel tertentu Anda akan menjalankan kembali analisis. Anda tidak dapat memperlakukan seluruh area studi sebagai bidang acak Anda karena Anda tidak memiliki pengambilan sampel berkelanjutan. Akibatnya, hanya memiliki celah sampel Anda, memiliki plot independen. Fungsi Kest yang saya panggil, secara default, menggunakan koreksi tepi "perbatasan". Ada opsi koreksi tepi lainnya yang tersedia. Saya berpendapat bahwa unit eksperimental Anda adalah celah kanopi dan harus dianalisis.
Jeffrey Evans
1
Dalam memikirkan hal ini sedikit lebih. Anda harus benar-benar menggunakan poligon yang mewakili setiap celah sebagai jendela Anda. Jika Anda mengelompokkan masalah Anda untuk mencerminkan unit eksperimental maka CSR dan K akan menjadi bias karena area tersebut tidak mencerminkan ukuran celah kanopi yang sebenarnya. Ini adalah masalah dalam rekomendasi saya dan @ Spacedman.
Jeffrey Evans
2
Perhatikan contoh saya yang diperluas hanya menggunakan grid kasar karena itu adalah cara yang cukup sederhana untuk membuat sesuatu dengan struktur yang kira-kira tepat. Topeng Anda akan terlihat seperti peta wilayah hutan terbuka Anda. Secara teknis salah untuk mencoba dan mendefinisikan topeng dari data!
Spacedman
1
@Spacedman. Saya suka pendekatan Anda dan ini tentu saja efisien. Kekhawatiran khusus saya adalah bahwa celah kanopi adalah unit eksperimental. Dalam pendekatan Anda, jika dua celah adalah proksimal, bandwidth dapat masuk akal termasuk pengamatan dari unit pengambilan sampel yang berbeda. Selain itu, statistik yang dihasilkan tidak boleh mencerminkan "kumpulan" unit eksperimental tetapi harus mewakili masing-masing unit dan menarik kesimpulan tentang proses spasial yang diambil dari pola umum di seluruh unit eksperimental. Jika diperlakukan secara global, ini merupakan proses intensitas non-stasioner yang, melanggar asumsi statistik.
Jeffrey Evans
4

Apa yang Anda miliki adalah pola titik dengan jendela yang merupakan sejumlah kecil wilayah poligon yang terputus.

Anda harus dapat menggunakan salah satu tes package:spatstatuntuk CSR selama Anda mengisinya dengan jendela yang benar. Ini bisa berupa sejumlah set (x, y) pasangan yang mendefinisikan setiap kliring atau matriks biner dari nilai (0,1) di atas spasi.

Pertama mari kita mendefinisikan sesuatu yang sedikit mirip dengan data Anda:

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)

dan mari kita berpura-pura membersihkan kita adalah sel persegi yang kebetulan seperti ini:

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)

Jadi kita bisa memplot fungsi-K dari titik-titik di jendela itu. Kami berharap ini menjadi non-CSR karena poin-poinnya tampak berkerumun di dalam sel. Perhatikan saya harus mengubah kisaran jarak menjadi kecil - dari urutan ukuran sel - jika fungsi-K akan dievaluasi jarak jarak ukuran pola keseluruhan.

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

Jika kita menghasilkan beberapa titik CSR di sel yang sama, kita dapat membandingkan plot fungsi-K. Yang ini harus lebih seperti CSR:

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

dua pola titik di windows tambal sulam

Anda tidak dapat benar-benar melihat titik-titik yang dikelompokkan dalam sel dalam pola pertama, tetapi jika Anda memplotnya sendiri di jendela grafis itu jelas. Titik-titik dalam pola kedua seragam di dalam sel (dan tidak ada di wilayah hitam) dan fungsi-K jelas berbeda dari Kpois(r), fungsi K CSR untuk data yang dikelompokkan dan serupa untuk data yang seragam.

Spacedman
sumber
2

Selain posting Andy:

Yang ingin Anda hitung adalah ukuran homogenitas spasial (sesuai hipotesis: "Apakah poin Anda terkelompok?") Seperti fungsi L dan K Ripley .

Posting blog ini menjelaskan caranya di R dengan cukup baik. Berdasarkan kode yang dijelaskan, saya pertama-tama akan memberi label pada masing-masing cluster dalam dataset Anda dan kemudian menghitung dalam satu lingkaran untuk setiap cluster amplop kritis melalui Ripley's K

Curlew
sumber
Saat ini saya menghapus jawaban saya. Beberapa analisis singkat menunjukkan bahwa secara oportunis mengidentifikasi plot berdasarkan K-berarti bias statistik lokal lebih terkelompok daripada yang disarankan secara kebetulan. Jawaban ini masih berlaku meskipun +1, (hanya pembuatan windows dari data lebih bermasalah daripada jawaban asli saya akan menyarankan).
Andy W