Trong trường hợp so sánh hai nhóm A và B dựa trên một đáp ứng, ta đi so sánh số trung bình `mu_A` và `mu_B` của hai nhóm bằng trắc nghiệm t (t-test) hay trắc nghiệm z (z-test). Khi đó ta kiểm định một cặp giả thuyết, trong đó giả thuyết không là:

    Ho : `mu_A=mu_B`

Tùy thuộc đặc điểm cụ thể của hai nhóm mà ta có nhiều cách thực hiện khác nhau. Thí dụ như khi hai mẫu này có phương sai của tổng thể chưa biết nhưng được giả định là bằng nhau, số phần tử `p_A` và `p_B` của hai nhóm không quá lớn, ta cần tính số thống kê `t` sau:

`t=(bar x_A-bar x_B) / ( sqrt(s_p^2(1/p_A+1/p_B)))`

(12)

trong đó `bar x_A` và `bar x_B` là trung bình của hai nhóm A và B, `s_p` là phương sai gộp của hai nhóm được xác định theo công thức:

`s_p^2=((p_A-1)s_A^2 + (p_B-1)s_B^2) / (p_A+p_B-2)`

(13)

Do `t` tuân theo phân phối t (phân phối Student) với độ tự do `(p_A+p_B-2)` nên ta sẽ so sánh số thống kê `t_o` của mẫu với giá trị tời hạn `t`* để có kết luận.

Trong dữ liệu đa biến, trung bình mỗi nhóm được thể hiện dưới dạng một vectơ `mb(mu)` (vectơ trung bình) gồm `n` thành phần, tương ứng với `n` biến. Như vậy để so sánh hai nhóm A và B ta kiểm định cặp giả thuyết sau:

    Ho : `mb(mu)_A=mb(mu)_B`

    Ho : `mb(mu)_A != mb(mu)_B`

Ta thực hiện tương tự cho trường hợp một biến như sau:

• Số liệu của mỗi phần tử `j` được thể hiện bằng một vectơ `mb(x)_j` gồm `n` thành phần tương ứng với `n` biến,
• Trung bình của `p` phần tử được thể hiện bằng vectơ `bar mb(x)` gồm `n` thành phần tương ứng với n biến, mỗi thành phần `i` là trung bình `bar x_i` của `p` giá trị `x_(ij)`:

  `bar mb(x)=1/p sum_(j=1)^p mb(x)_j`

  (14)

• tương ứng với hai nhóm, ta có hai ma trận hiệp phương sai `mb(S)_A` và `mb(S)_B`, là các ma trận vuông cấp `n`.
• Ma trận hiệp phương sai gộp `mb(S)_p` được xác định theo công thức:

  `mb(S)_p=((p_A-1)mb(S)_A + (p_B-1)mb(S)_B ) / (p_A+p_B-2)`

  (15)

• Để kiểm định cặp giả thuyết trên, ta đi tính số thống kê `T^2` sau:

  `T^2=(p_A\ p_B)/(p_A+p_B) (bar mb(x)_A - bar mb(x)_B)^T mb(S)_p^(-1) (bar mb(x)_A - bar mb(x)_B)`

  (16)

  hay :

  `T^2=(p_A\ p_B)/(p_A+p_B) D^2`

  (17)

  với :

  `D^2=(bar mb(x)_A - bar mb(x)_B)^T mb(S)_p^(-1) (bar mb(x)_A - bar mb(x)_B)`

  (16)

  `D` là khoảng cách suy rộng giữa hai nhóm A và B, còn được gọi là khoảng cách Mahalanobis.
• `T^2` có phân phối `T^2` (hay phân phối Hotelling `T^2`). Vì vậy để kiểm định giả thuyết ta so sánh giá trị `T_o^2` của mẫu với giá trị tới hạn `T^2`* (từ Bảng Hotelling `T^2`). Từ đó rút ra kết luận.

