Trong các phần về phân tích phương sai mà ta vừa xem xét, mỗi chủ thể chỉ có thể tham gia vào một nhóm, mỗi cây trong vườn chỉ từ cây con do một công ty cung cấp, chịu tác động của một chế độ tưới. Nhưng trong nhiều trường hợp, cùng một chủ thể lại tham gia vào nhiều nhóm.

Lấy thí dụ một công ty C muốn so sánh ảnh hưởng của 4 phương pháp trang trí cửa hàng (ký hiệu `A_1` đến `A_4`) đến lượng sản phẩm bán ra (`Y`). Để thử nghiệm, công ty C sử dụng 5 đại lý (ký hiệu `S_1` đến `S_5`). Kết quả của đợt thử nghiệm này được trình bày trên Bảng 1.

Qua thí dụ trên, ta thấy mỗi đại lý đều thử nghiệm cả 4 phương pháp trang trí, vì thể ta chỉ cần 5 đại lý thay vì 20 đại lý như trường hợp không lặp. Ngoài ra vì mỗi đại lý đều thử nghiệm cả 4 phương pháp trang trí nên kết quả so sánh sự khác biệt giữa các phương pháp trang trí sẽ hợp lý hơn.

Trong trường hợp tổng quát, bảng dữ liệu của sự phân nhóm bởi 1 yếu tố với chủ thể lặp và bố trí cân bằng có thể được trình bày như trên Bảng 2.

Bảng 1 Giá trị của đáp ứng `Y` cho `a` nhóm với chủ thể lặp

		Yếu tố A						Trung bình `bar y_S`
		Nhóm `A_1`	Nhóm `A_2`	. . .	Nhóm `A_i`	. . .	Nhóm `A_a`
Chủ thể	`S_1`	`y_(11)`	`y_(21)`	. . .	`y_(i1)`	. . .	`y_(a1)`	`bar y_(S1)`
	`S_2`	`y_(12)`	`y_(22)`	. . .	`y_(i2)`	. . .	`y_(a2)`	`bar y_(S2)`
	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	`S_j`	`y_(1j)`	`y_(2j)`	. . .	`y_(ij)`	. . .	`y_(aj)`	`bar y_(Sj)`
	. . .	. . .	. . .	. . .	. . .	. . .	. . .
	`S_n`	`y_(1n)`	`y_(2n)`		`y_(i n)`	. . .	`y_(an)`	`bar y_(Sn)`
Trung bình `bar y_A`		`bar y_(A1)`	`bar y_(A2)`	. . .	`bar y_(Ai)`	. . .	`bar y_(Aa)`	`bar y`

Trong Bảng 2, `y_(ij)` là đáp ứng `Y` của chủ thể `S_j` tương ứng với nhóm `A_i`, `bar y_(Ai)` là trung bình của tất cả `n` chủ thể trong nhóm `A_i`, `bar y_(Sj)` là trung bình của chủ thể `S_j` trong tất cả `a` nhóm, `bar y` là trung bình của `an` số liệu.

 

Theo lý thuyết của mô hình tuyến tính, người ta xem kết quả `y_(ij)` được thể hiện dưới dạng:

`y_(ij)=mu+tau_i+e_(ij)`(26)

trong đó `mu` là trung bình của `an` giá trị `y_(ij)` ; `tau_i` là tác động của nhóm `A_i`; `beta_j` là ảnh hưởng của chủ thể `S_j`; và `e_(ij)` là sai lệch.

Cũng như phân tích phương sai trong các trường hợp khác, ta lần lượt thực hiện các bước sau:

• tính các tổng phương sai `SS`,
• tính các trung bình phương sai `MS=(SS)/(df)` (`df` là độ tự do),
• tính (các) giá trị `F_o` rồi so sánh với giá trị tới hạn `F`* (hoặc tính giá trị `p`),
• kết luận.

Trước hết, tính tổng phương sai chung `SS_T` theo công thức:

`SS_T=sum_(i=1)^a sum_(j=1)^n (y_(ij)-bar y)^2`

(27)

Tổng phương sai chung `SS_T` có độ tự do là `an-1` và gồm hai thành phần:

• tổng phương sai `SS_S` thể hiện mức độ khác biệt giữa các chủ thể và được xác định bằng công thức:

  `SS_S=a sum_(j=1)^n (bar y_(Sj)-bar y)^2`

  (28)

  tổng phương sai này có độ tự do là `n-1`,
• tổng phương sai `SS_(wS)` (within subject) thể hiện sự thay đổi trong từng chủ thể giữa các nhóm và được xác định bằng công thức:

  `SS_(wS)=sum_(j=1)^n sum_(i=1)^a (y_(ij)-bar y_(Sj))^2`

  (29)

  tổng phương sai này có độ tự do là `n(a-1)`,

Người ta chứng minh được:

`SS_T=SS_S+SS_(wS)`(30)

Tổng phương sai `SS_(wS)` lại bao gồm hai thành phần:

• tổng phương sai nhóm (hay yếu tố) `SS_A` thể hiện mức độ khác biệt giữa các nhóm và được xác định bằng công thức:

  `SS_A=n sum_(i=1)^a (bar y_i-bar y)^2`

  (31)

  tổng phương sai này có độ tự do là `a-1`,
• tổng phương sai sai lệch `SS_E` thể hiện khác biệt giữa các nhóm trong một chủ thể; tổng phương sai này có độ tự do là `(n-1)(a-1)`,

Ta cũng có :

`SS_(wS)=SS_A+SS_E`(32)

Vậy :

`SS_T=SS_S+SS_A+SS_E`(33)

