Bảng Phân Công Công Việc

Stt MSSV Tên Công Việc
1 21110310 Nguyễn Ngọc Huynh 4.9, 4.8
2 21110446 Chu Nguyễn Gia Khánh 4.7, 4.10
3 21110319 Nguyễn Thị Anh Tuyết 4.5, 4.5

Bài 4.5

Quan trắc trên 24 con rùa cái và 24 con rùa đực.

câu (a) Kiểm định sự bằng nhau của hai vectơ trung bình tổng thể, với \(\alpha = 0.01\).

câu (b) Nếu giải thuyết trong câu (a) bị bác bỏ, tìm tổ hợp tuyến tính của các trung bình thành phần có ý nghĩa quan trọng trong việc bác bỏ \(H_0\).

câu (c) Tìm các khyoanrg tin cậy đồng thời cho hiệu các trung bình thành phần. So sánh với các khoảng Bonferroni.

Mô tả dữ liệu bài 4.5

Dữ liệu trong bài toán 4.5 liên quan đến 24 con rùa cái và 24 con rùa đực.

Mỗi con rùa được mô tả bởi 3 biến liên tục đại diện cho các chỉ số:

  • Length: chiều dài
  • Width: chiều ngang
  • Height: chiều cao

Tóm tắt đặc trưng dữ liệu:

  • Số biến phụ thuộc: 3
  • Số chiều của dữ liệu: 3 chiều (mỗi con rùa có 3 chỉ số)
  • Loại phân tích: MANOVA 1_factor

Trước tiên ta tiến hành nạp dữ liệu.

# Đọc dữ liệu
data <- read.table("exo4-5.dat", header = FALSE)

head(data)
##    V1 V2 V3     V4
## 1  98 81 38 female
## 2 103 84 38 female
## 3 103 86 42 female
## 4 105 86 42 female
## 5 109 88 44 female
## 6 123 92 50 female
n = dim(data)[1]
p = dim(data)[2]

layout(matrix(c(1, 2, 3), nrow = 1, ncol = 3 ))

hist(data[,1], main = "Lenght")
hist(data[,2], main = "Width")
hist(data[,3], main = "Height")

Nhận xét

Dựa vào histogram của từng biến V1, V2, V3 ta thấy rằng dữ liệu có xu hướng lệch phải, nên ta cần phải thực hiện biến đổi để đưa dữ liệu về dạng gần chuẩn hơn. Ta sẽ thực hiện biến đổi Box-Cox để dảm bảo các giả thuyết khi áp dụng MANOVA được thỏa mãn.

Ta tiến hành chuyển đổi cột giới tính về dạng factor và khảo sát histogram theo từng giới tính.

data_45 <- data
data_45$V4 <- as.factor(data_45$V4)

male = data_45[(which(data_45$V4 == 'male')),]
female = data_45[(which(data_45$V4 == 'female')),]
layout(matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3))

hist(male[,1], main = "Length - Male")
hist(male[,2], main = "Width - Male")
hist(male[,3], main = "Height - Male")

hist(male[,1], main = "Length - Female")
hist(male[,2], main = "Width - Female")
hist(male[,3], main = "Height - Female")

Nhận xét

Ta thấy rằng histogram của từng biến theo giới tính có dạng gần chuẩn.

Để đảm bảo cho các giả thuyết của mô hình MANOVA, ta thực hiện biến đổi Box_Cox nhiều chiều.

library(car)
## Loading required package: carData
boxcox <- powerTransform(cbind(data$V1, data$V2, data$V3), family = "bcPower")

summary(boxcox)
## bcPower Transformations to Multinormality 
##    Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1    0.7273           1      -0.2417       1.6962
## Y2   -0.4340           0      -1.4804       0.6125
## Y3   -0.4810           0      -1.5965       0.6345
## 
## Likelihood ratio test that transformation parameters are equal to 0
##  (all log transformations)
##                                LRT df      pval
## LR test, lambda = (0 0 0) 13.58573  3 0.0035269
## 
## Likelihood ratio test that no transformations are needed
##                                LRT df       pval
## LR test, lambda = (1 1 1) 18.95005  3 0.00027998
df_transformed <- as.data.frame(bcPower(cbind(data$V1, data$V2, data$V3), coef(boxcox)))

data_45 <- data
data_45[,-4] <- df_transformed

male = data_45[(which(data_45$V4 == 'male')),]
female = data_45[(which(data_45$V4 == 'female')),]

Ta tiến hành kiểm tra mô hình có thỏa các điều kiện: tuân theo phân phối chuẩn nhiều chiều của từng tổng thể, tính đồng nhất của ma trận hiệp phương sai, tính độc lập của các tổng thể.

Kiểm tra các điều kiện

Kiểm tra mô hình có tuân theo phân phối chuẩn nhiều chiều của từng tổng thể

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:car':
## 
##     recode
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(rstatix)
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:stats':
## 
##     filter
data_45 %>%
  group_by(V4) %>%
  shapiro_test(V1, V2, V3) %>%
  arrange(variable)
## # A tibble: 6 × 4
##   V4     variable statistic     p
##   <chr>  <chr>        <dbl> <dbl>
## 1 female V1           0.951 0.279
## 2 male   V1           0.971 0.698
## 3 female V2           0.978 0.846
## 4 male   V2           0.979 0.886
## 5 female V3           0.952 0.303
## 6 male   V3           0.964 0.531
library(ggpubr)
## Loading required package: ggplot2
layout(matrix(c(1, 2, 3, 4, 5, 6), nrow = 2, ncol = 3))

# QQ plot Length của rùa

ggqqplot(data_45, "V1", facet.by = "V4", ylab = "V1", ggtheme = theme_bw())

# QQ plot Width của rùa

ggqqplot(data_45, "V2", facet.by = "V4", ylab = "V2", ggtheme = theme_bw())

# QQ plot Height của rùa

ggqqplot(data_45, "V3", facet.by = "V4", ylab = "V3", ggtheme = theme_bw())

data_45 %>%
  select(V1, V2, V3) %>%
  mshapiro_test()
## # A tibble: 1 × 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.958  0.0840
male %>%
  select(V1, V2, V3) %>%
  mshapiro_test()
## # A tibble: 1 × 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.917  0.0509
female %>%
  select(V1, V2, V3) %>%
  mshapiro_test()
## # A tibble: 1 × 2
##   statistic p.value
##       <dbl>   <dbl>
## 1     0.945   0.208

Giả thuyết kiểm định tổng thể có phân phối chuẩn

  • \(H_0\) (Giả thuyết không): \(V_i, i=1,2,3\) ứng với giới tính rùa đực, rùa cái được lấy từ tổng thể có phân phối chuẩn.

  • \(H_1\) (Giả thuyết đối): \(V_i, i=1,2,3\) ứng với giới tính rùa đực, rùa cái được lấy từ tổng thể không có phân phối chuẩn.

Nhận xét Ta nhận thấy rằng p-value của từng biến \(V_i\) ứng với mỗi giới tính lớn hơn 0.05. Ta chấp nhận \(H_0\) với độ tin cậy 95%.

Giả thuyết kiểm định tổng thể có phân phối chuẩn nhiều chiều

Số chiều: 3

  • \(H_0\) (Giả thuyết không): Tổng thể có phân phối chuẩn nhiều chiều.

  • \(H_1\) (Giả thuyết đối): Tổng thể không có phân phối chuẩn nhiều chiều.

Nhận xét Thực hiện kiểm định mshapiro thu được kết quả p-value=0.08396443>0.05, ta chấp nhận \(H_0\) với độ tin cậy 95%. Điều này nghĩa là tổng thể tuân theo phân phối chuẩn nhiều chiều.

Giả thuyết kiểm định tổng thể có phân phối chuẩn nhiều chiều đối với từng tổng thể giới tính (rùa đục)

Số chiều: 3

  • \(H_0\) (Giả thuyết không): Tổng thể của rùa đực có phân phối chuẩn nhiều chiều.

  • \(H_1\) (Giả thuyết đối): Tổng thể của rùa đực không có phân phối chuẩn nhiều chiều.

Nhận xét Thực hiện kiểm định mshapiro thu được kết quả p-value=0.05091876>0.05, ta chấp nhận \(H_0\) với độ tin cậy 95%. Điều này nghĩa là tổng thể của rùa đực tuân theo phân phối chuẩn nhiều chiều.

Giả thuyết kiểm định tổng thể có phân phối chuẩn nhiều chiều đối với từng tổng thể giới tính (rùa cái)

Số chiều: 3

  • \(H_0\) (Giả thuyết không): Tổng thể của rùa cái có phân phối chuẩn nhiều chiều.

  • \(H_1\) (Giả thuyết đối): Tổng thể của rùa cái không có phân phối chuẩn nhiều chiều.

Nhận xét Thực hiện kiểm định mshapiro thu được kết quả p-value=0.2077751>0.05, ta chấp nhận \(H_0\) với độ tin cậy 95%. Điều này nghĩa là tổng thể của rùa cái tuân theo phân phối chuẩn nhiều chiều.

Kiểm tra mô hình có tính đồng nhất của ma trận hiệp phương sai

box_m(data_45[, c("V1", "V2", "V3")], data_45$V4)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      9.95   0.127         6 Box's M-test for Homogeneity of Covariance Matric…

Giả thuyết kiểm định tính đồng nhất của ma trận hiệp phương sai

  • \(H_0\) (Giả thuyết không): Các tổng thể có chung ma trận hiệp phương sai.

  • \(H_1\) (Giả thuyết đối): Các tổng thể không có chung ma trận hiệp phương sai.

Nhận xét Thực hiện kiểm định Box’s M-test thu được kết quả p-value=0.1266747>0.05, ta chấp nhận \(H_0\) với độ tin cậy 95%. Điều này nghĩa là các tổng thểcó chung ma trận hiệp phương sai.

Kiểm tra mô hình có tính độc lập

library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
layout(matrix(c(1, 2), nrow = 1, ncol = 2))

results <- data_45 %>%
  select(V1, V2, V3, V4) %>%
  group_by(V4) %>%
  doo(~ggpairs(.) + theme_bw(), result = "plots")

results$plots
## [[1]]

## 
## [[2]]

Nhận xét Có hiện tượng đa cộng tuyến giữa các biến V1, V2, V3 vì hệ số tương quan đều lớn hơn 0.9.

câu (a) Kiểm định sự bằng nhau của hai vectơ trung bình tổng thể, với \(\alpha = 0.01\).

Với 24 con rùa cái và 24 con rùa đực, ta thực hiện giả thuyết kiểm định

  • \(H_0\) (Giả thuyết không):

    \[ \mu_{Female} = \mu_{Male} \]

vectơ trung bình của các biến Length, Width, Height không khác nhau giữa hai nhóm.

  • \(H_1\) (Giả thuyết đối):

    \[ \mu_{Female} \ne \mu_{Male} \]

có ít nhất một thành phần trong vectơ trung bình khác nhau giữa hai nhóm.

library(car)

mod <- Manova(lm(cbind(V1, V2, V3) ~ V4, data = data_45))
summary(mod)
## 
## Type II MANOVA Tests:
## 
## Sum of squares and products for error:
##            V1          V2         V3
## V1 966.821055 2.989880028 3.99919536
## V2   2.989880 0.009778764 0.01248563
## V3   3.999195 0.012485627 0.01773154
## 
## ------------------------------------------
##  
## Term: V4 
## 
## Sum of squares and products for the hypothesis:
##            V1          V2         V3
## V1 433.166629 1.435843804 2.67721266
## V2   1.435844 0.004759479 0.00887432
## V3   2.677213 0.008874320 0.01654668
## 
## Multivariate Tests: V4
##                  Df test stat approx F num Df den Df     Pr(>F)    
## Pillai            1 0.6633494 28.89977      3     44 1.7559e-10 ***
## Wilks             1 0.3366506 28.89977      3     44 1.7559e-10 ***
## Hotelling-Lawley  1 1.9704389 28.89977      3     44 1.7559e-10 ***
## Roy               1 1.9704389 28.89977      3     44 1.7559e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nhận xét Sử dụng hàm Manova(), ta quan sát thấy các giá trị Pr(>F) đều nhỏ hơn 0.01 nên ta bác bỏ \(H_0\) với mức ý nghĩa 0.01. Do đó, các vectơ trung bình về các chỉ số Length, Width, Height là khác nhau giữa hai giới tính va điều này có nghĩa là giới tính của loài rùa của ảnh hưởng đến các chỉ số cơ thể.