So sánh số trung bình của nhiều nhóm khi chỉ có một biến độc lập (yếu tố) và một biến phụ thuộc (đáp ứng) chính là trường hợp phân tích phương sai một yếu tố thông thường, trong đó ta kiểm định cặp giả thuyết sau:

    Ho : Tất cả trung bình của các nhóm đều bằng nhau

    Ha : Có ít nhất hai số trung bình không bằng nhau

Để kiểm định cặp giả thuyết này, ta lần lượt thực hiện:

• Tính các tổng phương sai chung `SS_T`, tổng phương sai yếu tố `SS_A`, tổng phương sai sai số `SS_E`.
• Tính các trung bình tổng phương sai yếu tố `MS_A`, trung bình phương sai sai số `MS_E` bằng cách lấy tổng phương sai `SS` chia cho độ tự do `df``.
• Tính số thống kê `F` là tỷ số của `MS_A` và `MS_E`.
• Vì `F` có phân phối Fisher nên ta so sánh giá trị `F_o` của mẫu và giá trị tới hạn `F`*.
• Kết quả so sánh cho phép ta kết luận về số trung bình của các nhóm.

Xét một dữ liệu gồm `n` biến của một mẫu có `a` nhóm và mỗi nhóm đều có `p` phần tử (trường hợp cân bằng). Như trên đã trình bày, số liệu của mỗi phần tử `j` thuộc nhóm `k` được thể hiện bằng một vectơ `mb(x)_(kj)`. Như vậy dữ liệu này có thể trình bày trong Bảng 1.

Bảng 1 Dữ liệu của mẫu gồm `a` nhóm, mỗi nhóm `p` phần tử

		Nhóm `1`	Nhóm `2`	. . .	Nhóm `k`	. . .	Nhóm `a`
Phần tử	1	`mb(x)_(11)`	`mb(x)_(21)`	. . .	`mb(x)_(k1)`	. . .	`mb(x)_(a1)`
	2	`mb(x)_(12)`	`mb(x)_(22)`	. . .	`mb(x)_(k2)`	. . .	`mb(x)_(a2)`
	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	j	`mb(x)_(1j)`	`mb(x)_(2j)`	. . .	`mb(x)_(kj)`	. . .	`mb(x)_(aj)`
	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	p	`mb(x)_(1p)`	`mb(x)_(2p)`	. . .	`mb(x)_(kp)`	. . .	`mb(x)_(ap)`
Tổng		`mb(x)_(1•)`	`mb(x)_(2•)`	. . .	`mb(x)_(k•)`	. . .	`mb(x)_(a•)`
Trung bình		`bar mb(x)_1`	`bar mb(x)_2`	. . .	`bar mb(x)_k`	. . .	`bar mb(x)_a`

Ta có một số lưu ý liên quan đến bảng dữ liệu này:

• Không kể các ô tiêu đề, mỗi ô trong bảng này là một vectơ có `n` thành phần tương ứng với `n` biến.
• Ký hiệu "•" để diễn tả tổng, nghĩa là `mb(x)_(k•)` là một vec tơ có `n` thành phần, mỗi thành phần là tổng của `p` phần tử trong nhóm `k`.
• `bar mb(x)_k` là một vec tơ trung bình nhóm, có `n` thành phần, mỗi thành phần là trung bình của `p` phần tử trong nhóm `k`.
• `bar mb(x)` là vec tơ trung bình chung, có `n` thành phần, mỗi thành phần là trung bình của `ap` phần tử của dữ liệu.

Phân tích phương sai cho dữ liệu đa biến về thực chất cũng là kiểm định cặp giả thuyết:

    Ho : Tất cả vectơ trung bình của các nhóm đều bằng nhau

    Ha : Có ít nhất hai vectơ trung bình không bằng nhau

Để kiểm định cặp giả thuyết này, về nguyên tắc ta cũng lần lượt thực hiện các bước như khi phân tích phương sai một yếu tố. Tuy nhiên do tính chất phức tạp của xử lý đa biến nên phương thức thực hiện cũng có một số điểm khác biệt.

Trước hết ta tính hai ma trận sau :

• Tổng phương sai nghiệm thức `SS_A` mở rộng thành ma trận `mb(H)` được xác định bằng công thức:

  `mb(H)=p sum_(k=1)^a (bar mb(x)_k - bar mb(x)) (bar mb(x)_k - bar mb(x))^T`

  (19)

• Tổng phương sai sai số `SS_E` mở rộng thành ma trận `mb(E)` được xác định như sau:

  `mb(E)=sum_(k=1)^a sum_(j=1)^p (mb(x)_(kj) - bar mb(x)_k) (mb(x)_(kj) - bar mb(x)_k)^T`

  (20)

• Trong trường hợp cân bằng, độ tự do của `mb(H)` và `mb(E)` lần lượt là `(a-1`) và `a(p-1)`.

Ghi chú : Trong trường hợp không cân bằng, số phần tử `p_k` của các nhóm không giống nhau, các ma trận `mb(H)` và `mb(E)` được xác định như sau:

`mb(H)= sum_(k=1)^a p_k (bar mb(x)_k - bar mb(x)) (bar mb(x)_k - bar mb(x))^T `	(21)
`mb(E)=sum_(k=1)^a sum_(j=1)^(p_k) (mb(x)_(kj) - bar mb(x)_k) (mb(x)_(kj) - bar mb(x)_k)^T`	(22)

Sau khi xác định được `mb(H)` và `mb(E)`, ta sử dụng kiểm định dùng cho phân tích phương sai dữ liệu đa biến. Ta có một số cách khác nhau, nhưng cách thực hiện gần tương tự nhau: tính số thống kê của mẫu, so sánh với giá trị tới hạn, từ đó rút ra kết luận. Trong phần này, chúng tôi giới thiệu hai kiểm định là kiểm định Wilks và kiểm định Pillai.

Kiểm định Wilks

Trong kiểm định Wilks, còn gọi là kiểm định Wilks `Lambda`, số thống kê `Lambda` được xác định bằng công thức:

`Lambda=(|mb(E)|)/(|mb(E)+mb(H)|)`

(23)

Sau đó so sánh giá trị `Lambda_o` của mẫu với giá trị tới hạn `Lambda`* (từ bảng Wilks). Nếu `Lambda_o < Lambda`*, ta bác bỏ Ho, nếu không, ta chấp nhận Ho. Từ đó ta có kết luận về vectơ trung bình của các nhóm.

Kiểm định Pillai

Trong kiểm định Pillai, số thống kê `V_s` được xác định bằng công thức:

`V_s=tr[ (mb(E)+mb(H))^(-1) mb(H) ]`(24)

Sau đó so sánh giá trị `V_s0` của mẫu với giá trị tới hạn `V_s`* (từ bảng Pillai). Nếu `V_(so) > V_s`*, ta bác bỏ Ho, nếu không, ta chấp nhận Ho.

Khi xác định giá trị tới hạn `V_s`* từ bảng Pillai, ta cần sử dụng ba thông số sau:

`s=min(df_H,n)`	(25)
`m=(|df_H-n|-1)/2`	(26)
`N=(df_E-n-1)/2`	(27)

Sau khi chuyển dữ liệu vào R với tên xso, ta sử dụng lệnh str để xem qua cấu trúc của bảng dữ liệu, ta có kết quả sau:

> str(xso)
'data.frame':	150 obs. of  5 variables:
 $ epoch: Factor w/ 5 levels "1850BC","200BC",..: 4 4 4 4 4 4 4 4 4 4 ...
 $ mb   : int  131 125 131 119 136 138 139 125 131 134 ...
 $ bh   : int  138 131 132 132 143 137 130 136 134 134 ...
 $ bl   : int  89 92 99 96 100 89 108 93 102 99 ...
 $ nh   : int  49 48 50 44 54 56 48 48 51 51 ...

Ta thấy dữ liệu gồm 150 phần tử (150 xương sọ) và 5 biến gồm 1 biến độc lập, epoch, và 4 biến phụ thuộc mb (maximum breaths), bh (basibregmatic heights), bl (basialiveolar length) và nh (nasal heights) là các kích thước của xương sọ.

Ta sử dụng lệnh manova của R để phân tích phương sai cho trường hợp nhiều biến phụ thuộc. Khi sử dụng lệnh này, ta cần lưu ý một số điểm sau:

• Tương tự như các lệnh lm và aov, mối quan hệ giữa biến phụ thuộc và biến độc lập được đưa vào đối số dưới dạng     (biến phụ thuộc) ~ (biến độc lập)
• Các biến phụ thuộc được liên kết lại bằng lệnh cbind().

Vậy trong thí dụ này ta có thể sử dụng lệnh sau:

 kq <- manova(cbind(mb, bh, bl, nh) ~ epoch, data = xso)

Kết quả xử lý được lưu vào biến kq. Đó là một danh sách với 13 thành phần (dùng lệnh names):

> names(kq)
 [1] "coefficients"  "residuals"     "effects"       "rank"          "fitted.values"
 [6] "assign"        "qr"            "df.residual"   "contrasts"     "xlevels"      
[11] "call"          "terms"         "model"

Để xem xét kết quả ta dùng lệnh summary. Trong trường hợp phân tích phương sai đa biến, ta có thể dùng đối số test để sử dụng cách kiểm định theo yêu cầu. Đối số này có 4 giá trị: "Pillai" (mặc định), "Roy", "Wilks", và "Hotelling-Lawley". Thí dụ:

> summary(kq)
           Df  Pillai approx F num Df den Df    Pr(>F)    
epoch       4 0.35331    3.512     16    580 4.675e-06 ***
Residuals 145                                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

hay :

> summary(kq, test = "Wilks")
           Df   Wilks approx F num Df den Df   Pr(>F)    
epoch       4 0.66359   3.9009     16 434.45 7.01e-07 ***
Residuals 145                                            
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Dựa vào giá trị `p` (Pr(>F)) ở độ tin cậy 95%, ta thấy cả hai kiểm định Pillai và Wilks đều có kết quả tương tự là bác bỏ Ho, nghĩa là kích thước xương sọ của các thời kỳ không hoàn toàn giống nhau.

Để có thể xem xét chi tiết hơn sự khác biệt của từng biến độc lập, ta sử dụng lệnh summary.aov. Thí dụ:

> summary.aov(kq)
 Response mb :
             Df  Sum Sq Mean Sq F value    Pr(>F)    
epoch         4  502.83 125.707  5.9546 0.0001826 ***
Residuals   145 3061.07  21.111                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response bh :
             Df Sum Sq Mean Sq F value  Pr(>F)  
epoch         4  229.9  57.477  2.4474 0.04897 *
Residuals   145 3405.3  23.485                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response bl :
             Df Sum Sq Mean Sq F value    Pr(>F)    
epoch         4  803.3 200.823  8.3057 4.636e-06 ***
Residuals   145 3506.0  24.179                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response nh :
             Df Sum Sq Mean Sq F value Pr(>F)
epoch         4   61.2  15.300   1.507 0.2032
Residuals   145 1472.1  10.153

Dựa vào giá trị `p` của kết quả phân tích này với độ tin cậy 95%, ta thấy xương sọ giữa các thời kỳ có sự khác biệt về các kích thước mb, bh và bl; riêng sự khác biệt của nh là không có ý nghĩa.

Ta có thể sử dụng lệnh manova để so sánh hai vectơ trung bình bằng cách dùng đối số subset để liệt kê các nhóm muốn so sánh. Thí dụ ta muốn so sánh kích thước xương sọ tại hai thời kỳ 4000BC và AD150 thì ta dùng lệnh sau:

 kq2 <- manova(cbind(mb, bh, bl, nh) ~ epoch, data = xso,
             subset = epoch %in% c("4000BC", "AD150"))

Kết quả thu được là :

> summary(kq2)
          Df  Pillai approx F num Df den Df    Pr(>F)    
epoch      1 0.36182   7.7956      4     55 4.736e-05 ***
Residuals 58                                             
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1