Sẽ trở lại đây...

S

au một thời gian vắng bóng. Giờ đây tôi sẽ quay lại với cái đam mê của mình, tức là “tập tành, bắt chước...” những mẹo hay có liên quan đến R. Thú thật. Trong thời gian vừa qua, do hơi bận rộn cho việc nghiệm thu nhiệm vụ cũng như công việc cuối năm. Thực ra, hiện giờ không phải mình rảnh cho lắm, tuy nhiên, vẫn cố gắng tìm hiểu, tập tành, bắt chước... rồi cùng trao đổi, chia sẻ với những ai quan tâm để cùng nhau học hỏi, tiến bộ hơn trong vấn đề ứng dụng R trong học tập, công việc, nghiên cứu.

Thú thực một điều. Khi ứng dụng R trong việc xử lý, phân tích cũng như vẽ các biểu đồ trong các báo cáo nghiên cứu. Nhiều người vẫn không tin đó là những hình vẽ, biểu đồ... mà chính người “chấp bút” tạo nên, đặc biệt những ai chưa nghe đến “R”. Cứ nghĩ là “cóp nhặt” đâu đó, mà còn không trích dẫn nguồn. Tuy nhiên, không hẳn mình đã có đủ khả năng hiểu hết về ý nghĩa của các hình, biểu đồ mà mình tạo ra, nên việc giải thích, diễn giải ý nghĩa của các biểu đồ còn nhiều hạn chế. Thực ra, đó càng là “động lực” cho mình tìm hiểu và muốn hoàn thiện hơn.

Trong thời gian vừa rồi, khi đi xin ý kiến các bác chuyên gia có liên quan đến chuyên môn cũng bị hỏi về vấn đề “hình vẽ, biểu đồ” lấy nguồn ở đâu mà không “trích nguồn”. Có lẽ, việc giải thích, diễn giải ý nghĩa của các biểu đồ thực sự sẽ là thách thức đối với cá nhân trong thời gian tới, đặc biệt trong các hội đồng nghiệm thu các nhiệm vụ (toàn những cây đa cây đề), trước nay chúng ta hay quen với các loại bản đồ đơn giản (excel, spss...). 

Read More

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

Phân tích hậu định với R

Bài viết này xuất phát từ nhận định cũng như câu hỏi của Thầy Emu Trần (Trần Bình Đà) khi mình có đăng mấy hình ảnh về biểu đồ hộp thể hiện một số chỉ tiêu sinh trưởng của cây Keo lá liềm được thử nghiệm ở các CTTN khác nhau. Khi đó mình cũng mới dừng lại ở việc phân tích thống kê đơn giản, tức là dừng lại ở việc kiểm sự khác biệt có ý nghĩa thống kê hay không giữa các CTTN về một số chỉ tiêu nghiên cứu, mà chưa đi sâu kiểm tra xem giữa các các công thức nào (so sách cặp đôi) có sự sai khác có ý nghĩa thống kê.

Trong phạm vi bài viết này, mình giới thiệu cách phân tích hậu định trong phân tích phương sai (soạn theo bài giảng của GS Nguyễn Văn Tuấn), vừa là để trả lời câu hỏi rất hữu ích và ý nghĩa của Thầy Emu Trần. Nhân đây cũng xin được gửi lời cảm ơn sâu sắc đến GS Nguyễn Văn Tuấn đã có những bài giảng, tài liệu... giới thiệu về R, nghiên cứu khoa học rất hữu ích và là người truyền lửa về nghiên cứu khoa học cho thế hệ sau, người luôn trăn trở cho nền giáo dục, khoa học nước nhà. Và được xin gửi lời cảm ơn về tất cả, không hẳn chỉ là câu hỏi rất chi là hữu ích của Thầy Emu Trần (Emu Trần Theo thầy thấy, kết quả trên của em không có sự khác biệt có ý nghĩa cho cả các công thức và ĐC. Kết quả TB của ĐC có cao hơn các CT, nhưng giao động của số liệu thì không cho thấy có sự khác biệt. Em thử phân tích hậu định bằng Kruskal xem sao) (trích từ facebook Emu Trần).

Trước tiên mình vẽ hình mà bữa trước mình có đăng và từ đó gợi mở ra nhiều điều.

Codes để vẽ hình như sau:

> 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)+ theme_bw()+ theme_classic()+ xlab("CTTN")+ ylab("tree height, m")+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()

> Age1.2tp2

Kết quả kiểm tra thống kê:

Biểu đồ hộp trên cho ta cái nhìn tổng quan về chỉ tiêu sinh trưởng đường kính gốc của cây Keo lá liềm về một số thông tin như: Bách phân vị 75%, 50%, 25%, outlier (chấm đỏ)...

Kết quả cho CT3 như sau:

> summary(TP_CT3_1.2)          Local          Age        CTTN     stump_diameter  Cam_Duong  :  0   Min.   :1.2   CT 1:  0   Min.   :0.320   Le_Thuy    :  0   1st Qu.:1.2   CT 2:  0   1st Qu.:2.147   Trieu_Phong:110   Median :1.2   CT 3:110   Median :2.580                     Mean   :1.2   DC  :  0   Mean   :2.675                     3rd Qu.:1.2   DC  :  0   3rd Qu.:3.250                     Max.   :1.2              Max.   :4.650    tree_height     canopy_diameter   main_trunk     bough_50_cm    Min.   :0.3500   Min.   :0.500   Min.   :1.000   Min.   :0.000   1st Qu.:0.7000   1st Qu.:0.900   1st Qu.:2.000   1st Qu.:2.000   Median :0.9100   Median :1.140   Median :4.000   Median :2.000   Mean   :0.9543   Mean   :1.162   Mean   :3.555   Mean   :2.318   3rd Qu.:1.2000   3rd Qu.:1.350   3rd Qu.:5.000   3rd Qu.:3.000   Max.   :2.0000   Max.   :3.550   Max.   :8.000   Max.   :5.000     phan_than        song_chet stump_diameter_growth tree_height_growth  Min.   :0.0000   Min.   :1   Min.   :0.270         Min.   :0.2900     1st Qu.:0.0000   1st Qu.:1   1st Qu.:1.792         1st Qu.:0.5800     Median :0.0000   Median :1   Median :2.150         Median :0.7600     Mean   :0.1273   Mean   :1   Mean   :2.229         Mean   :0.7958     3rd Qu.:0.0000   3rd Qu.:1   3rd Qu.:2.710         3rd Qu.:1.0000     Max.   :1.0000   Max.   :1   Max.   :3.870         Max.   :1.6700     canopy_diameter_growth  litter_fall     Min.   :0.4200         Min.   :0.0000   1st Qu.:0.7500         1st Qu.:0.1000   Median :0.9500         Median :0.1200   Mean   :0.9675         Mean   :0.1265   3rd Qu.:1.1300         3rd Qu.:0.1600   Max.   :2.9600         Max.   :0.2300

Có thể diễn giải kết quả cho CT3 như sau: Chỉ tiêu đường kính gốc (bôi vàng), có 75% đối tượng cây Keo lá liềm điều tra có sinh trưởng đường kính gốc bằng 3,25 cm hoặc nhỏ hơn (3rd Qu.:3.250); có 50% đối tượng có sinh trưởng đạt 2,67 cm hoặc nhỏ hơn và có 25% đối tượng có sinh trưởng đường kính gốc đạt 2,147cm hoặc nhỏ hơn (1st Qu.:2.147).

Để  kiểm tra thống kê về chỉ tiêu sinh trưởng đường kính gốc cây Keo lá liềm ở 04 CTTN có khác nhau rõ rệt hay không (có ý nghĩa thống kê với mức độ tin cậy 95%)?  Ta sử dụng phân tích phương sai (ANOVA) để kiểm tra.

Trong R có thể sử dụng lệnh như sau:

> ao=aov(stump_diameter~CTTN) > summary(ao) Kết quả như sau:              Df Sum Sq Mean Sq F value  Pr(>F)   CTTN          3   7.67   2.555   4.413 0.00482 ** Residuals   243 140.70   0.579                   --- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

>

Kết quả cho thấy, trị số pvalue = 0,00482 < 0,05 (Ft = 4,413 > Flt = 2,06), tức là 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?

Để trả lời câu hỏi này, chúng ta có thể sử dụng phân tích hậu định, tức là phân tích hậu định trong phân tích phương sai.

Tuy nhiên, vấn đề đặt ra hiện nay là có nhiều phương pháp phân tích hậu định như: Fisher’s method (LSD), Bonferroni’s method, Duncan’s mutiple range test, Scheffe’, Tukey’s Honest Significant Difference, Dunnett’s test... Vậy phương pháp nào thích hợp nhất?

Trong phạm vi bài viết mình có sử dụng phương pháp Tukey’s Honest Significant Difference trong R để kiểm tra.

Kết quả như sau:

> ao=aov(stump_diameter~CTTN)

> TukeyHSD(ao)   Tukey multiple comparisons of means     95% family-wise confidence level Fit: aov(formula = stump_diameter ~ CTTN) $CTTN                 diff         lwr        upr     p adj CT 2-CT 1 -0.2674558 -0.63761365 0.10270205 0.2441079 CT 3-CT 1  0.0912766 -0.25175065 0.43430384 0.9014797 DC-CT 1   -0.3787234 -0.91386696 0.15642015 0.2614210 CT 3-CT 2  0.3587324  0.05906615 0.65839864 0.0116606 DC-CT 2   -0.1112676 -0.61970717 0.39717196 0.9420240 DC-CT 3   -0.4700000 -0.95904156 0.01904156 0.0646052

>

Kết quả cho thấy, chỉ có CT3-CT2 là có sự khác biệt có ý nghĩa thống kê với độ tin cậy 95% về chỉ tiêu sinh trưởng đường kính gốc. Các so sánh còn lại chưa có sự khác biệt có ý nghĩa thống kê (giá trị pvalue > 0,05). Theo đó, đường kính gốc ở CT3 cao hơn có ý nghĩa 0,36cm, dao động từ 0,06cm đến 0,66cm so với CT2. 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(ao), ordeder=T)

Kết quả hình trên có thể giải thích như sau: Trên biểu đồ cho thấy những cặp đôi CT nào không cắt đường giá trị 0,0 (nét đứt) thì cặp đôi so sánh đó có ý nghĩa thống kê, tức là lệch hẳn về phía (âm, dương) của đường giá trị. Ta thấy, chỉ có CT3-CT2 không cắt ngang đường nét đứt, tức là có sự khác biệt có ý nghĩa thống kê. Còn lại là không có sự khác biệt có ý nghĩa thống kê.

Read More