câu (b) Nếu giải thuyết trong câu (a) bị bác bỏ, tìm tổ hợp tuyến tính của các trung bình thành phần có ý nghĩa quan trọng trong việc bác bỏ \(H_0\).

library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:rstatix':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
lda_mod <- lda(V4 ~ V1 + V2 + V3, data = data_45)
cat("Tổ hợp tuyến tính từ LDA: \n")
## Tổ hợp tuyến tính từ LDA:
print(coef(lda_mod))
##             LD1
## V1    0.6209391
## V2   -3.6939251
## V3 -172.4978334

Nhận xét

  • V1 ứng với Length, V2 ứng với Width, V3 ứng với Height.

  • Height đóng góp lớn nhất do hệ số của Height khi lấy giá trị tuyệt đối thu được 172.4978334, điều này cho thấy Height là biến quan trọng nhất trong việc phân việc rùa cái và rùa đực.

  • Hướng ngược hệ số âm ở Height cho thấy sự khác biệt ngược chiều ở Height.

câu (c) Tìm các khyoanrg tin cậy đồng thời cho hiệu các trung bình thành phần. So sánh với các khoảng Bonferroni.

Các tham số

n1 <- nrow(female) # Số rùa cái
n2 <- nrow(male) # Số rùa đực
p <- 3 # Số biến
alpha <- 0.01

Các giá trị trung bình và ma trận hiệp phương sai

mean_f <- colMeans(female[, 1:3])
mean_m <- colMeans(male[, 1:3])
Sw_f <- cov(female[, 1:3])
Sw_m <- cov(male[, 1:3])

Ma trận hiệp phương sai gộp

Sw <- ((n1 - 1) * Sw_f + (n2 - 1) * Sw_m) / (n1 + n2 - 2)
Sw
##             V1           V2           V3
## V1 21.01784901 0.0649973919 0.0869390295
## V2  0.06499739 0.0002125818 0.0002714267
## V3  0.08693903 0.0002714267 0.0003854683

Hiệu trung bình

mean_diff <- mean_f - mean_m

Khoảng tin cậy đồng thời

F_critical <- qf(1 - alpha, p, n1 + n2 - p - 1)
T2_critical <- ((n1 + n2 -2) * p / (n1 + n2 - p - 1)) * F_critical
simul_critical <- sqrt(T2_critical * (1/n1 + 1/n2) * diag(Sw))

cat("Khoảng tin cậy đồng thời (alpha =", alpha, "): \n")
## Khoảng tin cậy đồng thời (alpha = 0.01 ):
for (i in 1:p) {
  lower <- mean_diff[i] - simul_critical[i]
  upper <- mean_diff[i] + simul_critical[i]
  cat(names(mean_diff)[i], ": [", lower, ", ", upper, "]\n")
}
## V1 : [ 1.170225 ,  10.84597 ]
## V2 : [ 0.004529491 ,  0.03530131 ]
## V3 : [ 0.01641509 ,  0.05785172 ]

Khoảng tin cậy Bonferroni

t_critical <- qt(1 - alpha / (2 * p), df = n1 + n2 - 2)
bonf_critical <- t_critical * sqrt((1/n1 + 1/n2) * diag(Sw))

cat("Khoảng tin cậy Bonferroni (alpha =", alpha, "): \n")
## Khoảng tin cậy Bonferroni (alpha = 0.01 ):
for (i in 1:p) {
  lower <- mean_diff[i] - bonf_critical[i]
  upper <- mean_diff[i] + bonf_critical[i]
  cat(names(mean_diff)[i], ": [", lower, ", ", upper, "]\n")
}
## V1 : [ 1.910385 ,  10.10581 ]
## V2 : [ 0.006883427 ,  0.03294738 ]
## V3 : [ 0.01958484 ,  0.05468196 ]

So sánh

  • V1 ứng với Length, V2 ứng với Width, V3 ứng với Height.

  • Các khoảng tin cậy Bonferroni đều nằm gọn trong khoảng tin cậy đồng thời.

  • Khoảng tin cậy đồng thời rộng hơn vì điều chỉnh đồng thời cho toàn bộ vectơ.

  • Cả hai loại khoảng tin cậy đều không chưa 0, khẳng định sự khác biệt trung bình giữa hai nhóm với mức ý nghĩa thống kê tại \(\alpha=0.01\).

  • Biến Length có độ rộng khoảng tin cậy đồng thời lớn nhất thể hiện độ biến thiên cao nhất, còn biến Height có khoảng tin cậy đồng thời hẹp nhất thể hiện độ ổn định tương đối.

Bài 4.6

Trong một nghiên cứu về thời gian nhận dạng số theo hai biểu diễn khác nhau: biểu diễn bằng từ và biểu diễn bằng ký số Ả Rập, người ta đưa ra hai con số có cùng tính chẵn/lẻ hoặc không cùng tính chẵn/lẻ và đo lường các biến X chỉ thời gian, tính bằng mili-giây, người được khảo sát trả lời đúng với

  • X1 dành cho trường hợp số không có cùng tính chẵn lẻ, diễn tả bằng từ;

  • X2 dành cho trường hợp số có cùng tính chẵn lẻ, diễn tả bằng từ;

  • X3 dành cho trường hợp số không có cùng tính chẵn lẻ, diễn tả bằng ký số Ả Rập;

  • X4 dành cho trường hợp số có cùng tính chẵn lẻ, diễn tả bằng ký số Ả Rập.

câu (a) Kiểm định sự bằng nhau giữa các phương thức, với \(\alpha=0.05\).

câu (b) Xây dựng các khoảng tin cậy đồng thời cho sự khác biệt giữa hai biểu diễn số, giữa có và không có cùng tính chẵn lẻ, và hiệu ứng tương tác giữa hai biểu diễn số và tính chẵn lẻ.

Mô tả dữ liệu bài 4.6

Dữ liệu trong bài toán 4.6 lnghiên cứu về thời gian nhận dạng số qua hai cách biểu diễn:biểu diễn bằng từ và biểu diễn bằng ký số Ả Rập, với 4 nhóm điều kiện trình bày khác nhau.

Dữ liệu bao gồm:

  • Biến số: X1, X2, X3, X4 chỉ thời gian phản hồi trung bình (đơn vị: ms) tương ứng với 4 biểu diễn khác nhau.

  • Điều kiện:

X1: không có cùng tính chẵn lẻ, diễn tả bằng từ.

X2: có cùng tính chẵn lẻ, diễn tả bằng từ.

X3: không có cùng tính chẵn lẻ, diễn tả bằng ký số Ả Rập.

X4: có cùng tính chẵn lẻ, diễn tả bằng ký số Ả Rập.

Tóm tắt đặc trưng dữ liệu:

  • Số biến phụ thuộc: 1
  • Số chiều của dữ liệu: 1 chiều
  • Loại phân tích: ANOVA 2_factor

Trước tiên ta tiến hành nạp dữ liệu.

# Đọc dữ liệu
data <- read.table("exo4-6.dat", header = FALSE)

colnames(data) = c("X1", "X2", "X3", "X4")

head(data)
##     X1    X2     X3  X4
## 1  869 860.5  691.0 601
## 2  995 875.0  678.0 659
## 3 1056 930.5  833.0 826
## 4 1126 954.0  888.0 728
## 5 1044 909.0  865.0 839
## 6  925 856.5 1059.5 797
data$SubjectID = factor(1:nrow(data))

Dài hóa và tách điều kiện

library(tidyr)
library(ggpubr)
library(dplyr)
long_data = pivot_longer(data,
                         cols = c("X1", "X2", "X3", "X4"),
                         names_to = "Condition",
                         values_to = "Response") %>%
  mutate(
    SubjectID = factor(rep(data$SubjectID, each = 4)),
    Parity = case_when(
      Condition %in% c("X1", "X3") ~ "Different",
      Condition %in% c("X2", "X4") ~ "Same"   
    ),
    Notation = case_when(
      Condition %in% c("X1", "X2") ~ "Word",
      Condition %in% c("X3", "X4") ~ "Digit"
    )
  )

Chuyển về dạng factor

long_data$Parity <- factor(long_data$Parity)
long_data$Notation <- factor(long_data$Notation)

head(long_data)
## # A tibble: 6 × 5
##   SubjectID Condition Response Parity    Notation
##   <fct>     <chr>        <dbl> <fct>     <fct>   
## 1 1         X1            869  Different Word    
## 2 1         X2            860. Same      Word    
## 3 1         X3            691  Different Digit   
## 4 1         X4            601  Same      Digit   
## 5 2         X1            995  Different Word    
## 6 2         X2            875  Same      Word

Kiểm tra mô hình có tuân theo phân phối chuẩn đơn biến

grouped = with(long_data, split(Response, interaction(Parity, Notation)))
normality_test = lapply(grouped, shapiro.test)

sapply(normality_test, function(x) x$p.value)
## Different.Digit      Same.Digit  Different.Word       Same.Word 
##      0.28534564      0.09384732      0.38360796      0.04634694

Giả thuyết kiểm định tổng thể có phân phối chuẩn đơn biến

  • \(H_0\) (Giả thuyết không): Biến Response trong từng điều kiện tuân theo phân phối chuẩn.

  • \(H_1\) (Giả thuyết đối): Biến Response trong từng điều kiện không tuân theo phân phối chuẩn.

Nhận xét Thực hiện kiểm định Shapiro-Wilk thu được các kết quả

  • Different.Digit p-value=0.28534564 > 0.05, chưa đủ bằng chứng để bác bỏ \(H_0\) nên dữ liệu tuân theo phân phối chuẩn .

  • Same.Digit p-value=0.09384732 > 0.05, chưa đủ bằng chứng để bác bỏ \(H_0\) nên dữ liệu tuân theo phân phối chuẩn.

  • Different.Word p-value=0.38360796 > 0.05, chưa đủ bằng chứng để bác bỏ \(H_0\) nên dữ liệu tuân theo phân phối chuẩn.

  • Same.Word p-value=0.04634694 < 0.05, có bằng chứng để bác bỏ \(H_0\) nên dữ liệu không tuân theo phân phối chuẩn.

Trong 4 điều kiện, có 3 nhóm có p-value > 0.05, chỉ riêng nhóm Same.Word với p-value=0.047 gần với 0.05 vi phạm giả định phân phối chuẩn. Nhưng đây là ANOVA đo lặp nên vẫn có thể tiếp tục với những vi phạm nhỏ về chuẩn.

Lưu ý: Trong trường hợp này bộ dữ liệu lặp được bỏ qua tính đồng nhất phương sai vì trong mỗi factor chỉ có 2 mức.

par(mfrow = c(2,2))
for (g in names(grouped)) {
  qqnorm(grouped[[g]], main = paste("Q-Q:", g))
  qqline(grouped[[g]])
}

library(rstatix)
long_data %>%
  group_by(Parity, Notation) %>%
  identify_outliers(Response)
## # A tibble: 1 × 7
##   Parity Notation SubjectID Condition Response is.outlier is.extreme
##   <fct>  <fct>    <fct>     <chr>        <dbl> <lgl>      <lgl>     
## 1 Same   Word     8         X2            1311 TRUE       FALSE

câu (a) Kiểm định sự bằng nhau giữa các phương thức, với \(\alpha=0.05\).

model = aov(Response ~ Parity * Notation + Error(SubjectID / (Parity*Notation)), data = long_data )

summary(model)
## 
## Error: SubjectID
##           Df  Sum Sq Mean Sq F value Pr(>F)
## Residuals 31 2256856   72802               
## 
## Error: SubjectID:Parity
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Parity     1 340570  340570   69.58 2.01e-09 ***
## Residuals 31 151725    4894                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: SubjectID:Notation
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## Notation   1 753608  753608   74.86 9.03e-10 ***
## Residuals 31 312087   10067                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: SubjectID:Parity:Notation
##                 Df Sum Sq Mean Sq F value Pr(>F)
## Parity:Notation  1   4022    4022   1.608  0.214
## Residuals       31  77557    2502

