Giới thiệu trang blog về “ứng dụng R trong lâm nghiệp”

Trong cái note này, mình giới thiệu đến quý bạn đọc, những ai quan tâm đến R, đặc biệt ứng dụng R trong nghiên cứu lâm nghiệp về một trang fanpage: https://www.facebook.com/ungdungRtronglamnghiep/ và một trang blog cá nhân: http://ungdungrtronglamnghiep.blogspot.com/ có liên quan đến R và ứng dụng R trong nghiên cứu lâm nghiệp mà cá nhân mình tự tìm hiểu, tập tành và ứng dụng đến những vấn đề có liên quan đến công việc (thực hiện các nhiệm vụ KHCN các cấp), tức là có sử dụng thống kê và R trong xử lý số liệu có liên quan đến lĩnh vực lâm - nông nghiệp.

Thực sự mà nói, bản thân cũng đang trong quá trình tự tìm hiểu, tập tành, bắt chước... cũng tạm gọi là biết chút xíu, qua đó ứng dụng vào những thực tế các vấn đề có liên quan. Vì vậy, trên quan điểm cùng nhau trao đổi, học hỏi lẫn nhau về những thông tin, mẹo hay, codes hay... về R, chứ không dám “khoe” với ai về bất kỳ một vấn đề gì, càng không dám “lên lớp” bất kỳ một ai, dù là nhỏ nhất cũng như không có mục đích nào khác.

Hy vọng, qua trang fanpage và trang blog cá nhân có liên quan đến R và ứng dụng R trong nghiên cứu lâm nghiệp, mọi người sẽ cùng trao đổi, học tập lẫn nhau và cùng hướng tới “ứng dụng sâu rộng R trong học tập, công việc, đặc biệt trong nghiên cứu lâm - nông nghiệp”.

Read More

Hoán chuyển dữ liệu

Kiểm định thống kê

Khi kiểm định ý nghĩa thống kê (test of signifcance), phần lớn các phép tính dựa vào giả định biến số phải là một biến số phân phối chuẩn (normal distributin). Do đó, một trong những việc quan trọng khi xem xét dữ kiện là phải kiểm định giả thuyết phân phối chuẩn của một biến số [1].

Ở Phân tích hậu định với R, mình phân tích hậu định trong phân tích phương sai nhằm đi sâu kiểm tra xem giữa các công thức thí nghiệm (CTTN), công thức nào (so sánh cặp đôi) có sự sai khác có ý nghĩa thống kê. Trong khi đó, thông thường các phân tích (không phải đa số) mới dừng ở việc phân tích thống kê đơn giản; nói cách khác, mới dừng lại ở việc kiểm tra sự khác nhau có ý nghĩa thống kê hay không giữa các CTTN về một chỉ tiêu nghiên cứu nào đó.

Trong phạm vi bài viết, mình thực hiện phân tích hậu định về chỉ tiêu nghiên cứu đường kính gốc (stump diameter) của cây Keo lá liềm trồng trên đất cát vùng ven biển Bắc Trung bộ. Tuy nhiên, trước khi phân tích chúng ta cần kiểm định giả thiết phân phối của biến stump diameter có phải tuân theo phân phối chuẩn hay không. Đến đây có 2 trường hợp xảy ra: (1) Biến stump diameter tuân theo phân phối chuẩn, việc phân tích hậu định diễn ra bình thường; (2) Biến stump diameter không tuân theo luật phân phối chuẩn, cần có bước trung gian, tức là hoán chuyển dữ liệu để khắc phục những giá trị bất thường như skewness (độ lệch)... giúp đưa biến stump diameter về gần với phân phối chuẩn và thực hiện các bước tiếp theo.

Trước tiên, chúng ta kiểm định giả thiết phân phối chuẩn của biến stump diameter:

Bạn có thể có cái nhìn tổng quan phân bố biến stump diameter qua hình sau với lệnh hist (x)

> hist(stump_diameter)

Qua hình trên có thể nhận thấy, biến stump diameter không tuân theo luật phân phối chuẩn (lệch về các giá trị nhỏ).

