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

DATA

Tần số xuất hiện của các vị trí nguồn gốc

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

Thực hiện MANOVA

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.

Đo kích thước của hiệu ứng

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.

Triển khai Phân tích Phân biệt Tuyến tính LDA

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()