Giả thuyết kiểm định tổng quát:

  • \(H_0\) (Giả thuyết không): Không có sự khác biệt trung bình giữa bốn phương thức đo lường.

  • \(H_1\) (Giả thuyết đối): Có ít nhât một phương thức tạo ra kết quả khác biệt.

Để kiểm định điều này, ta sử dụng ANOVA đo lặp với 2 nhân tố

-Parity: Different - Same

-Notation: Word - Digit

Giả thuyết kiểm định theo từng yếu tố (main effect - Parity):

  • \(H_0\) (Giả thuyết không): Không có sự khác biệt giữa các mức Parity.

  • \(H_1\) (Giả thuyết đối): Có sự khác biệt giữa các mức Parity.

Nhận xét sử dụng hàm aov(), ta quan sát thấy Pr(>F) = 2.01e-09 < 0.05 nên ta bác bỏ \(H_0\) với mức ý nghĩa 0.05. Tức là có sự khác biệt đáng kể giữa các mức Parity.

Giả thuyết kiểm định theo từng yếu tố (main effect - Notation):

  • \(H_0\) (Giả thuyết không): Không có sự khác biệt giữa các mức Notation

  • \(H_1\) (Giả thuyết đối): Có sự khác biệt giữa các mức Notation

Nhận xét sử dụng hàm aov(), ta quan sát thấy Pr(>F) = 9.03e-10 < 0.05 nên ta bác bỏ \(H_0\) với mức ý nghĩa 0.05. Tức là có sự khác biệt đáng kể giữa các mức Notation.

Giả thuyết kiểm định theo từng yếu tố (interation):

  • \(H_0\) (Giả thuyết không): Không có sự tương tác giữa Parity và Notation.

  • \(H_1\) (Giả thuyết đối): Có sự tương tác giữa Parity và Notation.

Nhận xét sử dụng hàm aov(), ta quan sát thấy Pr(>F) = 0.214 > 0.05 chưa đủ bằng chứng để bác bỏ \(H_0\) với mức ý nghĩa 0.05. Tức là không có sự tương tác giữa Parity và Notation.

Kết luận Có sự khác biệt đáng kể giữa các phương thức đo lường, tuy nhiên không có sự tương tác giữa hai yếu tố, nghĩa là hiệu ứng tương tác của một yếu tố không phụ thuộc vào mức độ của yếu tố còn lại.

câu (b) Xây dựng các khoảng tin cậy đồng thời cho sự khác biệt giữa hai biểu diễn số, giữa có và không có cùng tính chẵn lẻ, và hiệu ứng tương tác giữa hai biểu diễn số và tính chẵn lẻ.

Xây dựng khoảng tin cậy đồng thời cho sự khác biệt giữa hai biểu diễn số

score_word = long_data$Response[long_data$Notation == "Word"]
score_digit = long_data$Response[long_data$Notation == "Digit"]

group_w = as.matrix(score_word)
group_d = as.matrix(score_digit)

n_w = nrow(group_w)
n_d = nrow(group_d)
p = ncol(group_w)
df = n_w + n_d - 2

mean_diff_notation = colMeans(group_d) - colMeans(group_w)
Sp_notation = ((n_w - 1) * cov(group_w) + (n_d -1) * cov(group_d)) / df
stderr_notation = sqrt(diag(Sp_notation) * (1/n_w + 1/n_d))

f_crit = qf(1 - alpha, p, df - p + 1)
factor = sqrt(p * (df - 1) / (df - p + 1) * f_crit)
t_bonf = qt( 1 - alpha / (2 * p), df)

result_notation = data.frame(
    Lower_Hotelling = mean_diff_notation - factor * stderr_notation,
    Upper_Hotelling = mean_diff_notation + factor * stderr_notation,
    Lower_Bonferroni = mean_diff_notation - t_bonf * stderr_notation,
    Upper_Bonferroni = mean_diff_notation + t_bonf * stderr_notation,
    row.names = "Notation: Digit - Word"
)
print(result_notation)
##                        Lower_Hotelling Upper_Hotelling Lower_Bonferroni
## Notation: Digit - Word       -226.1902       -80.73168        -226.4805
##                        Upper_Bonferroni
## Notation: Digit - Word        -80.44134

Xây dựng khoảng tin cậy đồng thời giữa có và không có cùng tính chẵn lẻ

score_same = long_data$Response[long_data$Parity == "Same"]
score_diff = long_data$Response[long_data$Parity == "Different"]

group_same = as.matrix(score_same)
group_diff = as.matrix(score_diff)

n_same = nrow(group_same)
n_diff = nrow(group_diff)
df = n_same + n_diff - 2

mean_diff_parity = colMeans(group_diff) - colMeans(group_same)
Sp_parity = ((n_same - 1) * cov(group_same) + (n_diff - 1) * cov(group_diff)) / df
stderr_parity = sqrt(diag(Sp_parity) * (1 / n_same + 1 / n_diff))

result_parity = data.frame(
    Lower_Hotelling = mean_diff_parity - factor * stderr_parity,
    Upper_Hotelling = mean_diff_parity + factor * stderr_parity,
    Lower_Bonferroni = mean_diff_parity - t_bonf * stderr_parity,
    Upper_Bonferroni = mean_diff_parity + t_bonf * stderr_parity,
    row.names = "Parity: Different - Same"
)
print(result_parity)
##                          Lower_Hotelling Upper_Hotelling Lower_Bonferroni
## Parity: Different - Same        25.80314         180.525         25.49431
##                          Upper_Bonferroni
## Parity: Different - Same         180.8338

Khoảng tin cậy đồng thời hiệu ứng tương tác giữa hai biểu diễn số và tính chẵn lẻ

library(dplyr)

subject_means = long_data %>%
  group_by(SubjectID, Parity, Notation) %>%
  summarize(Mean = mean(Response), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = c(Parity, Notation), values_from = Mean)

interaction_diff = with(subject_means,
    (`Different_Digit` - `Different_Word`) - (`Same_Digit` - `Same_Word`))

group_interaction = as.matrix(interaction_diff)
n_inter = length(interaction_diff)
p = 1
df = n_inter - 1

mean_diff_inter = colMeans(group_interaction)
Sp_inter = var(group_interaction)
stderr_inter = sqrt(Sp_inter * (1 / n_inter))

f_crit = qf(1 - alpha, p, df - p + 1)
factor = sqrt(p * (df - 1) / (df - p + 1) * f_crit)
t_bonf = qt(1 - alpha / 2, df)

result_interaction = data.frame(
  Lower_Hotelling = mean_diff_inter - factor * stderr_inter,
  Upper_Hotelling = mean_diff_inter + factor * stderr_inter,
  Lower_Bonferroni = mean_diff_inter - t_bonf * stderr_inter,
  Upper_Bonferroni = mean_diff_inter + t_bonf * stderr_inter,
  row.names = "Interaction: (Digit - Word) - (Different - Same)"
)

print(result_interaction)
##                                                  Lower_Hotelling
## Interaction: (Digit - Word) - (Different - Same)        -25.3151
##                                                  Upper_Hotelling
## Interaction: (Digit - Word) - (Different - Same)        70.15885
##                                                  Lower_Bonferroni
## Interaction: (Digit - Word) - (Different - Same)        -26.10419
##                                                  Upper_Bonferroni
## Interaction: (Digit - Word) - (Different - Same)         70.94794

Bài 4.7

câu (a) Dùng Hotelling’s \(T^2\) để kiểm tra sự khác biệt trung bình giữa hai nhóm (DB vs WO) với 4 biến: Height, Weight, ArmLen, HandLen

Câu (b) Xây dựng khoảng tin cậy đồng thời 95% cho 4 biến: Height, Weight, ArmLen, HandLen

câu (c) Kiểm tra các điều kiện cho MANOVA có thỏa với dữ liệu này không ? Giải thích ?

Câu (d) Tìm tổ hợp tuyến tính tối ưu của 4 biến: Height, Weight, ArmLen, HandLen sao cho sai biệt trung bình giữa hai nhóm “DB” và “WO” là lớn nhất

Mô tả dữ liệu bài 4.7

Dữ liệu trong bài toán 4.7 liên quan đến các đặc trưng thể chất của các cầu thủ trong giải bóng bầu dục Mỹ (NFL). Mỗi cầu thủ được phân vào một trong hai vị trí thi đấu:

  • DB (Defensive Back): Hậu vệ phòng thủ
  • WO (Wideout): Cầu thủ tấn công ở biên

Mỗi cầu thủ được mô tả bởi 4 biến liên tục đại diện cho các chỉ số hình thể:

  • Height: Chiều cao (inch)
  • Weight: Cân nặng (pound)
  • ArmLen: Chiều dài tay (inch)
  • HandLen: Chiều dài bàn tay (inch)

Tóm tắt đặc trưng dữ liệu:

  • Số biến phụ thuộc: 4
  • Số nhóm (nhân tố): 2 vị trí (DB, WO)
  • Số chiều của dữ liệu: 4 chiều (mỗi cầu thủ có 4 chỉ số thể hình)
  • Loại phân tích: One-way MANOVA

câu (a) Dùng Hotelling’s \(T^2\) để kiểm tra sự khác biệt trung bình giữa hai nhóm (DB vs WO) với 4 biến: Height, Weight, ArmLen, HandLen

Giả thuyết kiểm định:

  • \(H_0\) (Giả thuyết không):

    \[ \mu_{DB} = \mu_{WO} \]

vector trung bình của các biến Height, Weight, ArmLen, HandLen không khác nhau giữa 2 nhóm.

  • \(H_1\) (Giả thuyết đối):

    \[ \mu_{DB} \ne \mu_{WO} \]

có ít nhất một thành phần trong vector trung bình khác nhau giữa hai nhóm.

# Đọc dữ liệu
nfl <- read.csv("NFL 2024.csv")

Đọc dữ liệu từ file CSV vào biến NFL

# Lọc dữ liệu theo hai nhóm Position = DB và WO
nfl_sub <- nfl[nfl$Position %in% c("DB", "WO"), ]
group <- factor(nfl_sub$Position)

Lọc lại dữ liệu, chỉ giữ lại các cầu thủ có Position là DB hoặc WO. Tạo biến định danh nhóm (group factor) cho phân tích thống kê (có 2 mức: DB và WO).

# Chọn các biến để phân tích
Y <- nfl_sub[, c("Height", "Weight", "ArmLen", "HandLen")]

Y <- nfl_sub[, c(“Height”, “Weight”, “ArmLen”, “HandLen”)].

# Tách nhóm
Y1 <- Y[group == "DB", ]
Y2 <- Y[group == "WO", ]

Tách ma trận Y thành hai nhóm theo Position.

# Kích thước mẫu và số biến
n1 <- nrow(Y1)
n2 <- nrow(Y2)
p <- ncol(Y)

Lưu lại: n1: số mẫu nhóm DB n2: số mẫu nhóm WO p: số biến phân tích (4 biến).

# Tính hiệu trung bình và ma trận hiệp phương sai gộp
dbar <- colMeans(Y1) - colMeans(Y2)
Spooled <- ((n1 - 1)*cov(Y1) + (n2 - 1)*cov(Y2)) / (n1 + n2 - 2)

Tính hiệu trung bình giữa hai nhóm (vector \(\bar{Y}_1 - \bar{Y}_2\)).

Tính ma trận hiệp phương sai gộp \(S_p\) dùng chung cho cả hai nhóm.

# Tính thống kê Hotelling's T² và chuyển sang F
T2 <- (n1 * n2) / (n1 + n2) * t(dbar) %*% solve(Spooled) %*% dbar
F0 <- as.numeric((n1 + n2 - p - 1)*T2 / ((n1 + n2 - 2)*p))

Tính Hotelling’s \(T^2\)** theo công thức: \[ T^2 = \frac{n_1 n_2}{n_1 + n_2} (\bar{Y}_1 - \bar{Y}_2)^T S_p^{-1} (\bar{Y}_1 - \bar{Y}_2) \]

Chuyển \(T^2\) sang thống kê \(F\) để so sánh với bảng phân phối \(F\):

\[ F_0 = \left( \frac{n_1 + n_2 - p - 1}{(n_1 + n_2 - 2)p} \right) \cdot T^2 \]

