Đây là bài thực hành về chủ đề xác định kiểu phân bố của các cá thể cây rừng trong không gian với phần mềm R.
Tiêu chuẩn Clark-Evans
Có 3 kiểu phân bố của cây rừng trong không gian: ngẫu nhiên, cụm hoặc đều.
Nguồn: https://docs.qgis.org
Để xác định kiểu phân bố của đối tượng nghiên cứu, chúng ta có thể sử dụng tiêu chuẩn Clark-Evans được phát triển bởi 2 nhà khoa học Philip Clark và Francis Evans năm 1954 (Clark, P.J. and Evans, F.C. (1954) Distance to nearest neighbour as a measure of spatial relationships in populations Ecology 35, 445–453.)
Tiêu chuẩn này được tính dựa vào tỉ lệ giữa giá trị trung bình quan sát được \(\bar{r}_{A}\) và giá trị trung bình mong đợi \(\bar{r}_{E}\) về khoảng cách của các cá thể nghiên cứu đến cá thể gần nhất của nó:
\[R=\frac{\bar{r}_{A}}{\bar{r}_{E}}\]
Để tính \(\bar{r}_{A}\), ta chỉ cần đo khoảng cách của từng cá thể đến cá thể gần nhất của nó trong ô mẫu nghiên cứu. Sau đó lấy giá trị trung bình của các khoảng cách đã đo được.
Còn \(\bar{r}_{E}\) được tính bằng mật độ cá thể theo công thức sau:
\[\bar{r}_{E}=\frac{1}{2\sqrt{\lambda}}=\frac{1}{2\sqrt{\frac{n}{S}}}\]
Với \(\lambda = \frac{n}{S}\)
- n: Số cá thể
- S: Diện tích khu vực nghiên cứu
Khi đó, nếu chỉ số R này
- R = 1 thì các cá thể được xem là phân bố ngẫu nhiên trong khu vực nghiên cứu
- R > 1 phân bố đều
- R < 1 phân bố cụm
Việc so sánh này cần được thực hiện bằng một kiểm định thống kê.
Thực hành
Mở tập tin bản đồ
Ví dụ chúng ta có 2 tập tin vector về khu vực ô mẫu nghiên cứu và các cá thể quan sát của loài trong ô mẫu đó. Chúng ta sử dụng readOGR() để mở các tập tin này trong R.
library(maptools)
library(sp)
library(rgdal)
myplot <- readOGR('https://raw.githubusercontent.com/nlxbach/data/main/myplot.geojson')## OGR data source with driver: GeoJSON
## Source: "https://raw.githubusercontent.com/nlxbach/data/main/myplot.geojson", layer: "myplot"
## with 1 features
## It has 1 fieldsmyspecies <- readOGR('https://raw.githubusercontent.com/nlxbach/data/main/myspecies.geojson')## OGR data source with driver: GeoJSON
## Source: "https://raw.githubusercontent.com/nlxbach/data/main/myspecies.geojson", layer: "myspecies2"
## with 100 features
## It has 2 fieldsĐây là phân bố các cá thể của loài nghiên cứu
plot(myplot)
points(myspecies, pch = 19)
Tính toán các giá trị
Ta sử dụng gói spatstat cho các tính toán này. Do đó ta cần chuyển dữ liệu về đúng định dạng Point Pattern để có thể sử dụng các lệnh của gói này.
library(spatstat)
# ô mẫu
myplotOwin <- as.owin(W = myplot)
# Loài
myspecies_ppp <- ppp(x = myspecies@coords[,1],
y = myspecies@coords[,2],
window = myplotOwin)Thực hiện các phép tính thủ công:
Tính giá trị \(\bar{r}_{A}\):
Đầu tiên ta cần tính khoảng cách từng cá thể đến cá thể khác gần nhất của nó. Có nghĩa là nếu trên thực địa, tại mỗi cây ta cần xác định 1 cây khác mà gần nó nhất rồi đo khoảng cách giữa 2 cây này. Nếu khảo sát trên 100 cá thể thì ta sẽ có 100 giá trị khoảng cách, sau đó ta tính giá trị trung bình từ 100 khoảng cách này. Trong R ta có thể tính với hàm nndist()
(nndistX <- nndist(myspecies_ppp))## [1] 1592.2969 656.6641 577.8883 1701.2331 2414.5308 2481.3491 1106.0310
## [8] 1641.0137 1254.0177 535.2841 1395.4634 1238.8146 1544.9810 722.8317
## [15] 1169.3213 517.0831 724.0788 786.1549 2195.5257 2499.4656 2534.8725
## [22] 735.5827 1771.0017 839.6179 1335.5047 577.8883 656.7571 3666.9263
## [29] 1871.7707 1305.3942 656.6641 1106.0310 1768.7605 2827.9720 1927.8402
## [36] 425.5805 722.8317 1335.5047 253.2553 535.2841 921.5949 1373.5472
## [43] 380.7128 1223.9652 1700.8107 1995.7779 1768.7605 985.6549 1373.5472
## [50] 1422.9000 1418.3159 253.2553 2268.9642 1600.9863 1039.6989 425.5805
## [57] 1223.9652 1425.0644 344.4207 2174.1745 1619.8098 634.7580 344.4207
## [64] 2087.1714 839.6179 1109.4159 924.6595 2174.1745 1927.8402 634.7580
## [71] 735.5827 1085.9412 1723.1577 1253.3906 924.6595 2324.4844 1037.9505
## [78] 985.6549 126.0070 126.0070 1700.8107 1418.3159 1913.8332 1297.7528
## [85] 1235.2553 2127.4401 1152.3458 517.0831 1995.7779 786.1549 921.5949
## [92] 1039.6989 724.0788 1189.6770 1091.6188 1297.7528 2003.6321 1139.5312
## [99] 380.7128 2889.8414\(\bar{r}_{A}\) là giá trị trung bình của các khoảng cách trên: 1283.731 mét
(Dobs <- mean(nndistX))## [1] 1283.731Tính giá trị \(\bar{r}_{E}\):
Diện tích khu vực nghiên cứu S: 681 028 249 \(m^{2}\)
(areaW <- area(myplotOwin))## [1] 681028249Số lượng cá thể điều tra: 100 cá thể
(npts <- npoints(myspecies_ppp))## [1] 100Khi đó áp dụng vào công thức tính \(\bar{r}_{E}\) ta được giá trị: 1304.826 mét
intensity <- npts/areaW
(Dpois <- 1/(2 * sqrt(intensity)))## [1] 1304.826Tỉ lệ Clark-Evans: 0,98
(Dobs/Dpois)## [1] 0.9838335Chúng ta vừa thực hiện các bước thủ công để tính tỉ lệ Clark-Evans. Tuy nhiên tỉ lệ này có thể được tính chỉ với một dòng lệnh bằng hàm clarkevans():
(clarkevans(myspecies_ppp, correction ="non"))## [1] 0.9838335Clark-Evans Test
Tuy nhiên làm sao để biết được giá trị 0,98 này thật sự lớn hơn 1 nhỏ hơn 1 hay bằng 1. Chúng ta cần áp dụng phương pháp kiểm tra thống kê dựa vào chỉ số z. Với giả thuyết H0 là các cá thể phân bố ngẫu nhiên trong khu vực nghiên cứu.
\[z=\frac{\bar{r}_{A} - \bar{r}_{E}}{S\bar{r}}\] Với \[S\bar{r}=\sqrt{\frac{4 - \pi}{4\pi}\frac{1}{n\lambda}}=\frac{0.2613616}{\sqrt{n\lambda}}\]
Nguồn: https://docs.qgis.org
Giả sử chọn \(\alpha = 5\%\) thì nếu chỉ số z nằm trong khoảng * -1,96 < z < 1,96 : ta chấp nhận giả thuyết H0, các cá thể phân bố ngẫu nhiên.
* z <= -1,96 : Các cá thể phân bố cụm
* z >= 1,96 : Các cá thể phân bố đều
Ta tính z như sau:
SE <- sqrt((4 - pi)/(4 * pi)/(npts^2/areaW))
(z <- (Dobs - Dpois)/SE)## [1] -0.3092755Với z = -0,3, ta chấp nhận giả thuyết H0: Các cá thể của loài nghiên cứu phân bố ngẫu nhiên trong không gian.
Hàm clarkevans.test
Tất cả các bước trên có thể được thực hiện đơn giản với lệnh clarkevans.test()
clarkevans.test(myspecies_ppp)##
## Clark-Evans test
## No edge correction
## Z-test
##
## data: myspecies_ppp
## R = 0.98383, p-value = 0.7571
## alternative hypothesis: two-sidedR = 0.98383 chính là tỉ lệ Clark-Evans. Hàm còn cung cấp giá trị p-value = 0.7571 để chúng ta có thể đưa ra nhận định. Trong trường hợp này là chấp nhận giả thuyết H0, tương tự như chúng ta đã tính toán thủ công.
Lời kết
Như các bạn đã thấy chúng ta chỉ cần thực hiện 1 dòng lệnh clarkevans.test() là đã có thể giải được bài toán về sự phân bố của loài trong không gian. Phương pháp này thường được áp dụng ngược lại để tìm khoảng cách tối ưu giữa 2 cá thể trong công tác trồng rừng, để các cá thể có thể phát triển tốt.
Chúc các bạn thực hành vui vẻ.