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ê.