# Ngưỡng tới hạn và p-value
alpha <- 0.05
F_crit <- qf(1 - alpha, df1 = p, df2 = n1 + n2 - p - 1)
p_val <- 1 - pf(F0, df1 = p, df2 = n1 + n2 - p - 1)

Tính ngưỡng tới hạn \(F_{\alpha; p,\; n_1 + n_2 - p - 1}\) với mức ý nghĩa \(\alpha = 0{,}05\)

Tính p-value của kiểm định để so sánh với \(\alpha\)

# In kết quả
cat("T² =", T2, "\n")
## T² = 14.47333
cat("F =", F0, "\n")
## F = 3.490627
cat("Ngưỡng tới hạn F_crit =", F_crit, "\n")
## Ngưỡng tới hạn F_crit = 2.483034
cat("p-value =", p_val, "\n")
## p-value = 0.0110447

Nhận xét

Vì giá trị thống kê kiểm định là \(F = 3.02 > F_{\text{crit}} = 2.483034\) và giá trị p-value là \(p = 0.0110447 < \alpha = 0.05\), ta bác bỏ giả thuyết không \(H_0\) ở mức ý nghĩa 5%.

Điều này cho thấy có bằng chứng thống kê đủ mạnh để kết luận rằng vector trung bình của các biến định lượng Height, Weight, ArmLen, HandLen không giống nhau giữa hai nhóm DB và WO.

Dữ liệu cho thấy có sự khác biệt có ý nghĩa thống kê giữa hai nhóm về ít nhất một trong các biến đã phân tích.

Câu (b) Xây dựng khoảng tin cậy đồng thời 95% cho 4 biến: Height, Weight, ArmLen, HandLen

# Lọc dữ liệu cho 4 biết
Y_sub <- Y[, c("Height", "Weight", "ArmLen", "HandLen")]
Y1_sub <- Y1[, c("Height", "Weight", "ArmLen", "HandLen")]
Y2_sub <- Y2[, c("Height", "Weight", "ArmLen", "HandLen")]

Chọn đúng 4 biến cần phân tích. Dữ liệu được chia thành:

Y1_sub: nhóm DB

Y2_sub: nhóm WO

# Kích thước mới
p_sub <- ncol(Y_sub)

# Hiệu trung bình
dbar_sub <- colMeans(Y1_sub) - colMeans(Y2_sub)

Tính hiệu trung bình mẫu cho từng biến:

\[ \bar{Y}_{DB} - \bar{Y}_{WO} \]

Đây là ước lượng điểm (estimate) cho hiệu trung bình cần xây dựng khoảng tin cậy.

# Ma trận hiệp phương sai gộp cho 4 biến
Spooled_sub <- ((n1 - 1)*cov(Y1_sub) + (n2 - 1)*cov(Y2_sub)) / (n1 + n2 - 2)
sp_diag <- diag(Spooled_sub)

Tính ma trận hiệp phương sai gộp \(S_p\) từ hai nhóm.

sp_diag lấy phần đường chéo của \(S_p\), chính là phương sai của từng biến.

# Hệ số hiệu chỉnh từ phân phối F
alpha <- 0.05
c_alpha <- (p_sub * (n1 + n2 - 2)) / (n1 + n2 - p_sub - 1) * qf(1 - alpha, p_sub, n1 + n2 - p_sub - 1)

Tính hệ số hiệu chỉnh \(c_\alpha\) dùng để mở rộng khoảng tin cậy, nhằm đảm bảo xác suất đúng đồng thời \(\geq 95\%\).

# Tính khoảng tin cậy đồng thời
CI_matrix <- matrix(NA, nrow = p_sub, ncol = 3)
colnames(CI_matrix) <- c("Estimate", "Lower CI", "Upper CI")
rownames(CI_matrix) <- colnames(Y_sub)

for (j in 1:p_sub) {
  se <- sqrt(c_alpha * (1/n1 + 1/n2) * sp_diag[j])
  CI_matrix[j, 1] <- dbar_sub[j]
  CI_matrix[j, 2] <- dbar_sub[j] - se
  CI_matrix[j, 3] <- dbar_sub[j] + se
}

Lặp qua từng biến để tính khoảng tin cậy:

\[ \text{Estimate} \pm \text{Sai số chuẩn hiệu chỉnh} = \bar{y}_{1j} - \bar{y}_{2j} \pm \left( c_\alpha \cdot \sqrt{ \left( \frac{1}{n_1} + \frac{1}{n_2} \right) \cdot s_{jj} } \right) \]

Lưu vào bảng CI_matrix gồm: Estimate, Lower CI, Upper CI.

# In kết quả
print(round(CI_matrix, 4))
##         Estimate Lower CI Upper CI
## Height   -1.3654  -2.8651   0.1343
## Weight   -6.1619 -15.6465   3.3228
## ArmLen   -0.5913  -1.3730   0.1903
## HandLen  -0.2676  -0.6030   0.0678

Câu (c) – Xây dựng khoảng tin cậy đồng thời 95% Bonferroni cho 4 biến:Height, Weight, ArmLen, HandLen

# Số biến
p <- length(dbar)
alpha <- 0.05
alpha_star <- alpha / (2 * p)  # Bonferroni adjustment
# Giá trị tới hạn t
t_crit <- qt(1 - alpha_star, df = n1 + n2 - 2)
# Ma trận phương sai (đường chéo)
sp_diag <- diag(Spooled)
# Tính khoảng tin cậy
CI_bonf <- matrix(NA, nrow = p, ncol = 3)
colnames(CI_bonf) <- c("Estimate", "Lower CI", "Upper CI")
rownames(CI_bonf) <- colnames(Y1)

for (j in 1:p) {
  se <- sqrt((1/n1 + 1/n2) * sp_diag[j])
  margin <- t_crit * se
  CI_bonf[j, 1] <- dbar[j]
  CI_bonf[j, 2] <- dbar[j] - margin
  CI_bonf[j, 3] <- dbar[j] + margin
}

  • Với mỗi biến \(j\):

    • Tính sai số chuẩn:

      \[ SE_j = \sqrt{ \left( \frac{1}{n_1} + \frac{1}{n_2} \right) \cdot s_{jj} } \]

    • Tính biên độ sai số hiệu chỉnh:
      \(\text{margin} = t_{\alpha^*} \cdot SE_j\)

    • Từ đó tạo khoảng tin cậy:
      \[ \hat{\delta}_j \pm \text{margin} \]

    print(round(CI_bonf, 4))
    ##         Estimate Lower CI Upper CI
## Height   -1.3654  -2.5582  -0.1726
## Weight   -6.1619 -13.7054   1.3817
## ArmLen   -0.5913  -1.2130   0.0304
## HandLen  -0.2676  -0.5344  -0.0009

Câu (d) Tìm tổ hợp tuyến tính tối ưu của 4 biến: Height, Weight, ArmLen, HandLen sao cho sai biệt trung bình giữa hai nhóm “DB” và “WO” là lớn nhất

Mục tiêu:

Tìm vector trọng số a sao cho tổ hợp:

\[ Z = \mathbf{a}^\top \mathbf{x} = a_1 \cdot \text{Height} + a_2 \cdot \text{Weight} + a_3 \cdot \text{ArmLen} + a_4 \cdot \text{HandLen} \]

có khả năng phân biệt mạnh nhất giữa hai nhóm.

# Tạo ma trận ngoài: d %*% t(d)
ddT <- dbar %*% t(dbar)

dbar là vector hiệu trung bình \(\mathbf{d} = \bar{y}_1 - \bar{y}_2\)

Tính tích ngoài: \(\mathbf{d}\mathbf{d}^\top\), gọi là ma trận phương sai giữa nhóm (B)

# Ma trận cần lấy trị riêng
M <- solve(Spooled) %*% ddT

Đây chính là ma trận cần phân tích trị riêng để giải bài toán tối ưu:

\[ M = S_p^{-1} \cdot \mathbf{d}\mathbf{d}^\top \]

# Tính trị riêng
eig <- eigen(M)
# Lấy trị riêng lớn nhất
lambda_max <- eig$values[1]
a_max <- eig$vectors[,1]  # Vector riêng ứng với lambda_max

Tính trị riêng và vector riêng của ma trận \(M\)

Lấy trị riêng lớn nhất (ứng với hướng phân biệt tốt nhất)

a_max: là tổ hợp tuyến tính tối ưu (trục phân biệt thứ nhất)

# Chuẩn hóa a để có tổng bình phương bằng 1
a_max <- a_max / sqrt(sum(a_max^2))
# Hiển thị
names(a_max) <- colnames(Y1)

Chuẩn hóa sao cho:

\[ \|\mathbf{a}\|^2 = \sum a_i^2 = 1 \]

Giúp đảm bảo vector a có độ dài 1, dễ so sánh và giải thích.

cat("Tổ hợp tuyến tính phân biệt hai nhóm DB và WO tốt nhất là:\n\n")
## Tổ hợp tuyến tính phân biệt hai nhóm DB và WO tốt nhất là:
# In hệ số tổ hợp tuyến tính đã chuẩn hóa
print(round(a_max, 4))
##  Height  Weight  ArmLen HandLen 
##  0.2931 -0.0119 -0.0434  0.9550
terms <- paste0(round(a_max, 4), " * ", names(a_max))
cat("\nZ = ", paste(terms, collapse = " + "), "\n\n")
## 
## Z =  0.2931 * Height + -0.0119 * Weight + -0.0434 * ArmLen + 0.955 * HandLen

Z = 0.2931 * Height + -0.0119 * Weight + -0.0434 * ArmLen + 0.955 * HandLen

Các hệ số đã được chuẩn hóa sao cho tổng bình phương bằng 1

Hệ số nào có giá trị tuyệt đối càng lớn thì biến tương ứng càng góp phần mạnh vào việc phân biệt hai nhóm.

Dấu của hệ số cho biết hướng phân biệt: nếu \(a_i > 0\), thì nhóm có giá trị biến lớn hơn sẽ có Z lớn hơn, và ngược lại.

Bài 4.8

Mô tả bài toán và dữ liệu