Kiểm tra các giá trị skewness và Kurtosis

> describeBy(stump_diameter)

  vars   n mean   sd median trimmed  mad min  max range skew kurtosis   se

1    1 391 3.71 1.44    3.3    3.56 1.19 1.4 11.6  10.2 1.15     2.05 0.07

Một biến có phân phối chuẩn khi giá trị skewness (độ lệch) tiến gần tới giá trị 0 và kurtosis tiến gần tới giá trị
3. Kết quả (bôi vàng) cho thấy, có thể nhận định biến stump diameter không tuân theo luật phân phối chuẩn.


Theo CTTN:

> describeBy(stump_diameter, group=CTTN)
group: CT 1
  vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
1    1 64 3.22 0.95      3    3.17 0.96 1.4   6   4.6 0.51     0.01 0.12
------------------------------------------------------------ 
group: CT 2
  vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
1    1 69  3.3 1.31    3.2    3.14 1.04 1.4 7.2   5.8  1.2     1.21 0.16
------------------------------------------------------------ 
group: CT 3
  vars   n mean   sd median trimmed  mad min  max range skew kurtosis  se
1    1 228 3.81 1.47    3.5    3.67 1.19 1.4 11.6  10.2 1.27     2.88 0.1
------------------------------------------------------------ 
group: DC
  vars  n mean   sd median trimmed  mad min max range skew kurtosis   se
1    1 30 4.88 1.57   4.75    4.93 1.93   2 7.3   5.3 -0.1    -1.26 0.29

 Kết quả theo CTTN cũng cho nhận định

rằng, biến stump diameter không tuân theo luật phân phối chuẩn.

Tuy nhiên, để kiểm định nghiêm chỉnh
(bằng chứng khoa học - dựa trên kết quả cụ thể) xem biến stump diameter có tuân
theo luật phân phối chuẩn hay không? chúng ta cần phải sử dụng kiểm định thống
kê. Trong R, chúng ta sử dụng hàm shapiro.test (x), cụ thể như
sau:



 


  
> shapiro.test(stump_diameter)
 
        Shapiro-Wilk normality test
 
data:  stump_diameter
W = 0.96199, p-value = 4.02e-06
 

 


Kết quả cho thấy, trị số p (p-value)
= 4,02e-06 << 0,05, chúng ta có thể kết luận rằng biến số stump diameter
không đáp ứng luật phân phối chuẩn.

Như vậy, biến stump diameter không
đáp ứng luật phân phối chuẩn. Vì vậy, trước khi thực hiện các phân tích thống
kê tiếp theo, chúng ta cần phải hoán chuyển dữ liệu nhằm khắc phục
các giá trị bất thường của skewness, giúp đưa phân phối của biến stump diameter
về gần với phân phối chuẩn. Hiện có nhiều cách để hoán chuyển dữ liệu khi biến số 

không tuân theo luật phân phối chuẩn
như: sử dụng hàm logarit để hoán chuyển dữ liệu, phương pháp hoán chuyển dữ liệu
dựa vào hàm lũy thừa (phương pháp Box-Cox)... Tuy nhiên, trong phạm vi bài viết,
mình giới thiệu phương pháp hoán chuyển dữ liệu dựa vào hàm logarit, còn phương
pháp Box-Cox ở dịp khác mình sẽ giới thiệu sau.

Trước khi đi vào phân tích, để có cái

nhìn tổng quát về biến stump diameter theo các CTTN, mình vẽ biểu đồ hộp
(boxplot) dưới đây:

>
Age1.2tp2=ggplot(data=Age1.2tp, aes(Age1.2tp$CTTN, y=stump_diameter))+
geom_boxplot(aes(fill=CTTN), outlier.colour="red", outlier.size=2.7,
notch=T)+ theme_bw()+ theme_classic()+ xlab("CTTN")+ ylab("stump
diameter, cm")+ggtitle("A. crassicarpa 14 months of age in Trieu
Phong")+ geom_rangeframe()+ theme_tufte()+
scale_y_continuous(breaks=extended_range_breaks()(Age1.2$stump_diameter))+
theme(legend.position="top")+ coord_flip()+ geom_jitter(alpha=.2,
size=2, shape=16, color="blue")



Để tiện cho việc so sánh, chúng ta
cùng phân tích cả 2 trường hợp. Trường hợp 1, kiểm định thống kê thông thường
khi chưa hoán chuyển dữ liệu. Trường hợp 2, hoán chuyển dữ liệu trước khi thực
hiện các phân tích thống kê tiếp theo.

Hoán chuyển dữ liệu bằng hàm logarit
trong R như sau:

> log.stump_diameter=log(stump_diameter)
> log.stump_diameter
 
  [1] 1.3350011 0.8329091 1.1314021 0.5877867 0.9162907 1.1939225 1.1939225
  [8] 0.7419373 1.1631508 1.2237754 1.0647107 1.6292405 0.6931472 1.0296194
 [15] 1.1939225 0.6931472 0.9162907 0.6931472 0.8329091 1.0986123 1.0986123
 [22] 0.9932518 0.9932518 1.0647107 0.7884574 0.7419373 0.3364722 0.5877867
 [29] 1.0647107 0.5877867 0.5306283 1.0986123 1.0647107 0.8754687 1.3350011
 [36] 1.0986123 1.4586150 1.3609766 1.0986123 1.2527630 1.5686159 0.9162907
 [43] 1.0986123 0.9932518 0.8754687 1.0647107 1.1631508 1.3083328 0.5306283
 [50] 0.8329091 1.0986123 1.1939225 1.1631508 0.9932518 1.3083328 0.6418539
 [57] 0.8329091 0.5877867 0.7884574 1.0986123 0.6418539 1.5040774 0.9555114
 [64] 0.8329091 1.6486586 0.7884574 0.8329091 0.8754687 1.0647107 0.6418539
 [71] 1.0647107 0.8329091 0.9162907 1.3350011 0.7884574 0.9162907 0.8329091
 [78] 0.6931472 1.0986123 1.0986123 1.0647107 1.0296194 1.1314021 1.2809338
 [85] 1.3083328 1.3862944 1.4586150 1.7227666 1.3862944 1.3609766 1.3862944
 [92] 0.9932518 0.9162907 0.9932518 0.6418539 0.9932518 1.3083328 0.8329091
 [99] 1.0986123 1.3083328 1.0647107 1.2527630 1.0296194 0.6418539 1.2527630
[106] 1.2809338 0.8754687 1.0296194 0.8329091 1.4586150 0.7419373 1.1631508
[113] 0.8754687 0.7884574 1.1939225 1.0986123 1.1939225 1.1631508 1.3350011
[120] 1.1631508 1.5475625 1.5686159 0.7884574 1.3350011 1.0647107 1.5040774
[127] 1.3350011 1.5040774 1.0647107 1.7917595 1.4109870 1.3609766 1.6292405
[134] 1.6292405 1.5040774 0.9555114 1.0986123 1.3609766 1.3350011 1.0647107
[141] 0.9555114 1.0647107 0.4700036 1.0647107 0.4700036 0.9555114 1.0647107
[148] 0.9162907 0.3364722 1.1939225 1.1939225 1.6094379 1.3083328 1.0986123
[155] 0.6931472 1.3083328 1.1631508 0.9162907 0.7884574 0.9932518 1.1939225
[162] 0.6931472 1.1939225 1.5475625 1.1631508 1.2527630 0.7884574 0.8754687
[169] 1.2527630 0.6931472 1.1631508 0.9555114 1.5892352 0.9932518 1.9021075
[176] 1.6677068 1.0296194 1.3609766 1.4109870 1.0986123 1.5260563 1.4586150
[183] 1.7749524 1.3350011 1.1631508 0.9932518 0.6931472 0.6418539 1.1631508
[190] 1.1939225 1.0647107 1.1939225 0.7884574 1.0986123 1.6677068 1.3350011
[197] 1.9169226 1.5040774 0.7419373 0.7884574 0.7884574 1.9740810 1.9459101
[204] 1.4109870 0.9932518 1.3350011 1.0647107 0.5877867 1.1939225 1.2809338
[211] 0.4054651 0.6418539 0.9162907 1.3609766 0.9932518 1.3350011 0.6418539
[218] 0.4700036 1.3083328 0.9932518 1.0986123 0.7884574 1.1631508 0.3364722
[225] 0.9162907 0.8754687 1.2809338 1.2527630 1.0986123 1.0986123 1.3083328
[232] 1.3862944 1.1939225 1.1631508 1.2527630 0.9162907 0.9932518 1.1631508
[239] 0.9162907 1.3083328 1.6292405 1.0986123 0.9932518 0.9932518 1.3083328
[246] 1.4586150 1.8718022 1.6677068 1.0647107 1.3609766 1.4586150 1.1939225
[253] 1.5686159 0.6418539 1.3083328 1.3862944 1.8082888 1.7227666 1.6863990
[260] 1.7404662 1.2527630 1.4586150 1.1631508 1.0647107 1.2809338 1.1631508
[267] 1.4816045 1.3083328 1.1631508 1.3083328 1.6677068 1.2527630 1.3350011
[274] 0.9162907 1.6486586 0.8329091 2.4510051 2.0794415 2.0149030 2.0541237
[281] 2.0541237 1.4350845 1.7227666 1.2527630 1.6292405 1.1939225 1.6486586
[288] 1.2527630 1.7749524 1.9021075 1.1939225 1.2809338 1.5260563 1.9169226
[295] 1.7578579 1.8562980 1.9169226 1.5040774 1.3609766 0.9555114 1.3609766
[302] 1.1314021 0.9162907 1.1631508 1.3350011 1.2809338 1.5040774 1.3862944
[309] 1.5686159 1.6094379 1.3350011 1.3609766 1.3350011 1.9021075 1.0647107
[316] 1.6486586 1.8718022 1.6677068 1.4586150 1.4109870 1.4350845 1.2237754
[323] 1.4586150 1.1939225 1.5686159 1.5686159 1.7749524 1.6677068 1.6094379
[330] 1.6677068 1.7749524 1.9878743 1.9021075 1.9459101 1.4109870 1.3862944
[337] 1.1939225 1.7404662 1.5892352 1.1314021 0.9932518 1.5040774 1.6677068
[344] 1.7227666 1.8082888 1.5892352 1.5040774 1.5892352 1.5686159 1.1314021
[351] 0.7419373 1.3609766 1.7749524 1.5260563 1.5260563 1.3862944 1.0986123
[358] 1.0296194 1.8245493 1.7917595 1.5260563 1.4586150 1.1939225 1.3083328
[365] 0.6931472 1.8405496 1.9021075 1.2237754 1.2527630 1.7047481 1.9315214
[372] 1.5892352 1.4586150 1.9459101 1.5260563 1.7404662 0.8329091 1.4586150
[379] 0.8754687 1.7404662 1.4586150 1.9878743 1.9459101 1.7227666 1.3350011
[386] 1.6863990 1.8405496 1.8245493 1.0986123 1.3083328 1.9459101

Sau khi hoán chuyển, chúng ta tiến

hành các phân tích thống kê thông thường

> hist(log.stump_diameter)





> summary(log.stump_diameter)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3365  0.9933  1.1940  1.2420  1.5040  2.4510 
 
> describeBy(log.stump_diameter)
  vars   n mean   sd median trimmed  mad  min  max range skew kurtosis   se
1    1 391 1.24 0.37   1.19    1.24 0.39 0.34 2.45  2.11 0.13     -0.3 0.02
 
> describeBy(log.stump_diameter, group=CTTN)
group: CT 1
  vars  n mean  sd median trimmed  mad  min  max range  skew kurtosis   se
1    1 64 1.13 0.3    1.1    1.14 0.31 0.34 1.79  1.46 -0.28    -0.12 0.04
------------------------------------------------------------ 
group: CT 2
  vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
1    1 69 1.13 0.37   1.16    1.11 0.37 0.34 1.97  1.64 0.27    -0.26 0.04
------------------------------------------------------------ 
group: CT 3
  vars   n mean   sd median trimmed  mad  min  max range skew kurtosis   se
1    1 228 1.27 0.36   1.25    1.27 0.38 0.34 2.45  2.11 0.16    -0.27 0.02
------------------------------------------------------------ 
group: DC
  vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
1    1 30 1.53 0.36   1.56    1.56 0.41 0.69 1.99  1.29 -0.6    -0.62 0.07



Phân tích tích phương sai (ANOVA) cho biến stump diameter sau khi đã hoán chuyển dữ liệu.

 > ao1=aov(log.stump_diameter~CTTN)
> summary(ao1)
             Df Sum Sq Mean Sq F value   Pr(>F)    
CTTN          3   4.44  1.4797   11.73 2.28e-07 ***
Residuals   387  48.83  0.1262                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Trường hợp chưa hoán chuyển dữ liệu:


> ao2=aov(stump_diameter~CTTN)
> summary(ao2)
             Df Sum Sq Mean Sq F value   Pr(>F)    
CTTN          3   70.1  23.379   12.33 1.02e-07 ***
Residuals   387  733.9   1.896                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Kết quả ở cả 2 trường hợp trên cho thấy, trị số p-value =
2,28e-07 (đã hoán chuyển dữ liệu) và p-value= 1,02e-07 đều << 0,05. Có sự khác biệt có ý nghĩa 

thống kê với mức độ tin cậy 95% giữa các CTTN về sinh trưởng đường kính gốc của cây Keo
lá liềm. Tuy nhiên, kết quả không cho ta biết sự khác biệt giữa công thức nào
với công thức nào? Ta có 06 nhóm: ĐC-CT3, ĐC-CT2, ĐC-CT1, CT3-CT1, CT3-CT2 và
CT2-CT1. Vậy câu hỏi đặt ra là sự khác biệt có ý nghĩa thống kê về giữa
nhóm công thức nào?

Sử dụng phương pháp Tukey’s Honest Significant Difference trong R để kiểm tra


> TukeyHSD(ao1)

  Tukey multiple comparisons of means
    95% family-wise confidence level
 Fit: aov(formula = log.stump_diameter ~ CTTN)
 
$CTTN
                   diff         lwr       upr     p adj
CT 2-CT 1 -0.0006330752 -0.15969504 0.1584289 0.9999996
CT 3-CT 1  0.1453260317  0.01567109 0.2749810 0.0209890
DC-CT 1    0.4013840991  0.19858391 0.6041843 0.0000031
CT 3-CT 2  0.1459591069  0.02002560 0.2718926 0.0156276
DC-CT 2    0.4020171743  0.20157575 0.6024586 0.0000022
DC-CT 3    0.2560580675  0.07805116 0.4340650 0.0013451
 
Kết quả cho thấy, CT2-CT1 là chưa có sự khác biệt có ý nghĩa thống kê
 với độ tin cậy 95% về chỉ tiêu stump diameter (p-value = 0,999 > 0,05). 
Các công thức cặp đôi còn lại (Ct3-CT1, DC-CT1, CT3-CT2, DC-CT2, DC-CT3) 
đề có sự khác biệt có ý nghĩa thống kê (p-value < 0,05). Ngoài ra, chúng ta 
có thể so sánh sự khác biệt đó bằng biểu đồ với lệnh sau:

> plot(TukeyHSD(ao1), ordeder=T)

Hy vọng sau ví dụ này, chúng ta biết

cách hoán chuyển dữ liệu trước khi tiến hành các phân tích thống kê tiếp theo
khi biến phân tích chưa đáp ứng luật phân phối chuẩn. Ở bài tiếp theo, mình xin
được giới thiệu phương pháp hoán chuyển dữ liệu dựa vào hàm lũy thừa (phương
pháp Box-Cox) mà mình có đề cập bên trên (chỉ là học theo, bắt chước thôi).

=================================

[1] Nguyễn Văn Tuấn (2014). Phân tích dữ liệu với R. Nxb Tổng hợp TPHCM, tr.142-143.

Read More