Từ đây, chúng ta có thể tiếp tục tính các trung bình phương sai `MS`, giá trị `F_o`, ... như các trường hợp phân tích phương sai khác.

 

Ta sẽ tiếp tục với thí dụ về khảo sát ảnh hưởng của phương pháp trang trí cửa hàng (yếu tố A) đến lượng sản phẩm bán được (`Y`). Ta có bảng số liệu sau:

Giả sử ta chỉ quan tâm tác động của yếu tố A đến đáp ứng `Y`, vì vậy ta chỉ cần `MS_A` và `MS_E`, hay chỉ cần tính `SS_(wS)` và `SS_A` mà không cần tính các tổng phương sai khác.

Từ các công thức (15) và (17) ta có:

`SS_(wS)=sum_(j=1)^n sum_(i=1)^a (y_(ij)-bar y_(Sj))^2= [(4-6)^2+(5-6)^2+...+(5-3)^2+(4-3)^2]=58`

`SS_A=n sum_(i=1)^a (bar y_(Ai)-bar y)^2=5xx[(2-4)^2+(3-4)^2+(5-4)^2+(6-4)^2]=50`

Từ (18) ta suy ra :

  `SS_E=SS_(wS)-SS_A=58-50=8`

Với các độ tự do `df_A=3` và `df_E=12` ta có :

  `MS_A=(SS_A)/(df_A)=50/3=16,667`

  `MS_E=(SS_E)/(df_E)=8/12=0,6667`

Vậy :   `F_(oA)=(MS_A)/(MS_E)=(16,667)/(0,6667)=25`

Với độ tin cậy 95% thì :   `F`*`=F_(0,05, 3, 12)=3,490`

Vì `F_(oA) > F`* nên ta kết luận yếu tố A có tác động đến `Y`, phương pháp trang trí cửa hàng có ảnh hưởng đến lượng sản phẩm bán được.

Sử dụng R

Trước hết, ta chuẩn bị dữ liệu chưa trong tập tin pp_trang_tri.csv, trong đó có ba biến là PP_TT, Dai_Ly, và San_Pham. Ta chuyển tập tin này vào R để tạo bảng dữ liệu có tên là sp.

Khi phân tích phương sai với chủ thể lặp, ta vẫn sử dụng hàm aov nhưng đưa thêm đối số Error(Dai_Ly/PP_TT) vào như đoạn lệnh sau:

 kql <- aov(San_Pham ~ PP_TT + Error(Dai_Ly/PP_TT), data = sp)
 summary(kql)

Ta thu được kết quả sau :

> summary(kql)

Error: Dai_Ly
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals  4     24       6               

Error: Dai_Ly:PP_TT
          Df Sum Sq Mean Sq F value  Pr(>F)    
PP_TT      3     50  16.667      25 1.9e-05 ***
Residuals 12      8   0.667                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ta thu được kết quả tương tự như đã thực hiện ở trên.

 

Cũng như khi phân tích phương sai các trường hợp thông thường, sau khi có kết luận bác bỏ giả thuyết không Ho, ta cần tiếp tục so sánh sự khác biệt giữa các nhóm. Khi tính toán thủ công, quy trình thực hiện tương tự như khi chủ thể không lặp, bạn có thể tham khảo thêm trong phần "thí nghiệm với yếu tố lặp".

Nếu ta sử dụng R, quá trình sẽ tương đối phức tạp hơn, ta cần sử dụng đến các phụ kiện nlme (linear and non linear mixed effect model) và multcomp (multiple comparison) theo trình tự như sau:

• Sử dụng hàm lme của phụ kiện nlme để xác lập mô hình cho dữ liệu, trong đó, ta đưa chủ thể vào đối số. Dùng một biến (thí dụ kq) để chứa kết quả quá trình xử lý này.
• Sử dụng hàm glht của phụ kiện multcomp để so sánh các nhóm từ kết quả kq ở trên.

Để minh họa rõ hơn các điểm trên, ta so sánh 4 phương pháp trang trí trong thí dụ trên. Trước hết xác lập mô hình bằng lệnh lme rồi lưu kết quả vào biến kqm bằng đoạn lệnh sau:

 library(nlme)
 kqm <- lme(San_Pham ~ PP_TT, random = ~1|Dai_Ly/PP_TT, data = sp)

Ta thấy hàm lme sử dụng dạng thức phân tích tương tự các hàm lm hay aov, nhưng khai báo thêm hai biến yếu tố và chủ thể trong đối số random.

Tiếp tục dùng hàm glht để so sánh các nhóm bằng cách sử dụng kết quả lưu trong biến kqm rồi dùng lệnh summary để xem kết quả so sánh theo đoạn lệnh sau:

 library(multcomp)
 kcomp <- glht(kqm, linfct = mcp(PP_TT="Tukey"))
 summary(kcomp)

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

> summary(kcomp)

	 Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: lme.formula(fixed = San_Pham ~ PP_TT, data = sp, random = ~1 | 
    Dai_Ly/PP_TT)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
A2 - A1 == 0   1.0000     0.5164   1.936    0.213    
A3 - A1 == 0   3.0000     0.5164   5.809   <0.001 ***
A4 - A1 == 0   4.0000     0.5164   7.746   <0.001 ***
A3 - A2 == 0   2.0000     0.5164   3.873   <0.001 ***
A4 - A2 == 0   3.0000     0.5164   5.809   <0.001 ***
A4 - A3 == 0   1.0000     0.5164   1.936    0.213    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

Ta thấy có 6 sự so sánh giữa 4 phương pháp trang trí đã được thực hiện, trong đó các sự khác biệt giữa `A_1` với `A_2` và `A_3` với `A_4` không có ý nghĩa thống kê, còn 4 sự khác biệt khác đều có ý nghĩa.