Bài toán này dựa trên nghiên cứu “Biomechanical Comparison of the FasT-Fix Meniscal Repair Suture System with Vertical Mattress Sutures and Meniscus Arrows” (Borden, J., Nyland, D.N.M.Caborn, D.Pienkowski (2003), The American Journal of Sports Medicine, Vol. 31, #3, pp. 374-378).

Mục tiêu là kiểm định xem ba phương pháp cố định sụn chêm khác nhau có tác động đến sự khác biệt về các đặc tính cơ học hay không. Các đặc tính này bao gồm:

  • Load at failure (N): Tải trọng tại thời điểm hỏng (ký hiệu Load).
  • Displacement (mm): Độ dịch chuyển (ký hiệu Displacement).
  • Stiffness (N/mm): Độ cứng (ký hiệu Stiffness).

Có 3 phương pháp được so sánh:

  1. Vertical Suture (Khâu dọc - Nhóm 1)
  2. Meniscus Arrow (Mũi tên sụn chêm - Nhóm 2)
  3. FasT-Fix (Nhóm 3)

Dữ liệu (meniscus.dat) gồm 6 mẫu thử nghiệm (đầu gối) cho mỗi phương pháp. Chúng ta sẽ sử dụng Phân tích phương sai đa biến một yếu tố (One-way MANOVA) để kiểm định sự khác biệt của vector trung bình của ba biến đáp ứng giữa ba phương pháp, với mức ý nghĩa \(\\alpha = 0.05\). Thống kê kiểm định được sử dụng là Wilks’ Lambda.

Tóm tắt đặc trưng dữ liệu:

  • Số biến phụ thuộc: 3 (Load, Displacement, Stiffness)

  • Số nhóm (nhân tố Method): 3 (Vertical Suture, Meniscus Arrow, FasT-Fix)

  • Số chiều của dữ liệu: 3 chiều (mỗi mẫu có 3 chỉ số)

  • Số mẫu mỗi nhóm: 6

  • Tổng số mẫu: 18

  • Loại phân tích: MANOVA 1_factor (One-way MANOVA)

    library(dplyr)
library(rstatix)
library(ggpubr)
library(car)
library(MVN) 
library(knitr)
# Tạo dataframe từ dữ liệu đã cho
data_text_4_8 <- "
1   97.3    16.9     8.3
1   106.4   20.2     7.2
1   118.2   20.1     6.3
1    99.7   15.7     7.3
1   106.5   13.9     8.7
1    84.2   14.9     8.7
2    44.9    7.9     4.7
2    46.1   12.5     6.1
2    59.3   15.5     5.0
2    35.5   10.2     5.8
2    50.7    8.9     6.6
2    56.8   13.3     8.4
3    88.0   18.0     8.0
3   119.8   18.5     8.3
3    65.8    9.2     7.6
3    82.9   18.8     6.4
3   149.9   22.8     8.2
3   117.1   17.5     7.7
"
# Đọc dữ liệu từ text connection
data_meniscus <- read.table(text = data_text_4_8, header = FALSE,
                           col.names = c("Method", "Load", "Displacement", "Stiffness"))

# Chuyển cột 'Method' thành factor
data_meniscus$Method <- factor(data_meniscus$Method,
                              labels = c("Vertical_Suture", "Meniscus_Arrow", "FasT_Fix"))

cat("Cấu trúc dữ liệu:\n")
## Cấu trúc dữ liệu:
str(data_meniscus)
## 'data.frame':    18 obs. of  4 variables:
##  $ Method      : Factor w/ 3 levels "Vertical_Suture",..: 1 1 1 1 1 1 2 2 2 2 ...
##  $ Load        : num  97.3 106.4 118.2 99.7 106.5 ...
##  $ Displacement: num  16.9 20.2 20.1 15.7 13.9 14.9 7.9 12.5 15.5 10.2 ...
##  $ Stiffness   : num  8.3 7.2 6.3 7.3 8.7 8.7 4.7 6.1 5 5.8 ...
cat("\n\n6 dòng đầu của dữ liệu:\n")
## 
## 
## 6 dòng đầu của dữ liệu:
head(data_meniscus)
##            Method  Load Displacement Stiffness
## 1 Vertical_Suture  97.3         16.9       8.3
## 2 Vertical_Suture 106.4         20.2       7.2
## 3 Vertical_Suture 118.2         20.1       6.3
## 4 Vertical_Suture  99.7         15.7       7.3
## 5 Vertical_Suture 106.5         13.9       8.7
## 6 Vertical_Suture  84.2         14.9       8.7
cat("\n\nSố lượng mẫu mỗi nhóm:\n")
## 
## 
## Số lượng mẫu mỗi nhóm:
print(table(data_meniscus$Method))
## 
## Vertical_Suture  Meniscus_Arrow        FasT_Fix 
##               6               6               6

Nhận xét: Dữ liệu đã được nạp thành công. Có 3 nhóm phương pháp, mỗi nhóm có 6 quan sát. Các biến Load, Displacement, Stiffness là số thực, biến Method đã được chuyển thành factor.

Kiểm tra các điều kiện cho MANOVA

MANOVA yêu cầu một số giả định:

  1. Các quan sát độc lập.
  2. Các vector quan sát trong mỗi nhóm tuân theo phân phối chuẩn nhiều chiều.
  3. Ma trận hiệp phương sai của các nhóm bằng nhau (đồng nhất).

Tính độc lập của các quan sát

Giả định này thường được đảm bảo bởi thiết kế nghiên cứu. Trong trường hợp này, 6 mẫu thử nghiệm (đầu gối) cho mỗi phương pháp có thể được coi là độc lập với nhau.

Phân phối chuẩn nhiều chiều

Với kích thước mẫu nhỏ (\(n\_k=6\) cho mỗi nhóm), việc kiểm định chính thức tính chuẩn nhiều chiều có thể không đáng tin cậy. Thay vào đó, chúng ta sẽ kiểm tra tính chuẩn đơn biến cho từng biến trong mỗi nhóm.

Kiểm tra chuẩn đơn biến (Shapiro-Wilk Test) và trực quan hóa (Q-Q plots, Histogram)

Giả thuyết kiểm định Shapiro-Wilk:

  • \(H\_0\): Dữ liệu của biến tuân theo phân phối chuẩn.
  • \(H\_1\): Dữ liệu của biến không tuân theo phân phối chuẩn. Chúng ta sẽ bác bỏ \(H\_0\) nếu p-value < 0.05.
normality_tests_4_8 <- data_meniscus %>%
  group_by(Method) %>%
  summarise(
    shapiro_Load = shapiro.test(Load)$p.value,
    shapiro_Displacement = shapiro.test(Displacement)$p.value,
    shapiro_Stiffness = shapiro.test(Stiffness)$p.value,
    .groups = 'drop'
  )
cat("Kết quả kiểm định Shapiro-Wilk (p-values):\n")
## Kết quả kiểm định Shapiro-Wilk (p-values):
kable(normality_tests_4_8, caption = "P-values từ kiểm định Shapiro-Wilk cho tính chuẩn đơn biến")
P-values từ kiểm định Shapiro-Wilk cho tính chuẩn đơn biến
Method shapiro_Load shapiro_Displacement shapiro_Stiffness
Vertical_Suture 0.8911563 0.3035952 0.3365944
Meniscus_Arrow 0.8349323 0.8373858 0.5495798
FasT_Fix 0.8000031 0.1268883 0.1209486

Nhận xét: Tất cả các p-value từ kiểm định Shapiro-Wilk cho từng biến trong từng nhóm đều lớn hơn 0.05. Do đó, không có bằng chứng nào cho thấy sự vi phạm giả định phân phối chuẩn đơn biến cho bất kỳ biến nào trong bất kỳ nhóm phương pháp nào ở mức ý nghĩa 5%.

# Trực quan hóa bằng Q-Q plots
par(mfrow = c(3, 3)) # Sắp xếp 3x3 plots

# Load
ggqqplot(data_meniscus, "Load", color = "Method", ggtheme = theme_bw()) + facet_wrap(~Method) + ggtitle("Q-Q Plot: Load by Method")

# Displacement
ggqqplot(data_meniscus, "Displacement", color = "Method", ggtheme = theme_bw()) + facet_wrap(~Method) + ggtitle("Q-Q Plot: Displacement by Method")

# Stiffness
ggqqplot(data_meniscus, "Stiffness", color = "Method", ggtheme = theme_bw()) + facet_wrap(~Method) + ggtitle("Q-Q Plot: Stiffness by Method")

Nhận xét: Hầu hết các điểm dữ liệu tương đối gần với đường thẳng và trong khoảng tin cậy -> tuân theo phân phối chuẩn.

# Trực quan hóa bằng Histograms
# Load
hist(data_meniscus$Load[data_meniscus$Method=="Vertical_Suture"], main="Load - Vertical Suture", xlab="Load")

hist(data_meniscus$Load[data_meniscus$Method=="Meniscus_Arrow"], main="Load - Meniscus Arrow", xlab="Load")

hist(data_meniscus$Load[data_meniscus$Method=="FasT_Fix"], main="Load - FasT-Fix", xlab="Load")

# Displacement
hist(data_meniscus$Displacement[data_meniscus$Method=="Vertical_Suture"], main="Displacement - Vertical Suture", xlab="Displacement")

hist(data_meniscus$Displacement[data_meniscus$Method=="Meniscus_Arrow"], main="Displacement - Meniscus Arrow", xlab="Displacement")

hist(data_meniscus$Displacement[data_meniscus$Method=="FasT_Fix"], main="Displacement - FasT-Fix", xlab="Displacement")

# Stiffness
hist(data_meniscus$Stiffness[data_meniscus$Method=="Vertical_Suture"], main="Stiffness - Vertical Suture", xlab="Stiffness")

hist(data_meniscus$Stiffness[data_meniscus$Method=="Meniscus_Arrow"], main="Stiffness - Meniscus Arrow", xlab="Stiffness")

hist(data_meniscus$Stiffness[data_meniscus$Method=="FasT_Fix"], main="Stiffness - FasT-Fix", xlab="Stiffness")

Kiểm tra chuẩn nhiều chiều (Multivariate Normality)

Các kiểm định chuẩn nhiều chiều (ví dụ: mshapiro.test từ gói MVN hoặc các kiểm định Mardia) thường yêu cầu cỡ mẫu lớn hơn nhiều so với số biến để có kết quả ổn định. Với \(N=18\) và \(p=3\), kiểm định này có thể không phù hợp. Chúng ta sẽ bỏ qua kiểm định này do cỡ mẫu nhỏ.

Tính đồng nhất của ma trận hiệp phương sai (Homogeneity of Covariance Matrices)

Sử dụng kiểm định Box’s M để kiểm tra giả định này.

Giả thuyết kiểm định Box’s M:

  • \(H\_0\): Các ma trận hiệp phương sai của các nhóm là bằng nhau.
  • \(H\_1\): Ít nhất hai ma trận hiệp phương sai của các nhóm là khác nhau. Chúng ta sẽ bác bỏ \(H\_0\) (tức là ma trận hiệp phương sai không đồng nhất) nếu p-value < \(\\alpha\_{BoxM}\) (thường chọn \(\\alpha\_{BoxM}\) nghiêm ngặt hơn, ví dụ 0.01 hoặc 0.001, do kiểm định này nhạy với vi phạm chuẩn).
# Kiểm tra bằng Box's M test sử dụng hàm từ rstatix
Y_4_8 <- data_meniscus[, c("Load", "Displacement", "Stiffness")]
# box_m_result_4_8 <- box_m(Y_4_8, data_meniscus$Method) # Cách 1: data frame và vector nhóm
box_m_result_4_8 <- box_m(data_meniscus[, c("Load", "Displacement", "Stiffness")], data_meniscus$Method)


cat("Kết quả kiểm định Box's M (sử dụng rstatix::box_m):\n")
## Kết quả kiểm định Box's M (sử dụng rstatix::box_m):
print(box_m_result_4_8)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      15.7   0.203        12 Box's M-test for Homogeneity of Covariance Matric…

Nhận xét (dựa trên kết quả giả định của Box’s M):

  • Kiểm định Box’s M cho ra giá trị thống kê chi-square và p-value.
  • p-value > 0.01 (mức \(\\alpha\) ), chúng ta không bác bỏ \(H\_0\) và có thể kết luận rằng giả định về tính đồng nhất của ma trận hiệp phương sai được thỏa mãn.

3. Thực hiện MANOVA

Chúng ta sẽ thực hiện MANOVA để kiểm tra xem có sự khác biệt ý nghĩa thống kê nào về vector trung bình của Load, Displacement, và Stiffness giữa ba phương pháp cố định hay không.

Giả thuyết kiểm định MANOVA:

Gọi \(\mathbf{\mu}_k\) là vector trung bình của \(p=3\) biến phụ thuộc (Load, Displacement, Stiffness) cho phương pháp thứ \(k\), với \(k \in \{1, 2, 3\}\).

  • Giả thuyết không (\(H_0\)): Các vector trung bình của ba biến phụ thuộc là như nhau qua cả ba phương pháp. \[ H_0: \mathbf{\mu}_1 = \mathbf{\mu}_2 = \mathbf{\mu}_3 \]

  • Giả thuyết đối (\(H_1\)): Tồn tại ít nhất một cặp phương pháp \((i, j)\) với \(i \neq j\) sao cho vector trung bình của các biến phụ thuộc của chúng khác nhau. \[ H_1: \text{Tồn tại ít nhất một } (i, j) \text{ sao cho } \mathbf{\mu}_i \neq \mathbf{\mu}_j \text{, với } i, j \in \{1, 2, 3\}, i \neq j \] (Nói cách khác, không phải tất cả các vector trung bình \(\mathbf{\mu}_k\) đều bằng nhau.)

Mức ý nghĩa được sử dụng cho kiểm định là \(\alpha = 0.05\).

# Thực hiện MANOVA
manova_model_4_8 <- manova(cbind(Load, Displacement, Stiffness) ~ Method, data = data_meniscus)

# Xem kết quả với thống kê Wilks' Lambda
summary_manova_4_8_wilks <- summary(manova_model_4_8, test = "Wilks")

cat("Kết quả MANOVA (sử dụng Wilks' Lambda):\n")
## Kết quả MANOVA (sử dụng Wilks' Lambda):
print(summary_manova_4_8_wilks)
##           Df  Wilks approx F num Df den Df   Pr(>F)   
## Method     2 0.2868   3.7582      6     26 0.007949 **
## Residuals 15                                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Phân tích kết quả MANOVA (dựa trên kết quả giả định): Bảng kết quả summary_manova_4_8_wilks sẽ cung cấp các thông tin sau cho yếu tố Method:

  • Df: Bậc tự do của giả thuyết (cho yếu tố Method). Thường là \(g-1 = 3-1 = 2\). (Lưu ý: R sẽ nhân với số biến p, nên Df trong output của Wilks cho tử số của F thường là \(p(g-1)\)).
  • Wilks: Giá trị thống kê Wilks’ Lambda. Giá trị càng nhỏ (gần 0) càng cho thấy sự khác biệt lớn giữa các nhóm.
  • approx F: Giá trị F xấp xỉ được chuyển đổi từ Wilks’ Lambda.
  • num Df: Bậc tự do tử số của F xấp xỉ. Sẽ là \(p(g-1) = 3 \\times (3-1) = 6\).
  • den Df: Bậc tự do mẫu số của F xấp xỉ.
  • Pr(>F): P-value của kiểm định.

Nhận xét:

  • p-value (Pr(>F)) < 0.05, chúng ta bác bỏ giả thuyết không \(H\_0\). Điều này có nghĩa là có bằng chứng thống kê cho thấy ít nhất một phương pháp cố định tạo ra một bộ đặc tính cơ học (Load, Displacement, Stiffness) trung bình khác biệt so với các phương pháp khác hay Có bằng chứng thống kê đủ để kết luận rằng ba phương pháp cố định sụn chêm (Vertical Suture, Meniscus Arrow, FasT-Fix) có tác động khác nhau đến tổ hợp các đặc tính cơ học trung bình (Tải trọng tại thời điểm hỏng, Độ dịch chuyển, và Độ cứng) với mức ý nghĩa 5%.

Bài 4.9

Mô tả dữ liệu

Dữ liệu bao gồm thời gian thực hiện nhiệm vụ trên thiết bị di động với:

  • Nhân tố A: Kích thước màn hình (5”, 7”)
  • Nhân tố B: Loại giao diện (1 = Easy, 2 = Difficult, 3 = Normal)

Biến phụ thuộc: Thời gian hoàn thành nhiệm vụ

Các điều kiện thí nghiệm:

  • Trt11 = 5”/Easy
  • Trt12 = 5”/Difficult
  • Trt13 = 5”/Normal
  • Trt21 = 7”/Easy
  • Trt22 = 7”/Difficult
  • Trt23 = 7”/Normal ## Chuẩn bị dữ liệu Đầu tiên, cần chuyển dữ liệu từ dạng rộng (wide) sang dạng dài (long) để phân tích ANOVA.
# Đọc dữ liệu
data <- read.csv("mobile_3d_tasktime_multi.csv", header = TRUE)

# Chuyển đổi từ dạng rộng sang dạng dài
library(tidyr)
library(dplyr)
data_long <- pivot_longer(data, 
                         cols = -SubjectID, 
                         names_to = "Treatment", 
                         values_to = "Time")

# Tách nhân tố từ biến Treatment
data_long$ScreenSize <- ifelse(substr(data_long$Treatment, 4, 4) == "1", "5-inch", "7-inch")

data_long$Interface <- case_when(
  substr(data_long$Treatment, 5, 5) == "1" ~ "Easy",
  substr(data_long$Treatment, 5, 5) == "2" ~ "Difficult",
  substr(data_long$Treatment, 5, 5) == "3" ~ "Normal"
)

# Chuyển các biến nhân tố thành factor
data_long$ScreenSize <- factor(data_long$ScreenSize)
data_long$Interface <- factor(data_long$Interface, levels = c("Easy", "Difficult", "Normal"))
data_long$SubjectID <- factor(data_long$SubjectID)

Kiểm tra dữ liệu và giả định thống kê

# Thống kê mô tả
summary_stats <- data_long %>%
  group_by(ScreenSize, Interface) %>%
  summarise(
    Mean = mean(Time),
    Median = median(Time),
    SD = sd(Time),
    Min = min(Time),
    Max = max(Time),
    n = n()
  )
summary_stats
## # A tibble: 6 × 8
## # Groups:   ScreenSize [2]
##   ScreenSize Interface  Mean Median    SD   Min   Max     n
##   <fct>      <fct>     <dbl>  <dbl> <dbl> <dbl> <dbl> <int>
## 1 5-inch     Easy       26.4   27.8 11.2  1.45   49.4    30
## 2 5-inch     Difficult  29.2   30.4 20.3  0.381  72.5    30
## 3 5-inch     Normal     21.3   22.4 17.8  0.400  51.8    30
## 4 7-inch     Easy       13.8   13.1  8.96 0.680  35.0    30
## 5 7-inch     Difficult  36.5   17.0 43.1  0.614 176.     30
## 6 7-inch     Normal     29.4   32.1 15.7  2.35   64.3    30
# Vẽ đồ thị hộp (boxplot) để kiểm tra phân phối
library(ggplot2)
ggplot(data_long, aes(x = Interface, y = Time, fill = ScreenSize)) +
  geom_boxplot() +
  facet_grid(. ~ ScreenSize) +
  labs(title = "Thời gian thực hiện nhiệm vụ theo kích thước màn hình và giao diện",
       y = "Thời gian (giây)") +
  theme_bw()

So sánh giữa kích thước màn hình

Màn hình 5-inch

  • Thời gian trung vị cao hơn ở giao diện Easy (~30 giây)
  • Biểu hiện nhất quán hơn với ít giá trị ngoại lai
  • Phạm vi phân bố tương đối thu hẹp ở giao diện Normal

Màn hình 7-inch

  • Hiệu suất tốt nhất với giao diện Easy (trung vị ~15 giây)
  • Xuất hiện nhiều giá trị ngoại lai cực đoan ở giao diện Difficult (>100 giây)
  • Độ phân tán lớn hơn đáng kể, đặc biệt ở giao diện Difficult

Phân tích theo loại giao diện

Giao diện Easy

  • Hiệu suất tốt hơn rõ rệt trên màn hình 7-inch (trung vị thấp hơn ~50%)
  • Độ phân tán (IQR) hẹp hơn trên màn hình 7-inch, cho thấy sự ổn định
  • Không xuất hiện giá trị ngoại lai đáng kể ở cả hai kích thước màn hình

Giao diện Difficult

  • Gây ra thách thức lớn nhất cho người dùng ở cả hai kích thước màn hình
  • Màn hình 7-inch xuất hiện nhiều giá trị ngoại lai cực đoan (lên đến ~175 giây)
  • Độ phân tán lớn ở cả hai kích thước, thể hiện qua hộp rộng và râu dài

Giả thuyết của kiểm định Levene

Giả thuyết không (H₀): Phương sai của các nhóm là bằng nhau (có tính đồng nhất phương sai) H₀: σ²₁ = σ²₂ = σ²₃ = σ²₄ = σ²₅ = σ²₆

Giả thuyết thay thế (H₁): Có ít nhất một nhóm có phương sai khác biệt so với các nhóm còn lại (không có tính đồng nhất phương sai) H₁: σ²ᵢ ≠ σ²ⱼ cho ít nhất một cặp i,j

# Kiểm tra tính chuẩn
shapiro_test <- data_long %>%
  group_by(ScreenSize, Interface) %>%
  summarise(p_value = shapiro.test(Time)$p.value)
## `summarise()` has grouped output by 'ScreenSize'. You can override using the
## `.groups` argument.
# Kiểm tra tính đồng nhất phương sai
library(car)
levene_test <- leveneTest(Time ~ ScreenSize * Interface, data = data_long)
levene_test
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value  Pr(>F)    
## group   5  5.4502 0.00011 ***
##       174                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# a) Kiểm định sự khác biệt giữa 6 phương thức với mức ý nghĩa alpha = 0.05
# Thực hiện ANOVA  với phép đo lặp
library(ez)

# Chuyển đổi chính xác cho ezANOVA
data_long$SubjectID <- as.factor(data_long$SubjectID)

anova_result <- ezANOVA(
  data = data_long,
  dv = .(Time),
  wid = .(SubjectID),
  within = .(ScreenSize, Interface),
  detailed = TRUE,
  type = 3
)

print(anova_result)
## $ANOVA
##                 Effect DFn DFd          SSn       SSd          F            p
## 1          (Intercept)   1  29 122768.42235 55641.736 63.9858590 8.030515e-09
## 2           ScreenSize   1  29     35.02228  4467.768  0.2273274 6.370859e-01
## 3            Interface   2  58   4932.88239 16326.815  8.7618796 4.731376e-04
## 4 ScreenSize:Interface   2  58   4130.05899 11843.620 10.1127621 1.707284e-04
##   p<.05          ges
## 1     * 0.5817075360
## 2       0.0003965611
## 3     * 0.0529206424
## 4     * 0.0446927721
## 
## $`Mauchly's Test for Sphericity`
##                 Effect         W            p p<.05
## 3            Interface 0.2362905 1.691402e-09     *
## 4 ScreenSize:Interface 0.2454200 2.875664e-09     *
## 
## $`Sphericity Corrections`
##                 Effect       GGe       p[GG] p[GG]<.05       HFe       p[HF]
## 3            Interface 0.5669868 0.004302888         * 0.5745205 0.004139415
## 4 ScreenSize:Interface 0.5699369 0.002280992         * 0.5778189 0.002174270
##   p[HF]<.05
## 3         *
## 4         *

Kết luận:

Dựa trên kết quả ANOVA bạn cung cấp, chúng ta có thể kiểm định sự khác biệt giữa các phương thức (được xác định bởi yếu tố Interface và sự tương tác giữa ScreenSize và Interface) với mức ý nghĩa alpha = 0.05.

Diễn giải kết quả:

  • ScreenSize (Kích thước màn hình):

    • Giả thuyết:

      • H0​: Không có sự khác biệt đáng kể về Time (biến phụ thuộc) giữa các mức của ScreenSize.

      • H1​: Có sự khác biệt đáng kể về Time giữa ít nhất hai mức của ScreenSize.

    • Kết quả: F(1,29)=0.2273274, p=0.6370859.

    • Kết luận: Vì p>0.05, chúng ta không bác bỏ giả thuyết H0​. Điều này có nghĩa là không có bằng chứng thống kê cho thấy ScreenSize có ảnh hưởng đáng kể đến Time. Kích thước hiệu ứng (ges) cũng rất nhỏ (0.0003965611)

    Interface (Giao diện):

    • Giả thuyết:

      • H0​: Không có sự khác biệt đáng kể về Time giữa các mức của Interface.

      • H1​: Có sự khác biệt đáng kể về Time giữa ít nhất hai mức của Interface.

    • Kiểm định Mauchly cho tính cầu (Sphericity): Kết quả cho thấy W=0.2362905, p=1.691402e−09 (p<0.05). Điều này có nghĩa là giả định về tính cầu đã bị vi phạm. Do đó, chúng ta cần sử dụng các giá trị p đã được hiệu chỉnh (Greenhouse-Geisser hoặc Huynh-Feldt).

    • Kết quả (đã hiệu chỉnh):

      • Sử dụng Greenhouse-Geisser (GG): p[GG]=0.004302888

      • Sử dụng Huynh-Feldt (HF): p[HF]=0.004139415.

    • Kết luận: Vì cả hai giá trị p đã hiệu chỉnh đều nhỏ hơn 0.05 (p[GG]<0.05 và p[HF]<0.05), chúng ta bác bỏ giả thuyết H0​. Điều này có nghĩa là có bằng chứng thống kê cho thấy Interface có ảnh hưởng đáng kể đến Time. Kích thước hiệu ứng (ges) là 0.0529206424, cho thấy một hiệu ứng ở mức vừa phải.

Hiệu ứng tương tác (Interaction Effect):

  • ScreenSize:Interface (Tương tác giữa Kích thước màn hình và Giao diện):

    • Giả thuyết:

      • H0​: Ảnh hưởng của ScreenSize lên Time là như nhau ở tất cả các mức của Interface (và ngược lại, ảnh hưởng của Interface lên Time là như nhau ở tất cả các mức của ScreenSize). Không có sự tương tác.

      • H1​: Có sự tương tác giữa ScreenSize và Interface đối với Time.

    • Kiểm định Mauchly cho tính cầu (Sphericity): Kết quả cho thấy W=0.2454200, p=2.875664e−09 (p<0.05). Giả định về tính cầu cũng bị vi phạm cho hiệu ứng tương tác. Chúng ta cần sử dụng các giá trị p đã được hiệu chỉnh.

    • Kết quả (đã hiệu chỉnh):

      • Sử dụng Greenhouse-Geisser (GG): p[GG]=0.002280992.

      • Sử dụng (HF): p[HF]=0.002174270.

    • Kết luận: Vì cả hai giá trị p đã hiệu chỉnh đều nhỏ hơn 0.05 (p[GG]<0.05 và p[HF]<0.05), chúng ta bác bỏ giả thuyết H0​. Điều này có nghĩa là có bằng chứng thống kê cho thấy có sự tương tác đáng kể giữa ScreenSize và Interface đối với Time. Kích thước hiệu ứng (ges) là 0.0446927721, cho thấy một hiệu ứng ở mức vừa phải.

Tóm tắt kết luận chung:

  • Không có bằng chứng cho thấy kích thước màn hình (ScreenSize) đơn lẻ ảnh hưởng đến thời gian thực hiện (Time).

  • Loại giao diện (Interface) có ảnh hưởng đáng kể đến thời gian thực hiện (Time).

  • Quan trọng hơn, có một sự tương tác đáng kể giữa kích thước màn hình (ScreenSize) và loại giao diện (Interface). Điều này có nghĩa là cách mà loại giao diện ảnh hưởng đến thời gian thực hiện có thể khác nhau tùy thuộc vào kích thước màn hình đang được sử dụng, và ngược lại.

b ) Tìm khoảng tin cậy đồng thời 95%

