Chương 1. • Nhập dữ liệu. Trước hết, chúng ta cần phải nhập dữ liệu vào R. Bước thứ nhất là báo cho R biết rằng chúng ta có ba nhóm sinh viên học 3 giáo sư A, B và C với nhóm tương ứng (1, 2 và 3), nhóm 1 gồm 4 người, nhóm 2 có 4 người, và nhóm 3 có 4 người: group<-c(1,1,1,1,2,2,2,2,3,3,3,3) Để phân tích phương sai, chúng ta phải định nghĩa biến group là một nhân tố - factor.factor(group) 15 Bước kế tiếp, chúng ta nhập số điểm của 3 nhóm sinh viên đạt được >diem<-c(79,86,94,89,71,77,81,82,82,68,70,76) Đưa hai biến group và điểm vào một dataframe và gọi là data: > data <- data.frame (group, diem) > attach(data) Sau khi đã có dữ liệu sẵn sàng, chúng ta dùng hàm lm( ) để phân tích phương sai như sau: > analysis <- lm (diem ~ group) Trong hàm trên chúng ta cho R biết biến diem là một hàm số của group.
Gọi kết quả phân tích là analysis. • Kết quả phân tích phương sai. Tiếp theo, chúng ta dùng lệnh anova để biết kết quả phân tích: Analysis of Variance Table Response: diem Df Sum Sq Mean Sq F value Pr(>F) group 2 354.1 ‘ ’ 1 Trong kết quả trên, có ba cột: Df (degrees of freedom) là bậc tự do; Sum Sq là tổng bình phương (sum of squares), Mean Sq là trung bình bình phương (mean square); F value là giá trị kiểm định F. Để biết thêm chi tiết kết quả phân tích, chúng ta dùng lệnh summary như sau: summary(analysis) Call: lm(formula = diem ~ group) Residuals: Min 1Q Median 3Q Max -8.0 16 Coefficients: Estimate Std.
Error t value Pr(>|t|) (Intercept) 87.1 ‘ ’ 1 Residual standard error: 5.981 on 9 degrees of freedom Multiple R-squared: 0.957 on 2 and 9 DF, p-value: 0.03537 Tính toán trên R cho ta kết quả tương tự trong Ví dụ 1. Chúng tôi sẽ phân tích sâu hơn các kết quả của R in ra trong ví dụ tiếp theo. Cho số liệu trong Bảng 2.1 dưới đây so sánh độ galactose trong 3 nhóm bệnh nhân: nhóm 1 gồm 9 bệnh nhân với bệnh Crohn; nhóm 2 gồm 11 bệnh nhân với bệnh viêm ruột kết (colitis); và nhóm 3 gồm 20 đối tượng không có bệnh (gọi là nhóm đối chứng). Câu hỏi đặt ra là độ galactose giữa 3 nhóm bệnh nhân có khác nhau hay không? Gọi giá trị trung bình của ba nhóm là µ1 , µ2 , µ3 và nói theo ngôn ngữ của kiểm định giả thiết thì giả thiết đảo là: H0 : µ1 = µ2 = µ3.
Giả thiết chính là HA : có một khác biệt giữa 3 H j (j = 1, 2, 3) Bảng 2.1 Nhóm 1: Nhóm 2: Nhóm 3: Bệnh CRONH Bệnh viêm ruột kết Đối chứng 1343 1264 1809; 2850 1393 1314 1926; 2964 1420 1399 2283; 2973 1641 1605 2384; 3171 1897 2385 2447; 3257 2160 2511 2479; 3271 2169 2514 2495; 3288 2279 2767 2525; 3358 2890 2827 2541; 3643 2895 2769; 3657 3011 n=9 n=11 n=20 Trung bình: 1910 Trung bình: 2226 Trung bình: 2804 Độ lệch chuẩn: 516 Độ lệch chuẩn: 727 Độ lệch chuẩn: 527 17 Như đã biết, ta có thể thực hiện so sánh từng cặp phương sai bằng phương pháp kiểm định t trong chương trình đại học giữa nhóm 1 và 2, nhóm 2 và 3, và nhóm 1 và 3. Tuy nhiên, cách làm này không hợp lí, vì có ba phương sai khác nhau. Cách thích hợp nhất để so sánh này là phân tích phương sai. Phân tích phương sai có thể ứng dụng để so sánh nhiều nhóm cùng một lúc (simultaneous comparisons).1 Mô hình phân tích phương sai Để minh họa cho phương pháp phân tích phương sai, chúng ta phải dùng kí hiệu.
Gọi độ galactose của bệnh nhân i thuộc nhóm j ( j = 1, 2, 3) là xi j. Mô hình phân tích phương sai phát biểu rằng: xi j = µ + αi + εi j (2.1) Tức là, giá trị galactose của bất cứ bệnh nhân nào bằng giá trị trung bình của toàn quần thể (µ) cộng/trừ cho ảnh hưởng của nhóm j được đo bằng hệ số ảnh hưởng αi , và sai số εi j. Một giả định khác là εi j phải tuân theo luật phân phối chuẩn với trung bình 0 và phương sai σ2. Hai thông số cần ước tính là µ và αi.
Cũng như phân tích hồi qui tuyến tính, hai thông số này được ước tính bằng phương pháp bình phương nhỏ nhất; tức là tìm ước số µ̂ và α̂i , sao cho 2 ∑ xi j − µ̂ − α̂ i nhỏ nhất. Quay lại với số liệu nghiên cứu trên, chúng ta có những tóm tắt thống kê như sau: Bảng 2.2 Nhóm Số đối tượng (n j ) Trung bình Phương sai 1 - Crohn n1 = 20 x1 = 1910 s21 = 265944 2 - Viêm ruột kết n2 = 11 x2 = 2226 s22 = 473387 3- Đối chứng n3 = 20 x3 = 2804 s23 = 277500 Toàn bộ mẫu n = 20 x = 2444 Dễ thấy xi j = x + (x j − x) + (xi j − x j ) (2.2) trong đó, x là số trung bình của toàn mẫu, và x j là số trung bình của nhóm j. Nói cách khác, phần (x − x j ) phản ánh độ khác biệt (hay cũng có thể gọi là hiệu số) giữa trung bình từng nhóm và trung bình toàn mẫu, và phần (x − x j ) phản ánh hiệu số giữa một galactose của một đối tượng và số trung bình của từng nhóm. Từ đó, chúng ta có các nguồn dao động như sau: • Tổng bình phương cho toàn bộ mẫu là: SST = ∑ ∑(xi j − x)2 = (1343 − 2444)2 + (1393 − 2444)2 + · · · + (3657 − 2444)2 i j = 12133923 18 • Tổng bình phương phản ánh độ khác nhau giữa các nhóm: SSE = ∑ ∑(xi − x)2 = ∑ n j (xi − x)2 i j j = 9(1910 − 2444)2 + 11(2226 − 2444)2 + 20(2804 − 2444)2 = 5681168 • Tổng bình phương phản ánh độ dao động trong mỗi nhóm: SSF = ∑ ∑(xi j − x j )2 = ∑(n j − 1)s2j i j j = (9 − 1)(265944) + (11 − 1)(473387) + (20 − 1)(277500) = 12133922 SSF được tính từ mỗi bệnh nhân cho 3 nhóm, cho nên trung bình bình phương cho từng nhóm (mean square - MSE) là: MSE = SSF/(N − k) = 12133922/(40 − 3) = 327944 và trung bình bình phương giữa các nhóm là: MSF = SSE/(k − 1) = 5681168/(3 − 1) = 2841810 Trong đó N là tổng số bệnh nhân (N = 40) của ba nhóm, và k = 3 là số nhóm bệnh nhân.
Nếu có sự khác biệt giữa các nhóm, thì chúng ta kì vọng rằng MSF sẽ lớn hơn MSE. Thành ra, để kiểm tra giả thiết, chúng ta có thể dựa vào kiểm định F: F = MSF/MSE = 8.3) Với bậc tự do k − 1 và N − k. Các số liệu tính toán trên đây có thể trình bày trong một bảng phân tích phương sai (ANOVA table) như sau: Bảng 2.3 Trung bình Nguồn biến thiên Bậc tự do Tổng bình phương Kiểm định F bình phương Khác biệt giữa các nhóm 2 5681168 2841810 8,6655 Khác biệt trong 37 12133923 327944 từng nhóm Tổng số 39 12133923 2.2 Phân tích phương sai đơn giản với R Ta thấy, các tính toán trên tương đối phức tạp, và tốn khá nhiều thời gian. Tuy nhiên với R, các tính toán đó có thể làm trong vòng 1 giây, sau khi dữ liệu đã được chuẩn bị đúng cách.1 Nhập dữ liệu Trước hết, chúng ta cần phải nhập dữ liệu vào R.
Bước thứ nhất là báo cho R biết rằng chúng ta có ba nhóm bệnh nhân (1, 2 và 3), nhóm 1 gồm 9 người, nhóm 2 có 11 người, và nhóm 3 có 20 người: > group <- c(1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,2,2, 3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3,3) Để phân tích phương sai, chúng ta phải định nghĩa biến group là một nhân tố - factor.factor(group) Bước kế tiếp, chúng ta nhập số liệu galactose cho từng nhóm như định nghĩa trên (gọi object là galactose): > galactose <- c(1343,1393,1420,1641,1897, 2160,2169,2279, 2890,1264,1314,1399,1605,2385,2511,2514, 2767,2827,2895,3011, 1809,2850, 1926,2964, 2283,2973,2384,3171,2447,3257,2479,3271, 2495,3288,2525,3358,2541,3643,2769,3657) Đưa hai biến group và galactose vào một dataframe và gọi là data: > data <- data.frame (group, galactose) > attach(data) Sau khi đã có dữ liệu sẵn sàng, chúng ta dùng hàm lm( ) để phân tích phương sai như sau: > analysis <- lm (galactose ~ group) Trong hàm trên chúng ta cho R biết biến galactose là một hàm số của group. Gọi kết quả phân tích là analysis.2 Kết quả phân tích phương sai Tiếp theo, chúng ta dùng lệnh anova để biết kết quả phân tích: > anova (analysis) Analysis of Variance Table Response: galactose Df Sum Sq Mean Sq F value Pr(>F) group 2 5683620 2841810 8.1 ‘ ’ 1 20 Trong kết quả trên, có ba cột: Df (degrees of freedom) là bậc tự do; Sum Sq là tổng bình phương (sum of squares), Mean Sq là trung bình bình phương (mean square); F value là giá trị kiểm định F xác định trong (2.3); và Pr(> F) là trị số P liên quan đến kiểm định F. Dòng group trong kết quả trên có nghĩa là bình phương giữa các nhóm (between-groups) và residual là bình phương trong mỗi nhóm (within-group). Ở đây, chúng ta có: SSE = 5683620 và MSF = 2841810 và: MSF = 2841810 và MSF = 327944.00082 có nghĩa là tín hiệu cho thấy có sự khác biệt về độ galactose giữa ba nhóm.3 Ước số Để biết thêm chi tiết kết quả phân tích, chúng ta dùng lệnh summary như sau: > summary(analysis) Call: lm(formula = galactose ~ group) Residuals: Min 1Q Median 3Q Max -995.8 Coefficients: Estimate Std.
Error t value Pr(>|t|) (Intercept) 1910.1 ‘ ’ 1 Residual standard error: 572.7 on 37 degrees of freedom Multiple R-squared: 0.666 on 2 and 37 DF, p-value: 0.0008191 Theo kết quả trên đây, intercept chính là µ̂ trong mô hình (2. Nói cách khác µ̂ = 1910 và sai số chuẩn là 190, 9. Để ước tính thông số α̂ j , R đặt α̂1 = 0, và α̂2 = α̂2 − α̂1 = 316, 3, với sai số chuẩn là 257, và kiểm định t = 316.3/257 = 1, 229 với trị số p = 0, 2268, Nói cách khác, so với nhóm 1 (bệnh nhân Crohn), bệnh nhân viêm ruột kết có độ galactose trung bình cao hơn 257, nhưng độ khác biệt này không có ý nghĩa thống kê. 21 Tương tự, α̂2 = α̂2 − α̂1 = 894, 3, với sai số chuẩn là 229.9, kiểm định t = 894, 3/7229, 9 = 3, 89, và trị số p = 0, 00040.
So với bệnh nhân Crohn, nhóm đối chứng có độ galactose cao hơn 894, và mức độ khác biệt này có ý nghĩa thống kê.2 So sánh nhiều nhóm và điều chỉnh trị số p Xét k nhóm, chúng ta có ít nhất là k(k − 1)/2 so sánh. Ví dụ trên có 3 nhóm, cho nên tổng số so sánh khả dĩ là 3 (giữa nhóm 1 và 2, nhóm 1 và 3, và nhóm 2 và 3).