Tương tự như các trường hợp đã khảo sát, ta chuyển tập tin chieu-cao-cay.csv vào R và có bảng dữ liệu ccao. Trong bảng này ba biến mà ta quan tâm là đáp ứng Chieu_Cao và hai yếu tố CD_Tuoi và Cung_Cap. Ta sử dụng hàm aov để khảo sát tác động của chế độ tưới và nhà cung cấp bằng đoạn lệnh sau:

 kq2 <- aov(Chieu_Cao ~ Cung_Cap*CD_Tuoi, data = ccao)
 summary(kq2)

Kết quả thu được như sau :

> summary(kq2)
                 Df Sum Sq Mean Sq F value   Pr(>F)    
Cung_Cap          1   4417    4417  39.102 1.83e-06 ***
CD_Tuoi           2  10574    5287  46.809 5.21e-09 ***
Cung_Cap:CD_Tuoi  2    994     497   4.402   0.0235 *  
Residuals        24   2711     113                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Ta thấy kết quả trên có dạng gần tương tự như Bảng 3. Dựa vào giá trị `p` (cột "Pr(>F)") ta có thể kết luận rằng cả nhà cung cấp và chế độ tưới đều có ảnh hưởng đến chiều cao cây ở độ tin cậy 95%, ngoài ra giữa hai yếu tố này có sự tương tác.

Cũng như trường hợp một yếu tố, ta có thể phân tích sâu hơn bằng cách dùng lệnh TukeyHSD để so sánh các nhóm với nhau bằng đoạn lệnh:

 kt2 <- TukeyHSD(kq2, ordered = TRUE)
 kt2

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

> kt2
  Tukey multiple comparisons of means
    95% family-wise confidence level
    factor levels have been ordered

Fit: aov(formula = Chieu_Cao ~ Cung_Cap * CD_Tuoi, data = ccao)

$Cung_Cap
        diff      lwr      upr   p adj
A-B 24.26667 16.25725 32.27608 1.8e-06

$CD_Tuoi
                 diff       lwr      upr     p adj
Trung_Binh-It    20.5  8.630659 32.36934 0.0006734
Nhieu-It         45.9 34.030659 57.76934 0.0000000
Nhieu-Trung_Binh 25.4 13.530659 37.26934 0.0000503

$`Cung_Cap:CD_Tuoi`
                          diff        lwr      upr     p adj
A:It-B:It                 30.0   9.217255 50.78275 0.0020055
B:Trung_Binh-B:It         31.4  10.617255 52.18275 0.0012006
A:Trung_Binh-B:It         39.6  18.817255 60.38275 0.0000595
B:Nhieu-B:It              43.6  22.817255 64.38275 0.0000142
A:Nhieu-B:It              78.2  57.417255 98.98275 0.0000000
B:Trung_Binh-A:It          1.4 -19.382745 22.18275 0.9999382
A:Trung_Binh-A:It          9.6 -11.182745 30.38275 0.7102448
B:Nhieu-A:It              13.6  -7.182745 34.38275 0.3590393
A:Nhieu-A:It              48.2  27.417255 68.98275 0.0000029
A:Trung_Binh-B:Trung_Binh  8.2 -12.582745 28.98275 0.8229452
B:Nhieu-B:Trung_Binh      12.2  -8.582745 32.98275 0.4753775
A:Nhieu-B:Trung_Binh      46.8  26.017255 67.58275 0.0000046
B:Nhieu-A:Trung_Binh       4.0 -16.782745 24.78275 0.9903911
A:Nhieu-A:Trung_Binh      38.6  17.817255 59.38275 0.0000855
A:Nhieu-B:Nhieu           34.6  13.817255 55.38275 0.0003699

Bảng trên cho ta một số kết quả sau :

• Trong phần "$Cung_Cap", ta thấy khoảng ước lượng cho độ chênh lệch giữa hai nhà cung cấp là [16,25725, 32,27608] không chứa 0, hay `p = 0,0000018`. Điều này chứng tỏ có sự khác biệt giữa hai nhà cung cấp.
• Tương tự, trong phần "CD_Tuoi" các kết quả cũng cho thấy chế độ tưới có ảnh hưởng đến chiều cao cây.
• Trong phần "Cung_Cap:CD_Tuoi", ta thấy có 15 sự so sánh giữa 6 nhóm. Trong đó có 6 trường hợp mà sự khác biệt giữa các nhóm không có ý nghĩa thống kê ở khoảng tin cậy 95% (khi khoảng tin cậy chứa 0 hay `p< 0,05`)

Sự khác biệt giữa 6 nhóm trên có thể trình bày ở dạng biểu đồ khi ta dùng đoạn lệnh sau:

 par(mar = c(5, 12, 4, 2))
 plot(kt2, las = 1, col = "red")
 abline(v = 0, lty = "64", col = "blue", lwd = 2)

Trong đoạn lệnh trên, câu lệnh thứ hai là quan trọng nhất, dùng để vẽ. Hai câu lệnh còn lại điều chỉnh một số chi tiết của biểu đồ để tăng tính mỹ thuật. Kết quả được trình bày trên Hình 1.

Hình 1 So sánh tác động của các nhóm đến chiều cao cây

Trên Hình 1, trục hoành thể hiện chênh lệch của trung bình các nhóm; đường "0" thẳng đứng màu xanh tương ứng với độ chênh lệch bằng 0. Trên mỗi dòng của phần chú giải ở bên trái của Hình 1, nhóm bên trái có chiều cao cây lớn hơn. Như vậy cây thuộc nhóm "B-It" có chiều cao thấp nhất, nhóm "A-Nhieu" có chiều cao lớn nhất. Khoảng ước lượng cho hiệu số của trung bình chiều cao giữa các nhóm được biểu thị bằng các đoạn màu đỏ. Trên Hình 1, ta thấy có 6 đoạn màu đỏ cắt đường "0" màu xanh, như vậy có 6 cặp nhóm mà sự khác biệt về chiều cao không có ý nghĩa thống kê (thí dụ giữa hai nhóm A-Trung_Binh và A-It).

Ta có thể dùng hàm plot phối hợp với hàm Effect của phụ kiện effects để có một nhận định khái quát về tác động của hai yếu tố đến chiều cao cây bằng đoạn lệnh sau:

 library(effects)
 plot(Effect(c("Cung_Cap", "CD_Tuoi"), kq2))

Ta có Hình 2.

Hình 2 Tác động của các yếu tố đến chiều cao cây