# Sử dụng mô hình tuyến tính hỗn hợp (linear mixed model)
library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(emmeans)
## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:GGally':
## 
##     pigs
lmm <- lmer(Time ~ ScreenSize * Interface + (1|SubjectID), data = data_long)
summary(lmm)
## Linear mixed model fit by REML ['lmerMod']
## Formula: Time ~ ScreenSize * Interface + (1 | SubjectID)
##    Data: data_long
## 
## REML criterion at convergence: 1518.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6591 -0.4767 -0.0332  0.3846  7.0307 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  SubjectID (Intercept) 282.3    16.8    
##  Residual              225.1    15.0    
## Number of obs: 180, groups:  SubjectID, 30
## 
## Fixed effects:
##                                     Estimate Std. Error t value
## (Intercept)                           26.447      4.112   6.431
## ScreenSize7-inch                     -12.659      3.874  -3.268
## InterfaceDifficult                     2.793      3.874   0.721
## InterfaceNormal                       -5.110      3.874  -1.319
## ScreenSize7-inch:InterfaceDifficult   19.925      5.478   3.637
## ScreenSize7-inch:InterfaceNormal      20.698      5.478   3.778
## 
## Correlation of Fixed Effects:
##             (Intr) ScrS7- IntrfD IntrfN SS7-:ID
## ScrnSz7-nch -0.471                             
## IntrfcDffcl -0.471  0.500                      
## InterfcNrml -0.471  0.500  0.500               
## ScrnSz7-:ID  0.333 -0.707 -0.707 -0.354        
## ScrnSz7-:IN  0.333 -0.707 -0.354 -0.707  0.500
#(i) Khoảng tin cậy đồng thời 95% cho 5" so với 7"
emmeans_size <- emmeans(lmm, ~ ScreenSize)
## NOTE: Results may be misleading due to involvement in interactions
conf_int_size <- confint(emmeans_size)
print(conf_int_size)
##  ScreenSize emmean   SE   df lower.CL upper.CL
##  5-inch       25.7 3.45 36.1     18.7     32.7
##  7-inch       26.6 3.45 36.1     19.6     33.6
## 
## Results are averaged over the levels of: Interface 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95
# So sánh 5" vs 7"
size_contrast <- contrast(emmeans_size, method = "pairwise", adjust = "bonferroni")
print(size_contrast)
##  contrast            estimate   SE  df t.ratio p.value
##  (5-inch) - (7-inch)   -0.882 2.24 145  -0.394  0.6938
## 
## Results are averaged over the levels of: Interface 
## Degrees-of-freedom method: kenward-roger
#Phân tích chi tiết hơn: điều chỉnh theo từng nhóm
# Khoảng tin cậy cho 5" vs 7" ở mỗi loại giao diện
emmeans_size_by_interface <- emmeans(lmm, ~ ScreenSize | Interface)
size_by_interface_contrast <- contrast(emmeans_size_by_interface, method = "pairwise", adjust = "bonferroni")
print(size_by_interface_contrast)
## Interface = Easy:
##  contrast            estimate   SE  df t.ratio p.value
##  (5-inch) - (7-inch)    12.66 3.87 145   3.268  0.0014
## 
## Interface = Difficult:
##  contrast            estimate   SE  df t.ratio p.value
##  (5-inch) - (7-inch)    -7.27 3.87 145  -1.876  0.0627
## 
## Interface = Normal:
##  contrast            estimate   SE  df t.ratio p.value
##  (5-inch) - (7-inch)    -8.04 3.87 145  -2.075  0.0397
## 
## Degrees-of-freedom method: kenward-roger
# Khoảng tin cậy cho Difficult vs Easy ở mỗi kích thước màn hình
emmeans_interface_by_size <- emmeans(lmm, ~ Interface | ScreenSize)
interface_by_size_contrast <- contrast(emmeans_interface_by_size, method = "pairwise", adjust = "bonferroni")
print(interface_by_size_contrast)
## ScreenSize = 5-inch:
##  contrast           estimate   SE  df t.ratio p.value
##  Easy - Difficult      -2.79 3.87 145  -0.721  1.0000
##  Easy - Normal          5.11 3.87 145   1.319  0.5677
##  Difficult - Normal     7.90 3.87 145   2.040  0.1295
## 
## ScreenSize = 7-inch:
##  contrast           estimate   SE  df t.ratio p.value
##  Easy - Difficult     -22.72 3.87 145  -5.865  <.0001
##  Easy - Normal        -15.59 3.87 145  -4.024  0.0003
##  Difficult - Normal     7.13 3.87 145   1.841  0.2032
## 
## Degrees-of-freedom method: kenward-roger 
## P value adjustment: bonferroni method for 3 tests
# Vẽ biểu đồ tương tác
interaction_plot <- ggplot(summary_stats, aes(x = Interface, y = Mean, color = ScreenSize, group = ScreenSize)) +
  geom_line(size = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = Mean - SD, ymax = Mean + SD), width = 0.2) +
  labs(title = "Tương tác giữa kích thước màn hình và loại giao diện",
       y = "Thời gian trung bình (giây)",
       x = "Loại giao diện") +
  theme_bw() +
  theme(legend.position = "bottom")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
