MANOVA là viết tắt của Multivariate ANOVA hoặc Multivariate Analysis
Of Variance. Đó là một phần mở rộng của ANOVA thông thường.Ý tưởng chung
là giống nhau (sử dụng thử nghiệm để xác định ảnh hưởng của các biến độc
lập đối với biến phụ thuộc trong nghiên cứu hồi quy), nhưng kiểm tra
MANOVA phải bao gồm ít nhất hai biến phụ thuộc để phân tích sự khác biệt
giữa nhiều nhóm (nhân tố) của biến độc lập.
MANOVA trong R sử dụng kiểm tra Pillai’s Trace cho các tính toán, sau đó
được chuyển đổi thành thống kê F khi muốn kiểm tra mức độ quan trọng của
sự khác biệt trung bình của nhóm.
Giá trị rỗng và giả thuyết thay thế:
- H0: Các vectơ trung bình của nhóm giống nhau cho tất cả các nhóm hoặc
chúng không khác nhau đáng kể.
- H1: Ít nhất một trong các vectơ trung bình của nhóm khác với các vectơ
còn lại.
MANOVA trong R sẽ không cho biết nhóm nào khác với nhóm còn lại, nhưng điều đó dễ xác định thông qua post-hoc test, sử dụng Phân tích phân biệt tuyến tính (LDA) để trả lời câu hỏi này sau.
VD:
Nghiên cứu kiểm tra xem khu vực sống (các tỉnh thành khác nhau) và các
giống tiêu khác nhau có ảnh hưởng đến hàm lượng độ ẩm, tro tổng,
piperin, tinh dầu.
Biến độc lập là các tỉnh thành khác nhau
Biến phụ thuộc hàm lượng độ ẩm, tro tổng, piperin, tinh dầu
ViTri.freq <- table(hoalyTieu$ViTri)
ViTri.freq
##
## BRVT DakLak DakNong DongNai GiaLai QuangTri
## 6 5 7 9 9 30
barplot(ViTri.freq,
horiz = F,
col = rainbow(length(ViTri.freq)))
Các biến phụ thuộc
box_DoAm2 <- ggplot(hoalyTieu , aes(ViTri, DoAm, fill = ViTri)) + geom_boxplot(show.legend = F)
box_TroTong2 <- ggplot(hoalyTieu , aes(ViTri, TroTong, fill = ViTri)) + geom_boxplot(show.legend = F)
box_Piperine2 <- ggplot(hoalyTieu , aes(ViTri, Piperine, fill = ViTri)) + geom_boxplot(show.legend = F)
box_TinhDau2 <- ggplot(hoalyTieu , aes(ViTri, TinhDau, fill = ViTri)) + geom_boxplot(show.legend = F)
grid.arrange(box_DoAm2, box_TroTong2, box_Piperine2, box_TinhDau2, ncol = 2, nrow = 2)
Nhìn có vẻ mẫu tiêu ở Quảng Trị các chỉ tiêu hóa lý khác biệt với các mẫu còn lại
by( hoalyTieu$DoAm, hoalyTieu$ViTri, summary)
## hoalyTieu$ViTri: BRVT
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.932 10.490 11.123 11.218 11.708 13.947
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakLak
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.96 13.57 14.53 13.87 14.58 14.72
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakNong
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 10.85 11.46 12.28 12.43 13.08 14.77
## ------------------------------------------------------------
## hoalyTieu$ViTri: DongNai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.598 10.288 11.333 11.192 11.869 13.844
## ------------------------------------------------------------
## hoalyTieu$ViTri: GiaLai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.94 14.13 14.64 14.40 14.80 15.21
## ------------------------------------------------------------
## hoalyTieu$ViTri: QuangTri
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.01 11.70 12.11 12.57 13.28 15.67
by( hoalyTieu$TroTong, hoalyTieu$ViTri, summary)
## hoalyTieu$ViTri: BRVT
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.011 3.488 3.617 3.799 4.297 4.577
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakLak
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.374 3.443 3.628 3.802 4.208 4.356
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakNong
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.088 3.549 3.654 3.637 3.772 4.077
## ------------------------------------------------------------
## hoalyTieu$ViTri: DongNai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.050 3.257 3.469 3.514 3.615 4.339
## ------------------------------------------------------------
## hoalyTieu$ViTri: GiaLai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.086 3.327 3.535 3.550 3.713 4.098
## ------------------------------------------------------------
## hoalyTieu$ViTri: QuangTri
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.381 3.887 4.160 4.113 4.322 4.644
by( hoalyTieu$Piperine, hoalyTieu$ViTri, summary)
## hoalyTieu$ViTri: BRVT
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.504 3.786 4.200 3.957 4.459 4.647
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakLak
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.315 4.591 4.709 4.672 4.716 5.028
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakNong
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.663 4.780 4.997 5.132 5.480 5.748
## ------------------------------------------------------------
## hoalyTieu$ViTri: DongNai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.415 3.670 3.833 3.823 4.024 5.177
## ------------------------------------------------------------
## hoalyTieu$ViTri: GiaLai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.412 4.572 4.903 4.662 4.996 5.316
## ------------------------------------------------------------
## hoalyTieu$ViTri: QuangTri
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.203 4.629 4.842 4.998 5.453 5.980
by( hoalyTieu$TinhDau, hoalyTieu$ViTri, summary)
## hoalyTieu$ViTri: BRVT
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.400 2.450 2.650 2.683 2.775 3.200
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakLak
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.4 2.6 2.6 2.6 2.6 2.8
## ------------------------------------------------------------
## hoalyTieu$ViTri: DakNong
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.200 3.200 3.200 3.229 3.500 3.800
## ------------------------------------------------------------
## hoalyTieu$ViTri: DongNai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.400 2.800 3.000 3.022 3.200 4.000
## ------------------------------------------------------------
## hoalyTieu$ViTri: GiaLai
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.800 2.800 3.200 3.156 3.600 4.000
## ------------------------------------------------------------
## hoalyTieu$ViTri: QuangTri
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.60 2.80 3.00 3.07 3.20 3.90
dependent_vars <- cbind(hoalyTieu$DoAm, hoalyTieu$TroTong, hoalyTieu$Piperine, hoalyTieu$TinhDau)
independent_var <- hoalyTieu$ViTri
manova_model <- manova(dependent_vars ~ independent_var, data = hoalyTieu)
summary(manova_model)
## Df Pillai approx F num Df den Df Pr(>F)
## independent_var 5 1.1509 4.8472 20 240 7.026e-10 ***
## Residuals 60
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Tóm tắt thử nghiệm MANOVA trong R
Mô hình MANOVA kết quả báo cáo thống kê thử nghiệm Pillai là 1.1509 và p-value dưới 0,05 , có nghĩa là có thể bác bỏ giả thuyết rỗng một cách an toàn để ủng hộ giả thuyết thay thế - ít nhất một vectơ trung bình của nhóm khác với phần còn lại.
Pillai càng gần với 1, bằng chứng cho thấy biến giải thích có ảnh hưởng có ý nghĩa thống kê đến giá trị của các biến phản hồi càng mạnh
Thống kê Pillai được sử dụng để xác định ý nghĩa của sự khác biệt của các vectơ trung bình bằng cách so sánh nó với giá trị tới hạn.
Giá trị Bình phương Eta một phần (Partial Eta Squared) đo lường tác động của biến độc lập đối với các biến phụ thuộc. Nếu giá trị là 0,14 hoặc lớn hơn, có thể nói kích thước hiệu ứng lớn.
eta_squared(manova_model)
## # Effect Size for ANOVA (Type I)
##
## Parameter | Eta2 (partial) | 95% CI
## -----------------------------------------------
## independent_var | 0.29 | [0.16, 1.00]
##
## - One-sided CIs: upper bound fixed at [1.00].
Giá trị Bình phương Eta một phần cho thử nghiệm MANOVA
Giá trị là 0,29, có nghĩa là kích thước hiệu ứng lớn. Đây là một cách tốt để kiểm tra lại kết quả tóm tắt của thử nghiệm MANOVA, nhưng làm thế nào có thể thực sự biết vectơ trung bình của nhóm nào khác với các vectơ còn lại? Đây là lúc sử dụng post-hoc test ̣(bài kiểm tra sau giờ học)
Sử dụng Phân tích phân biệt tuyến tính (LDA), tìm kiếm sự kết hợp tuyến tính của các đối tượng đặc điểm phân tách tốt nhất hai hoặc nhiều nhóm.
hoaly_lda <- lda(independent_var ~ dependent_vars, CV = F)
hoaly_lda
## Call:
## lda(independent_var ~ dependent_vars, CV = F)
##
## Prior probabilities of groups:
## BRVT DakLak DakNong DongNai GiaLai QuangTri
## 0.09090909 0.07575758 0.10606061 0.13636364 0.13636364 0.45454545
##
## Group means:
## dependent_vars1 dependent_vars2 dependent_vars3 dependent_vars4
## BRVT 11.21782 3.798650 3.957417 2.683333
## DakLak 13.87011 3.801897 4.671953 2.600000
## DakNong 12.42780 3.637114 5.132358 3.228571
## DongNai 11.19231 3.513844 3.823422 3.022222
## GiaLai 14.39786 3.550371 4.661612 3.155556
## QuangTri 12.56708 4.112806 4.997888 3.070000
##
## Coefficients of linear discriminants:
## LD1 LD2 LD3 LD4
## dependent_vars1 0.53285105 -0.3615833 0.3604275 -0.2673852
## dependent_vars2 0.82292084 2.0397214 0.9229628 -1.5788975
## dependent_vars3 1.20822134 0.1221841 -0.7097773 1.2275978
## dependent_vars4 -0.04013983 -0.5606075 -1.7225504 -1.8430943
##
## Proportion of trace:
## LD1 LD2 LD3 LD4
## 0.5892 0.2442 0.1265 0.0401
Kết quả Phân tích Phân biệt Tuyến tính
lda_df <- data.frame(
ViTri = hoalyTieu[, "ViTri"],
lda = predict(hoaly_lda)$x
)
lda_df
## ViTri lda.LD1 lda.LD2 lda.LD3 lda.LD4
## 1 DakLak 0.28138526 -1.47223031 0.98295526 0.066070934
## 2 DakLak 1.34859935 0.55066509 1.95737343 -0.655140459
## 3 DakLak -0.46862105 0.01194365 0.27054198 1.334421848
## 4 QuangTri 0.14666951 2.24228761 0.51680782 -0.287876849
## 5 QuangTri -0.27438705 1.36012617 0.03607972 -0.373859980
## 6 QuangTri -0.06827178 -0.19202593 -1.74340226 -1.392790635
## 7 QuangTri -0.21777864 0.31412223 -0.75614330 0.122012949
## 8 QuangTri 0.44050858 1.25724077 -0.77072883 0.923457848
## 9 QuangTri -0.83064433 0.52926386 -0.40309352 -0.562429272
## 10 QuangTri 2.05477102 1.51738291 0.31543978 0.636608443
## 11 QuangTri 0.89369792 0.64687753 0.90548352 -0.224357863
## 12 QuangTri 1.46472651 0.68570912 -0.02632234 -0.556273426
## 13 QuangTri 1.50885532 1.18656451 -0.67604714 1.072900051
## 14 QuangTri 1.73419140 -1.14175587 1.32103626 -0.124960459
## 15 QuangTri 0.91251283 -1.53137879 0.70455226 -1.114358237
## 16 QuangTri 0.52738908 1.62615885 -0.19427408 -0.797816031
## 17 QuangTri 1.35698003 1.82403480 0.06225836 0.732169437
## 18 QuangTri 0.46522444 -1.11536276 -0.17032296 -0.140730818
## 19 DakNong -0.06379719 0.41535638 -1.58864287 1.265278659
## 20 DakNong 1.36846547 -1.25560795 0.07654933 -0.220494214
## 21 DakNong 0.52688364 0.13646172 -1.18456302 0.931362719
## 22 DakNong -0.69674059 0.57719841 -1.67326575 -1.406906649
## 23 DakNong -0.08845049 0.07535985 0.91914515 2.256564900
## 24 DakNong -0.24707643 -2.29826950 -1.73089127 -0.479919954
## 25 GiaLai 2.13000751 -0.31956150 0.26774802 -0.535739011
## 26 GiaLai -0.63276949 -0.79879406 3.48421384 0.589233804
## 27 GiaLai 1.29980668 -2.52936519 -1.44528869 -1.296208519
## 28 GiaLai 1.82102618 -0.96357050 -0.87835352 -1.478796330
## 29 GiaLai 0.04449079 -2.23336735 -0.37766718 0.326181705
## 30 QuangTri -0.35832157 0.68032975 -0.03194849 -0.990905105
## 31 QuangTri 0.39705570 0.49125548 0.20773082 -1.084321281
## 32 QuangTri 1.62798018 1.39165465 -1.17995511 -0.857342096
## 33 QuangTri -0.25738086 0.82252860 -0.82289126 -1.655305167
## 34 QuangTri -1.27723745 0.06481162 -0.86230985 0.072608470
## 35 QuangTri 0.03373039 0.60508046 -0.29837526 -0.513931202
## 36 QuangTri -1.58303237 -0.44604563 -0.70585333 0.780672899
## 37 BRVT -1.78180325 0.69465133 0.29154901 1.674638403
## 38 BRVT -1.69687539 0.08481684 0.26741394 0.589096165
## 39 BRVT -0.94740627 -0.68198727 -0.80046689 0.260782595
## 40 DongNai -2.63598552 0.22331759 0.52100266 1.275978034
## 41 DongNai 0.40249405 -0.94515995 -2.21504155 -0.875330854
## 42 DongNai -2.51639664 -0.45035724 -0.87258931 0.781318747
## 43 DongNai -2.70284063 -1.41992072 0.36004172 0.173247736
## 44 DongNai -1.18618451 1.16187479 0.37647041 -1.871432036
## 45 BRVT -1.14152428 -1.97783043 1.48617643 0.893342430
## 46 DongNai -2.38292567 -1.38755104 2.09154150 -2.111262325
## 47 DongNai -2.51639664 -0.45035724 -0.87258931 0.781318747
## 48 BRVT -0.26077390 2.27124549 0.86476647 -0.003677793
## 49 BRVT -4.02974169 2.53768348 1.22576604 -2.344059702
## 50 DongNai -3.60296852 0.62390551 -1.29446848 -0.049581622
## 51 DongNai -1.45713189 0.03460208 0.06209662 0.507555100
## 52 DakLak 1.85747701 0.23407132 1.57950059 0.063823622
## 53 DakLak 0.17912667 -0.97587430 0.96567138 1.664810817
## 54 QuangTri 1.00625990 -0.53908375 0.75582086 0.311388941
## 55 QuangTri 0.89235056 1.20316040 0.78883126 0.654765586
## 56 QuangTri 0.13028347 0.83930006 -0.98638498 1.248595705
## 57 QuangTri 2.34274024 1.13188014 0.17883911 0.991506899
## 58 QuangTri 0.47187382 0.67921236 0.70853588 -0.067568662
## 59 QuangTri 1.76048477 0.68420323 -0.11587991 -0.362152348
## 60 QuangTri 1.07085777 0.71770040 -0.65219553 1.054413342
## 61 QuangTri 0.90974308 -0.99392907 1.72205079 -0.908799275
## 62 DakNong 1.06617283 -0.73755056 -1.33014647 1.485219304
## 63 GiaLai 0.11473419 -1.48897452 1.49950814 -0.672455978
## 64 GiaLai 0.74623731 -1.54067624 -1.01802450 -0.538248857
## 65 GiaLai 0.96817832 -0.92044817 -0.02892943 -0.260668044
## 66 GiaLai -0.38050769 -1.32700319 -0.06644198 1.294354214
LDA dataset
Không có sự khác biệt đáng kể khi so sánh chỉ tiêu hóa lý của các mẫu tiêu khác vùng.
ggplot(lda_df) +
geom_point(aes(x = lda.LD1, y = lda.LD2, color = ViTri), size = 4) +
theme_classic()
ggplot(lda_df) +
geom_point(aes(x = lda.LD1, y = lda.LD3, color = ViTri), size = 4) +
theme_classic()
ggplot(lda_df) +
geom_point(aes(x = lda.LD1, y = lda.LD4, color = ViTri), size = 4) +
theme_classic()
ggplot(lda_df) +
geom_point(aes(x = lda.LD2, y = lda.LD3, color = ViTri), size = 4) +
theme_classic()
ggplot(lda_df) +
geom_point(aes(x = lda.LD2, y = lda.LD4, color = ViTri), size = 4) +
theme_classic()
ggplot(lda_df) +
geom_point(aes(x = lda.LD3, y = lda.LD4, color = ViTri), size = 4) +
theme_classic()