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

Ở bài 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, nghĩa 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?

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) đều 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 TP HCM, tr.142-143.