print(interaction_plot)

Bài 4.10

câu (a) Thành lập bảng MANOVA một nhân tố

câu (b) Xây dựng khoảng tin cậy đồng thời 95% cho các trung bình tổng thể qua ba thời kỳ

câu (c) Kiểm tra các điều kiện cho MANOVA có thỏa với dữ liệu này không ? Giải thích ?

Mô tả dữ liệu bài 4.10

Dữ liệu trong bài toán 4.10 bao gồm các đặc trưng hình thái của hộp sọ người được thu thập từ ba thời kỳ lịch sử khác nhau. Mỗi mẫu hộp sọ được mô tả bởi bốn biến định lượng như sau:

  • \(X_1\): Chiều rộng hộp sọ (Cranial Breadth)
  • \(X_2\): Chiều cao hộp sọ (Cranial Height)
  • \(X_3\): Chiều dài cơ bản hộp sọ (Basial Length)
  • \(X_4\): Chiều cao cơ bản hộp sọ (Basial Height)

Ngoài ra, mỗi mẫu còn được gán với một biến phân loại:

  • Group: Biến định danh cho thời kỳ khảo cổ học (gồm 3 mức: 1, 2 và 3)
  • Số biến phụ thuộc: 4
  • Số nhóm (nhân tố): 3 thời kỳ
  • Số chiều của dữ liệu: 4 chiều (vì mỗi cá thể được biểu diễn bằng 4 đặc trưng)
  • Loại phân tích: One-way MANOVA (với nhân tố là thời kỳ Group)

câu (a) Thành lập bảng MANOVA một nhân tố

Ta kiểm định xem có sự khác biệt thống kê trong tổ hợp trung bình của 4 đặc trưng hình thái hộp sọ giữa ba thời kỳ khảo sát hay không

\[ \begin{aligned} H_0 &: \mu_1 = \mu_2 = \mu_3 \\ H_1 &: \text{Tồn tại ít nhất một } \mu_i \ne \mu_j \end{aligned} \]

  • \(\mu_i\): vector trung bình 4 biến hình thái hộp sọ (X1–X4) tại thời kỳ thứ \(i\)
  • \(H_0\): Không có sự khác biệt nào giữa tổ hợp trung bình của 4 biến giữa 3 thời kỳ khảo sát
  • \(H_1\): Có ít nhất một thời kỳ có tổ hợp trung bình khác với thời kỳ còn lại
#  Đọc dữ liệu
data <- read.table("exo4-10.dat", header = FALSE)
colnames(data) <- c("X1", "X2", "X3", "X4", "Group")
data$Group <- factor(data$Group)

Đọc dữ liệu từ file .dat, không có dòng tiêu đề (header = FALSE)

Gán tên biến: - X1 → Chiều rộng hộp sọ

  • X2 → Chiều cao hộp sọ

  • X3 → Chiều dài cơ bản hộp sọ

  • X4 → Chiều cao cơ bản hộp sọ

  • Group → Thời kỳ khảo sát (1, 2, 3)

Chuyển Group thành kiểu phân loại (factor) để dùng cho mô hình

# kiểm tra chuẩn Shapiro-Wilk từng biến
cat("Kiểm tra chuẩn (Shapiro-Wilk) từng biến trong từng nhóm:\n")
## Kiểm tra chuẩn (Shapiro-Wilk) từng biến trong từng nhóm:
for (g in levels(data$Group)) {
  cat("\nThời kỳ", g, ":\n")
  group_data <- subset(data, Group == g)[, 1:4]
  
  for (j in 1:4) {
    shapiro_result <- shapiro.test(group_data[, j])
    cat(" - Biến", colnames(group_data)[j], 
        ": p-value =", round(shapiro_result$p.value, 4), "\n")
  }
}
## 
## Thời kỳ 1 :
##  - Biến X1 : p-value = 0.8603 
##  - Biến X2 : p-value = 0.2536 
##  - Biến X3 : p-value = 0.6282 
##  - Biến X4 : p-value = 0.6772 
## 
## Thời kỳ 2 :
##  - Biến X1 : p-value = 0.0403 
##  - Biến X2 : p-value = 0.5799 
##  - Biến X3 : p-value = 0.7495 
##  - Biến X4 : p-value = 0.3762 
## 
## Thời kỳ 3 :
##  - Biến X1 : p-value = 0.0314 
##  - Biến X2 : p-value = 0.9031 
##  - Biến X3 : p-value = 0.8849 
##  - Biến X4 : p-value = 0.3881

Ta kiểm tra từng biến đơn lẻ trong từng nhóm bằng Shapiro-Wilk.

Vòng lặp ngoài: duyệt qua từng nhóm (Group 1 → 3)

Vòng lặp trong: duyệt qua từng biến X1 → X4

Hàm shapiro.test() kiểm tra giả thuyết:

\[ \begin{aligned} H_0 &: \text{Biến có phân phối chuẩn} \\ H_1 &: \text{Biến không chuẩn} \end{aligned} \]

Nếu p-value > 0.05 → chấp nhận giả định chuẩn

Hầu hết các biến đều thỏa mãn giả định phân phối chuẩn đơn biến trong từng nhóm.

Chỉ có biến X1 trong thời kỳ 2 và thời kỳ 3 vi phạm chuẩn đơn biến (p-value < 0.05).

Tuy nhiên, không có nhóm nào vi phạm chuẩn ở nhiều biến cùng lúc, và Shapiro-Wilk tương đối nhạy khi n nhỏ nên có thể chấp nhận được khi kiểm định.

# Thực hiện kiểm định MANOVA 
manova_model <- manova(cbind(X1, X2, X3, X4) ~ Group, data = data)
summary_manova <- summary(manova_model, test = "Wilks")

cat("\nKết quả MANOVA (Wilks’ Lambda):\n")
## 
## Kết quả MANOVA (Wilks’ Lambda):
print(summary_manova)
##           Df  Wilks approx F num Df den Df  Pr(>F)  
## Group      2 0.8301   2.0491      8    168 0.04358 *
## Residuals 87                                        
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  • Hàm manova() tạo mô hình phân tích phương sai đa biến:

\[ H_0: \mu_1 = \mu_2 = \mu_3 \]

  • summary(..., test = "Wilks") : sử dụng thống kê Wilks’ Lambda (phổ biến nhất)

  • print() in bảng kết quả: có thống kê Wilks, F-approx, df và p-value

Vì p-value = 0.0436 < 0.05, ta bác bỏ giả thuyết không \(H_0\) ở mức ý nghĩa 5%.

⟹ Kết luận: Có sự khác biệt có ý nghĩa thống kê giữa ít nhất hai trong ba thời kỳ khảo sát về tổ hợp trung bình các đặc

trưng hình thái hộp sọ (bao gồm chiều rộng, chiều cao, chiều dài cơ bản và chiều cao cơ bản của hộp sọ).

câu (b) Xây dựng khoảng tin cậy đồng thời 95% cho các trung bình tổng thể qua ba thời kỳ

Sử dụng công thức Hotelling’s \(T^2\) để xây dựng khoảng tin cậy đồng thời 95% cho từng biến trong mỗi nhóm, kiểm soát xác suất sai toàn phần là 0.05.

# Đọc dữ liệu
data <- read.table("exo4-10.dat", header = FALSE)
colnames(data) <- c("X1", "X2", "X3", "X4", "Group")
data$Group <- factor(data$Group)

Gán tên các cột: 4 biến định lượng và Group (nhóm thời kỳ)

Chuyển Group thành factor để xử lý theo nhóm

# Thông số
p <- 4  # số biến
g <- length(unique(data$Group))  # số nhóm
n_i <- table(data$Group)
N <- nrow(data)
alpha <- 0.05

Xác định số biến (p = 4)

Số nhóm (g = 3)

Bảng đếm số lượng mẫu từng nhóm (n_i)

Kích thước mẫu tổng (N)

Mức ý nghĩa để dựng khoảng tin cậy (95%)

# Trung bình theo từng nhóm
group_means <- aggregate(. ~ Group, data = data, mean)

Tính trung bình của từng biến trong từng nhóm

Lưu vào group_means, sẽ dùng làm trung tâm khoảng tin cậy

# Ma trận hiệp phương sai gộp
Spooled <- matrix(0, nrow = p, ncol = p)
for (i in levels(data$Group)) {
  subset_i <- data[data$Group == i, 1:4]
  ni <- nrow(subset_i)
  Si <- cov(subset_i)
  Spooled <- Spooled + (ni - 1) * Si
}
Spooled <- Spooled / (N - g)

Duyệt từng nhóm:

  • Lấy tập con theo nhóm
  • Tính ma trận hiệp phương sai mẫu của nhóm \(S_i\)
  • Cộng vào tổng có trọng số: \(\sum (n_i - 1) S_i\)

Chia cho \(N - g\) để ra ma trận hiệp phương sai gộp \(S_p\)

# Thống kê tới hạn F
F_crit <- qf(1 - alpha, df1 = p, df2 = N - p - g + 1)
c_alpha <- (p * (N - g)) / (N - p - g + 1) * F_crit

Tìm ngưỡng F tới hạn từ phân phối F:

\[ F_{1 - \alpha}(p,\; N - p - g + 1) \]

Tính hệ số hiệu chỉnh:

\[ c_\alpha = \frac{p(N - g)}{N - p - g + 1} \cdot F_{1 - \alpha} \]

Đây là hệ số dùng để điều chỉnh sai số chuẩn theo Hotelling’s \(T^2\)

# Tính khoảng tin cậy cho từng nhóm
for (i in 1:g) {
  cat("\nThời kỳ", i, "\n")
  mean_vec <- as.numeric(group_means[i, 2:5])
  ni <- as.numeric(n_i[i])
  se <- sqrt(diag(Spooled) / ni * c_alpha)
  lower <- mean_vec - se
  upper <- mean_vec + se
  result <- data.frame(
    Biến = colnames(data)[1:4],
    Trung_bình = round(mean_vec, 2),
    CI_thấp = round(lower, 2),
    CI_cao = round(upper, 2)
  )
  print(result)
}
## 
## Thời kỳ 1 
##    Biến Trung_bình CI_thấp CI_cao
## X1   X1     131.37  128.72 134.02
## X2   X2     133.60  130.85 136.35
## X3   X3      99.17   96.26 102.08
## X4   X4      50.53   48.71  52.35
## 
## Thời kỳ 2 
##    Biến Trung_bình CI_thấp CI_cao
## X1   X1     132.37  129.72 135.02
## X2   X2     132.70  129.95 135.45
## X3   X3      99.07   96.16 101.98
## X4   X4      50.23   48.41  52.05
## 
## Thời kỳ 3 
##    Biến Trung_bình CI_thấp CI_cao
## X1   X1     134.47  131.82 137.12
## X2   X2     133.80  131.05 136.55
## X3   X3      96.03   93.12  98.94
## X4   X4      50.57   48.75  52.39

Lặp qua từng thời kỳ (1 → 3):

  • Tính trung bình của nhóm i: mean_vec

  • Tính sai số chuẩn hiệu chỉnh:

\[ SE_j = \sqrt{ \frac{c_\alpha \cdot s_{jj}}{n_i} } \]

  • Tạo khoảng tin cậy:

\[ \left[ \bar{x}_j - SE_j,\; \bar{x}_j + SE_j \right] \]

câu (c) Kiểm tra các điều kiện cho MANOVA có thỏa với dữ liệu này không ? Giải thích ?

# Đọc dữ liệu
data <- read.table("exo4-10.dat", header = FALSE)
colnames(data) <- c("X1", "X2", "X3", "X4", "Group")
data$Group <- factor(data$Group)

Gán tên cho các cột

Chuyển Group thành biến phân loại (factor) để dùng đúng trong phân tích theo nhóm

#  Kiểm tra chuẩn từng biến bằng Shapiro-Wilk
cat("Kiểm tra chuẩn (Shapiro-Wilk) từng biến trong từng nhóm:\n")
## Kiểm tra chuẩn (Shapiro-Wilk) từng biến trong từng nhóm:
for (g in levels(data$Group)) {
  cat("\nThời kỳ", g, ":\n")
  group_data <- subset(data, Group == g)[, 1:4]
  for (j in 1:4) {
    shapiro_result <- shapiro.test(group_data[, j])
    cat(" - Biến", colnames(group_data)[j], 
        ": p-value =", round(shapiro_result$p.value, 4), "\n")
  }
}
## 
## Thời kỳ 1 :
##  - Biến X1 : p-value = 0.8603 
##  - Biến X2 : p-value = 0.2536 
##  - Biến X3 : p-value = 0.6282 
##  - Biến X4 : p-value = 0.6772 
## 
## Thời kỳ 2 :
##  - Biến X1 : p-value = 0.0403 
##  - Biến X2 : p-value = 0.5799 
##  - Biến X3 : p-value = 0.7495 
##  - Biến X4 : p-value = 0.3762 
## 
## Thời kỳ 3 :
##  - Biến X1 : p-value = 0.0314 
##  - Biến X2 : p-value = 0.9031 
##  - Biến X3 : p-value = 0.8849 
##  - Biến X4 : p-value = 0.3881

Duyệt qua từng nhóm (Group 1 → 3):

Trong mỗi nhóm, kiểm định từng biến \(X_1\) đến \(X_4\) xem có phân phối chuẩn không bằng Shapiro-Wilk.

Hàm shapiro.test() trả về p-value:

  • Nếu p > 0.05: chấp nhận \(H_0\), biến có thể xem là chuẩn
  • Nếu p < 0.05: bác bỏ \(H_0\), biến không chuẩn

Nhận xét :

Ta đã kiểm định từng biến bằng Shapiro-Wilk trong từng nhóm.

Chỉ có biến X1 ở nhóm 2 và nhóm 3 có p-value < 0.05, tức là không hoàn toàn chuẩn đơn biến.

Thấy rằng :

Mỗi nhóm có n = 30 (tương đối lớn).

Vi phạm chỉ xảy ra ở một biến trong mỗi nhóm.

→ Giả định chuẩn đa biến vẫn được xem là chấp nhận được trong thực hành.

# So sánh ma trận hiệp phương sai từng nhóm
cat("\n So sánh ma trận hiệp phương sai của từng nhóm:\n")
## 
##  So sánh ma trận hiệp phương sai của từng nhóm:
for (g in levels(data$Group)) {
  group_data <- subset(data, Group == g)[, 1:4]
  cat("\nCovariance matrix of group", g, ":\n")
  print(round(cov(group_data), 2))
}
## 
## Covariance matrix of group 1 :
##       X1    X2    X3    X4
## X1 26.31  4.15  0.45  7.25
## X2  4.15 19.97 -0.79  0.39
## X3  0.45 -0.79 34.63 -1.92
## X4  7.25  0.39 -1.92  7.64
## 
## Covariance matrix of group 2 :
##       X1    X2    X3   X4
## X1 23.14  1.01  4.77 1.84
## X2  1.01 21.60  3.37 5.62
## X3  4.77  3.37 18.89 0.19
## X4  1.84  5.62  0.19 8.74
## 
## Covariance matrix of group 3 :
##       X1    X2    X3    X4
## X1 12.12  0.79 -0.77  0.90
## X2  0.79 24.79  3.59 -0.09
## X3 -0.77  3.59 20.72  1.67
## X4  0.90 -0.09  1.67 12.60

  • Tính và in ra ma trận hiệp phương sai (covariance matrix) của 4 biến trong từng nhóm

  • Mục đích: so sánh bằng mắt xem các ma trận có khác nhau rõ rệt không

  • Nếu các ma trận tương đồng → thỏa giả định đồng nhất phương sai

Kết luận :

Qua việc so sánh ma trận hiệp phương sai của từng nhóm, ta nhận thấy các ma trận này có giá trị phương sai và tương quan tương đối

tương đồng. Mặc dù có một vài khác biệt nhỏ về độ lớn, nhưng chúng không quá chênh lệch và không có sự sai khác vượt mức cho phép.

Vì vậy, có thể kết luận rằng giả định đồng nhất phương sai giữa các nhóm được thỏa mãn ở mức chấp nhận được, và không ảnh hưởng đáng

kể đến độ tin cậy của kiểm định MANOVA.

# Kiểm tra kích thước mẫu và phân loại nhóm
cat("\n Kiểm tra số mẫu và cấu trúc nhóm:\n")
## 
##  Kiểm tra số mẫu và cấu trúc nhóm:
print(table(data$Group))
## 
##  1  2  3 
## 30 30 30
cat("\nTổng số mẫu:", nrow(data), "\n")
## 
## Tổng số mẫu: 90
cat("Số biến:", ncol(data) - 1, "\n")
## Số biến: 4

table(data$Group): đếm số mẫu trong từng nhóm

nrow(data): tổng số quan sát toàn bộ

ncol(data) - 1: số biến định lượng (trừ cột Group)

Kiểm tra xem:

Mỗi nhóm có số mẫu đủ lớn (thường ≥ 30) để đảm bảo ổn định.

Nhận xét :

Các nhóm trong bộ dữ liệu có số lượng mẫu bằng nhau và đều có 30 mẫu, đảm bảo độ ổn định

cho kiểm định MANOVA, dù cho giả định chuẩn có thể bị vi phạm nhẹ.

Ngoài ra, cấu trúc của dữ liệu cũng đúng định dạng: 4 biến định lượng và 1 biến nhóm phân loại.

⟹ Dữ liệu đáp ứng tốt yêu cầu về thiết kế và kích thước mẫu để thực hiện MANOVA.

Kết luận :

Tuy có một vài vi phạm nhẹ về giả định phân phối chuẩn đơn biến nhưng do kích thước mẫu đủ lớn và các vi phạm không nghiêm trọng nên các điều kiện cho MANOVA thỏa với dữ liệu này.