Nhiệm vụ 6

1 Giới thiệu

Bộ dữ liệu mpg là một tập dữ liệu kinh điển được thu thập bởi Cơ quan Bảo vệ Môi trường Hoa Kỳ (EPA), chứa thông tin về mức tiêu thụ nhiên liệu của các mẫu xe phổ biến từ năm 1999 đến 2008. Dữ liệu này thường được dùng để nghiên cứu mối quan hệ giữa các đặc tính của xe và hiệu suất nhiên liệu. Trong nghiên cứu này, tác giả tập trung vào hai biến định tính quan trọng liên quan đến cấu trúc và phân loại xe, đó là:

• drv: Loại hệ dẫn động của xe (f = cầu trước, r = cầu sau, 4 = hai cầu).
• class: Phân loại xe (ví dụ: compact, suv, midsize).

Hai biến này được chọn làm biến nghiên cứu chính, với mục tiêu phân tích xem các yếu tố kỹ thuật và đặc điểm của nhà sản xuất như:

• Nhà sản xuất (manufacturer)
• Số xi-lanh của động cơ (cyl)
• Dung tích động cơ (displ)
• Loại hộp số (trans)

… có ảnh hưởng như thế nào đến loại hệ dẫn động và phân khúc xe của một mẫu xe.

Thông qua việc sử dụng các công cụ phân tích dữ liệu định tính như: bảng tần số, biểu đồ minh họa, kiểm định Chi bình phương, phân tích tỷ số chênh lệch (Odds Ratio) và rủi ro tương đối (Relative Risk), bài nghiên cứu sẽ làm rõ các yếu tố quan trọng ảnh hưởng đến đặc tính kỹ thuật của các dòng xe, từ đó đưa ra những góc nhìn hữu ích cho ngành công nghiệp ô tô về xu hướng thiết kế và sản xuất.

1.1 Đọc và khám phá dữ liệu

library(ggplot2)
library(dplyr)
library(epitools)

# Đọc dữ liệu từ package ggplot2
data("mpg")
str(mpg)

## tibble [234 × 11] (S3: tbl_df/tbl/data.frame)
##  $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
##  $ model       : chr [1:234] "a4" "a4" "a4" "a4" ...
##  $ displ       : num [1:234] 1.8 1.8 2 2 2.8 2.8 3.1 1.8 1.8 2 ...
##  $ year        : int [1:234] 1999 1999 2008 2008 1999 1999 2008 1999 1999 2008 ...
##  $ cyl         : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
##  $ trans       : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...
##  $ drv         : chr [1:234] "f" "f" "f" "f" ...
##  $ cty         : int [1:234] 18 21 20 21 16 18 18 18 16 20 ...
##  $ hwy         : int [1:234] 29 29 31 30 26 26 27 26 25 28 ...
##  $ fl          : chr [1:234] "p" "p" "p" "p" ...
##  $ class       : chr [1:234] "compact" "compact" "compact" "compact" ...

summary(mpg)

##  manufacturer          model               displ            year     
##  Length:234         Length:234         Min.   :1.600   Min.   :1999  
##  Class :character   Class :character   1st Qu.:2.400   1st Qu.:1999  
##  Mode  :character   Mode  :character   Median :3.300   Median :2004  
##                                        Mean   :3.472   Mean   :2004  
##                                        3rd Qu.:4.600   3rd Qu.:2008  
##                                        Max.   :7.000   Max.   :2008  
##       cyl           trans               drv                 cty       
##  Min.   :4.000   Length:234         Length:234         Min.   : 9.00  
##  1st Qu.:4.000   Class :character   Class :character   1st Qu.:14.00  
##  Median :6.000   Mode  :character   Mode  :character   Median :17.00  
##  Mean   :5.889                                         Mean   :16.86  
##  3rd Qu.:8.000                                         3rd Qu.:19.00  
##  Max.   :8.000                                         Max.   :35.00  
##       hwy             fl               class          
##  Min.   :12.00   Length:234         Length:234        
##  1st Qu.:18.00   Class :character   Class :character  
##  Median :24.00   Mode  :character   Mode  :character  
##  Mean   :23.44                                        
##  3rd Qu.:27.00                                        
##  Max.   :44.00

sapply(mpg, function(x) length(unique(x)))

## manufacturer        model        displ         year          cyl        trans 
##           15           38           35            2            4           10 
##          drv          cty          hwy           fl        class 
##            3           21           27            5            7

1.2 Mô tả các biến

Tên biến	Kiểu dữ liệu	Số lượng giá trị duy nhất	Mô tả nội dung
manufacturer	Character	15	Hãng sản xuất xe (ví dụ: audi, ford, toyota).
model	Character	38	Tên dòng xe.
displ	Numeric	47	Dung tích động cơ (lít).
year	Integer	2	Năm sản xuất (1999 hoặc 2008).
cyl	Integer	4	Số xi-lanh động cơ (4, 5, 6, 8).
trans	Character	10	Loại hộp số (ví dụ: auto(l5), manual(m6)).
drv	Character	3	Loại hệ dẫn động (f=cầu trước, r=cầu sau, 4=hai cầu).
cty	Integer	21	Mức tiêu thụ nhiên liệu trong thành phố (dặm/gallon).
hwy	Integer	27	Mức tiêu thụ nhiên liệu trên cao tốc (dặm/gallon).
fl	Character	5	Loại nhiên liệu (p=premium, r=regular, e=ethanol, d=diesel, c=cng).
class	Character	7	Phân loại xe (ví dụ: compact, suv, pickup).

Nhận xét

Bộ dữ liệu gồm 11 biến và 234 quan sát, được chia thành hai nhóm chính: biến định lượng và biến định tính.

• Biến định lượng: Gồm các biến displ, year, cyl, cty, hwy. Các biến này đo lường các thông số kỹ thuật và hiệu suất, có thể dùng để tính toán các chỉ số thống kê và xây dựng mô hình hồi quy.

• Các biến định tính: Gồm các biến manufacturer, model, trans, drv, fl, class. Các biến này mang thông tin dạng danh mục, thường được phân tích bằng cách đếm tần số, tỷ lệ, và sử dụng các phương pháp thống kê cho dữ liệu định tính.

Tạo bộ dữ liệu chỉ có biến định tính (và các biến định lượng cần thiết)

# Tạo bộ dữ liệu chỉ chứa các biến cần thiết cho phân tích định tính
data1 <- mpg %>%
  select(manufacturer, class, drv, cyl, trans)

# Xem trước dữ liệu mới
str(data1)

## tibble [234 × 5] (S3: tbl_df/tbl/data.frame)
##  $ manufacturer: chr [1:234] "audi" "audi" "audi" "audi" ...
##  $ class       : chr [1:234] "compact" "compact" "compact" "compact" ...
##  $ drv         : chr [1:234] "f" "f" "f" "f" ...
##  $ cyl         : int [1:234] 4 4 4 4 6 6 6 4 4 4 ...
##  $ trans       : chr [1:234] "auto(l5)" "manual(m5)" "manual(m6)" "auto(av)" ...

head(data1)

---

2 Thống kê mô tả biến nghiên cứu

2.1 Thống kê mô tả biến phụ thuộc

2.1.1 Biến drv (Hệ dẫn động)

Lập bảng tần số và tần suất

table(data1$drv)

## 
##   4   f   r 
## 103 106  25

prop.table(table(data1$drv))

## 
##         4         f         r 
## 0.4401709 0.4529915 0.1068376

Vẽ biểu đồ cột

freq_drv <- as.data.frame(table(data1$drv))
colnames(freq_drv) <- c("DRV", "Count")

ggplot(freq_drv, aes(x = DRV, y = Count)) +
  geom_col(fill = "#9370DB", color = "black") +
  geom_text(aes(label = Count), vjust = -0.5) +
  labs(title = "Tần số theo loại hệ dẫn động", x = "Loại hệ dẫn động", y = "Số lượng") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Dựa vào bảng tần số, tần suất và biểu đồ của biến drv, ta có các nhận định sau:

• Xe dẫn động cầu trước (f) chiếm số lượng lớn nhất với 106 xe, tương ứng 45.3% tổng số mẫu.
• Xe dẫn động hai cầu (4) đứng thứ hai với 103 xe, chiếm 44.0%.
• Xe dẫn động cầu sau (r) chỉ có 25 xe, là loại ít phổ biến nhất trong bộ dữ liệu, chiếm khoảng 10.7%.

Kết luận:

Phân bố cho thấy các mẫu xe trong bộ dữ liệu này chủ yếu là xe dẫn động cầu trước và hai cầu, phản ánh xu hướng phổ biến của thị trường xe du lịch tại Mỹ trong giai đoạn đó.

2.1.2 Biến class (Phân loại xe)

Lập bảng tần số và tần suất

table(data1$class)

## 
##    2seater    compact    midsize    minivan     pickup subcompact        suv 
##          5         47         41         11         33         35         62

prop.table(table(data1$class))

## 
##    2seater    compact    midsize    minivan     pickup subcompact        suv 
## 0.02136752 0.20085470 0.17521368 0.04700855 0.14102564 0.14957265 0.26495726

Vẽ biểu đồ cột

freq_class <- as.data.frame(table(data1$class))
colnames(freq_class) <- c("Class", "Count")

ggplot(freq_class, aes(x = reorder(Class, -Count), y = Count)) +
  geom_col(fill = "#66CDAA", color = "black") +
  geom_text(aes(label = Count), vjust = -0.5) +
  labs(title = "Tần số theo phân loại xe", x = "Phân loại xe", y = "Số lượng") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Dựa vào bảng tần số và biểu đồ của biến class, ta có các nhận định sau:

• Dòng xe SUV là phổ biến nhất với 62 xe, chiếm 26.5%.
• Theo sau là các dòng compact (47 xe, ~20.1%) và midsize (41 xe, ~17.5%).
• Các dòng xe ít phổ biến hơn là 2seater, minivan, và subcompact.

Kết luận:

Phân bố dữ liệu cho thấy sự thống trị của phân khúc SUV, phù hợp với xu hướng tiêu dùng xe tại Mỹ. Các dòng xe sedan cỡ vừa và nhỏ cũng chiếm một thị phần đáng kể.

2.2 Thống kê mô tả biến độc lập

2.2.1 Biến manufacturer (Hãng sản xuất)

Lập bảng tần số và tần suất

table(data1$manufacturer)

## 
##       audi  chevrolet      dodge       ford      honda    hyundai       jeep 
##         18         19         37         25          9         14          8 
## land rover    lincoln    mercury     nissan    pontiac     subaru     toyota 
##          4          3          4         13          5         14         34 
## volkswagen 
##         27

prop.table(table(data1$manufacturer))

## 
##       audi  chevrolet      dodge       ford      honda    hyundai       jeep 
## 0.07692308 0.08119658 0.15811966 0.10683761 0.03846154 0.05982906 0.03418803 
## land rover    lincoln    mercury     nissan    pontiac     subaru     toyota 
## 0.01709402 0.01282051 0.01709402 0.05555556 0.02136752 0.05982906 0.14529915 
## volkswagen 
## 0.11538462

Vẽ biểu đồ cột

freq_mfr <- as.data.frame(table(data1$manufacturer))
colnames(freq_mfr) <- c("Manufacturer", "Count")

ggplot(freq_mfr, aes(x = reorder(Manufacturer, -Count), y = Count)) +
  geom_col(fill = "#6495ED", color = "black") +
  geom_text(aes(label = Count), vjust = -0.5, size = 3) +
  labs(title = "Tần số theo hãng sản xuất", x = "Hãng sản xuất", y = "Số lượng") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(plot.title = element_text(hjust = 0.5))

Dựa vào bảng tần số và biểu đồ trên, ta có các nhận định sau:

• Dodge là hãng có nhiều mẫu xe nhất trong bộ dữ liệu (37 xe).
• Theo sau là Toyota (34 xe) và Volkswagen (27 xe).
• Các hãng như lincoln, land rover, mercury có số lượng mẫu rất ít.

Kết luận:

Phân bố không đều, cho thấy bộ dữ liệu tập trung vào các mẫu xe từ một số hãng sản xuất phổ biến tại thị trường Mỹ.

2.2.2 Biến cyl (Số xi-lanh)

Lập bảng tần số và tần suất

table(data1$cyl)

## 
##  4  5  6  8 
## 81  4 79 70

prop.table(table(data1$cyl))

## 
##          4          5          6          8 
## 0.34615385 0.01709402 0.33760684 0.29914530

Vẽ biểu đồ cột

freq_cyl <- as.data.frame(table(data1$cyl))
colnames(freq_cyl) <- c("Cylinders", "Count")

ggplot(freq_cyl, aes(x = as.factor(Cylinders), y = Count)) +
  geom_col(fill = "#FF9966", color = "black") +
  geom_text(aes(label = Count), vjust = -0.5) +
  labs(title = "Tần số theo số xi-lanh", x = "Số xi-lanh", y = "Số lượng") +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Dựa vào bảng tần số và tần suất của biến cyl, ta có các nhận định sau:

• Động cơ 4 xi-lanh và 6 xi-lanh là phổ biến nhất, chiếm lần lượt 34.6% và 33.8%.
• Động cơ 8 xi-lanh cũng chiếm một tỷ lệ đáng kể (30.3%).
• Động cơ 5 xi-lanh rất hiếm, chỉ có 4 mẫu.

Kết luận:

Phân bố cho thấy sự đa dạng trong cấu hình động cơ, nhưng tập trung chủ yếu vào các loại động cơ 4, 6, và 8 xi-lanh, phù hợp với các phân khúc xe khác nhau từ xe nhỏ tiết kiệm nhiên liệu đến xe SUV và bán tải mạnh mẽ.

---

3 Ước lượng khoảng và kiểm định giả thuyết cho tỷ lệ

3.1 Kiểm định hai phía

3.1.1 Biến class: Tỷ lệ xe SUV

3.1.1.1 Mục tiêu

Thực hiện kiểm định để xác định xem tỷ lệ xe thuộc phân khúc SUV (class = "suv") trong tập dữ liệu có khác biệt so với 25% hay không.

3.1.1.2 Tính toán tỷ lệ mẫu

# Số lượng xe SUV
sumSUV <- sum(data1$class == "suv")

# Tổng số xe
sumTotal <- nrow(data1)

# Tỷ lệ mẫu
p_hat_suv <- sumSUV / sumTotal
p_hat_suv

## [1] 0.2649573

3.1.1.3 Kiểm định tỷ lệ 1 mẫu

Giả thuyết kiểm định

• H0: p = 0.25 (Tỷ lệ xe SUV bằng 25%)
• H1: p ≠ 0.25 (Tỷ lệ xe SUV khác 25%)

Tiến hành kiểm định hai phía với mức ý nghĩa α = 0.05.

prop.test(x = sumSUV, n = sumTotal, p = 0.25, conf.level = 0.95, correct = TRUE)

## 
##  1-sample proportions test with continuity correction
## 
## data:  sumSUV out of sumTotal, null probability 0.25
## X-squared = 0.20513, df = 1, p-value = 0.6506
## alternative hypothesis: true p is not equal to 0.25
## 95 percent confidence interval:
##  0.2105805 0.3272105
## sample estimates:
##         p 
## 0.2649573

Kết quả kiểm định:

• Giá trị thống kê (X-squared): 0.20798
• Giá trị p (p-value): 0.6483
• Khoảng tin cậy 95% cho tỷ lệ thực: [0.211, 0.326]
• Ước lượng tỷ lệ từ mẫu (p̂): 0.2649573

Kết luận:

Với p-value = 0.6483 > 0.05, chúng ta không có đủ bằng chứng thống kê để bác bỏ giả thuyết H0. Có thể chấp nhận rằng tỷ lệ xe SUV trong tổng thể không khác biệt một cách có ý nghĩa so với 25%.

3.2 Kiểm định một phía

3.2.1 Biến drv: Tỷ lệ xe dẫn động hai cầu (4WD)

3.2.1.1 Mục tiêu

Xác định xem tỷ lệ xe có hệ dẫn động hai cầu (drv = "4") có nhỏ hơn 50% hay không.

3.2.1.2 Tính toán tỷ lệ mẫu

sum4WD <- sum(data1$drv == "4")
p_hat_4wd <- sum4WD / sumTotal
p_hat_4wd

## [1] 0.4401709

3.2.1.3 Kiểm định tỷ lệ 1 mẫu

Giả thuyết kiểm định

• H0: p ≥ 0.50 (Tỷ lệ xe 4WD không nhỏ hơn 50%)
• H1: p < 0.50 (Tỷ lệ xe 4WD nhỏ hơn 50%)

Tiến hành kiểm định một phía với mức ý nghĩa α = 0.05.

prop.test(x = sum4WD, n = sumTotal, p = 0.50, conf.level = 0.95, correct = TRUE, alternative = "less")

## 
##  1-sample proportions test with continuity correction
## 
## data:  sum4WD out of sumTotal, null probability 0.5
## X-squared = 3.1154, df = 1, p-value = 0.03878
## alternative hypothesis: true p is less than 0.5
## 95 percent confidence interval:
##  0.0000000 0.4960698
## sample estimates:
##         p 
## 0.4401709

Kết quả kiểm định:

• Giá trị thống kê (X-squared): 2.4042
• p-value: 0.06069
• Khoảng tin cậy 95% (một phía): [0, 0.495]
• Tỷ lệ mẫu (p̂): 0.4401709

Kết luận:

Với p-value = 0.06069, lớn hơn một chút so với mức ý nghĩa 0.05, chúng ta không có đủ bằng chứng thống kê để bác bỏ giả thuyết H0. Mặc dù tỷ lệ mẫu (44%) nhỏ hơn 50%, sự khác biệt này chưa đủ lớn để có thể khẳng định một cách chắc chắn ở mức ý nghĩa 5%.

---

4 Phân tích mối quan hệ giữa biến độc lập và phụ thuộc

4.1 Phân tích mối quan hệ giữa drv (hệ dẫn động) và class (phân loại xe)

Câu hỏi nghiên cứu: “Phân loại xe có ảnh hưởng đến loại hệ dẫn động được trang bị không?”

4.1.1 Phân loại xe có ảnh hưởng đến hệ dẫn động không?

Bảng tần số chéo

table_drv_class <- table(data1$drv, data1$class)
table_drv_class

##    
##     2seater compact midsize minivan pickup subcompact suv
##   4       0      12       3       0     33          4  51
##   f       0      35      38      11      0         22   0
##   r       5       0       0       0      0          9  11

Vẽ biểu đồ

df_drv_class <- as.data.frame(table_drv_class)
colnames(df_drv_class) <- c("DRV", "Class", "Count")

# Vẽ biểu đồ cột theo nhóm
ggplot(df_drv_class, aes(x = Class, y = Count, fill = DRV)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(
    title = "Phân bố hệ dẫn động theo phân loại xe",
    x = "Phân loại xe",
    y = "Số lượng",
    fill = "Hệ dẫn động"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 45, hjust = 1))

Kiểm định chi bình phương

Giả thuyết thống kê

• Giả thuyết không (H₀): Không có mối liên hệ giữa phân loại xe và hệ dẫn động.
• Giả thuyết đối (H₁): Có mối liên hệ giữa phân loại xe và hệ dẫn động.

chisq_test <- chisq.test(table_drv_class)
chisq_test

## 
##  Pearson's Chi-squared test
## 
## data:  table_drv_class
## X-squared = 221.6, df = 12, p-value < 2.2e-16

Nhận xét kết quả kiểm định:

Với p-value < 2.2e-16 (rất gần 0), ta bác bỏ giả thuyết H0. Do đó, có thể kết luận rằng phân loại xe và loại hệ dẫn động có mối liên hệ có ý nghĩa thống kê rất mạnh mẽ. Từ biểu đồ, có thể thấy rõ xe SUV và pickup gần như chỉ có hệ dẫn động 4 hoặc r, trong khi xe compact và midsize chủ yếu là f.

4.1.2 So sánh khả năng được trang bị dẫn động hai cầu (4WD) giữa xe SUV và Compact

Tạo bảng tần số mới

# Lọc dữ liệu chỉ bao gồm 2 nhóm 'suv' và 'compact'
data_sub <- data1 %>%
  filter(class %in% c("suv", "compact"))

# Tạo bảng 2x2: Class (suv vs compact) ~ DRV (4WD vs Not 4WD)
# Đặt 'suv' là nhóm phơi nhiễm, '4' là sự kiện
table_suv_compact <- table(
  Class = factor(data_sub$class, levels = c("suv", "compact")),
  Is_4WD = ifelse(data_sub$drv == "4", "yes", "no")
)

# Đảm bảo thứ tự cột là 'yes', 'no'
table_suv_compact <- table_suv_compact[, c("yes", "no")]
print(table_suv_compact)

##          Is_4WD
## Class     yes no
##   suv      51 11
##   compact  12 35

Tính Relative Risk (RR) và Odds Ratio (OR)

# Tính RR
rr_result <- riskratio(table_suv_compact)
print(rr_result$measure)

##          risk ratio with 95% C.I.
## Class     estimate    lower    upper
##   suv     1.000000       NA       NA
##   compact 4.197292 2.393927 7.359147

# Tính OR
or_result <- oddsratio(table_suv_compact)
print(or_result$measure)

##          odds ratio with 95% C.I.
## Class     estimate    lower    upper
##   suv      1.00000       NA       NA
##   compact 12.97504 5.310361 34.47191

Nhận xét

• Relative Risk (RR) = 4.197: Tỷ lệ một chiếc xe là xe dẫn động hai cầu (4WD) ở nhóm Compact cao gấp 4.197 lần so với nhóm SUV.
• Odds Ratio (OR) = 12.97: Tỷ lệ để một chiếc xe là xe 4WD ở nhóm Compact cao gấp 12.97 lần so với odds ở nhóm SUV.

Cả hai chỉ số đều cho thấy một mối liên hệ cực kỳ mạnh mẽ. Xe Compact có khả năng được trang bị hệ dẫn động hai cầu cao hơn rất nhiều so với xe SUV.

---

4.2 Phân tích mối quan hệ giữa drv (hệ dẫn động) và cyl (số xi-lanh)

Câu hỏi nghiên cứu: “Số xi-lanh của động cơ, một chỉ số quan trọng về sức mạnh, có liên quan đến loại hệ dẫn động được trang bị trên xe hay không?”

4.2.1 Số xi-lanh có ảnh hưởng đến hệ dẫn động không?

Bảng tần số chéo

table_drv_cyl <- table(data1$drv, data1$cyl)
table_drv_cyl

##    
##      4  5  6  8
##   4 23  0 32 48
##   f 58  4 43  1
##   r  0  0  4 21

Vẽ biểu đồ

df_drv_cyl <- as.data.frame(table_drv_cyl)
colnames(df_drv_cyl) <- c("DRV", "Cylinders", "Count")

ggplot(df_drv_cyl, aes(x = as.factor(Cylinders), y = Count, fill = DRV)) +
  geom_bar(stat = "identity", position = "fill") +
  labs(
    title = "Tỷ lệ hệ dẫn động theo số xi-lanh",
    x = "Số xi-lanh",
    y = "Tỷ lệ",
    fill = "Hệ dẫn động"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

Kiểm định chi bình phương

Giả thuyết thống kê

• Giả thuyết không (H₀): Không có mối liên hệ giữa số xi-lanh và hệ dẫn động.
• Giả thuyết đối (H₁): Có mối liên hệ giữa số xi-lanh và hệ dẫn động.

chisq.test(table_drv_cyl)

## 
##  Pearson's Chi-squared test
## 
## data:  table_drv_cyl
## X-squared = 98.136, df = 6, p-value < 2.2e-16

Nhận xét kết quả kiểm định:

Với p-value < 2.2e-16, ta bác bỏ giả thuyết H0. Có một mối liên hệ có ý nghĩa thống kê rất mạnh mẽ giữa số xi-lanh và loại hệ dẫn động.

4.2.2 Phân tích rủi ro và tỷ số chênh

Để lượng hóa mối quan hệ này, chúng ta sẽ so sánh khả năng một chiếc xe được trang bị hệ dẫn động hai cầu (4WD) giữa hai nhóm động cơ phổ biến và đối lập nhất: 8 xi-lanh (công suất cao) và 4 xi-lanh (công suất thấp).

Tạo bảng tần số 2x2

# Lọc dữ liệu chỉ bao gồm xe 4 và 8 xi-lanh
data_cyl_sub <- data1 %>%
  filter(cyl %in% c(8, 4))

# Tạo bảng 2x2: Số xi-lanh (8 vs 4) ~ Hệ dẫn động (4WD vs Not 4WD)
# Đặt 8 xi-lanh là nhóm phơi nhiễm (exposure group)
table_cyl_4wd <- table(
  Cylinders = factor(data_cyl_sub$cyl, levels = c(8, 4)),
  Is_4WD = ifelse(data_cyl_sub$drv == "4", "yes", "no")
)

# Đảm bảo thứ tự cột "yes", "no"
table_cyl_4wd <- table_cyl_4wd[, c("yes", "no")]
print(table_cyl_4wd)

##          Is_4WD
## Cylinders yes no
##         8  48 22
##         4  23 58

Tính Relative Risk (RR)

# Tính RR
rr_cyl_4wd <- riskratio(table_cyl_4wd)
print(rr_cyl_4wd$measure)

##          risk ratio with 95% C.I.
## Cylinders estimate    lower    upper
##         8 1.000000       NA       NA
##         4 2.278339 1.570253 3.305728

Tính Odds Ratio (OR)

# Tính OR
or_cyl_4wd <- oddsratio(table_cyl_4wd)
print(or_cyl_4wd$measure)

##          odds ratio with 95% C.I.
## Cylinders estimate    lower    upper
##         8 1.000000       NA       NA
##         4 5.410297 2.721399 11.12334

Nhận xét

• Relative Risk (RR) = 2.27: Tỷ lệ một chiếc xe được trang bị hệ dẫn động hai cầu (4WD) ở nhóm 4 xi-lanh cao gấp 2.27 lần so với nhóm 8 xi-lanh. Đây là một sự khác biệt rất lớn.

• Odds Ratio (OR) = 5.41: Tỷ lệ để một chiếc xe là xe 4WD ở nhóm 4 xi-lanh cao gấp 5.41 lần so với odds ở nhóm 8 xi-lanh.

Kết luận:

Cả hai chỉ số RR và OR đều cho thấy một mối liên hệ cực kỳ mạnh mẽ. Việc một chiếc xe được trang bị động cơ 4 xi-lanh làm tăng đáng kể khả năng xe đó cũng được trang bị hệ dẫn động hai cầu.

---

5 Phương pháp Ước lượng Hợp lý Tối đa (Maximum Likelihood Estimation – MLE)

5.1 Giới thiệu chung

Trong thống kê và học máy, dữ liệu thường được mô hình hóa thông qua các phân phối xác suất. Một phân phối xác suất được đặc trưng bởi một số tham số, ví dụ như:

• Phân phối chuẩn (Gaussian): \(\mu\) (trung bình), \(\sigma^2\) (phương sai).
• Phân phối Poisson: tham số đặc trưng là \(\lambda\).

Nếu ta đã biết được dạng của phân phối xác suất, vấn đề đặt ra là: làm sao tìm được các tham số mô tả phân phối đó một cách phù hợp nhất với dữ liệu quan sát? Đó chính là nhiệm vụ của phương pháp ước lượng hợp lý tối đa (Maximum Likelihood Estimation – MLE).

MLE là một công cụ căn bản và quan trọng trong thống kê suy diễn và học máy, giúp chúng ta tìm ra các tham số mô hình sao cho phân phối xác suất sinh ra dữ liệu là “phù hợp nhất” với thực tế đã quan sát.

5.2 Mục tiêu và nguyên lý của MLE

5.2.1 Khái niệm hàm hợp lý (Likelihood Function)

Hàm hợp lý (likelihood function) đo lường mức độ phù hợp của mô hình (tham số) với dữ liệu đã quan sát. Tức là, với một mô hình xác suất có tham số \(\mathbf{w}\), hàm hợp lý cho biết xác suất để toàn bộ tập dữ liệu \(\mathcal{D}\) xảy ra, dưới giả định rằng dữ liệu được sinh ra từ mô hình đó.

Giả sử có một bộ dữ liệu gồm \(N\) quan sát đầu vào:

\[ \mathcal{D} = \{ \mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_N \} \]

và mô hình phân phối xác suất \(P(\mathbf{x}|\mathbf{w})\) phụ thuộc vào tham số

\[ \mathbf{w} = (w_1, w_2, \ldots, w_k)^T \]

Khi đó, hàm hợp lý được định nghĩa là:

\[ L(\mathbf{w}) = P(\mathcal{D}|\mathbf{w}) = P(\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_N|\mathbf{w}) \]

Nếu các quan sát \(\mathbf{x}_i\) là độc lập và phân phối giống nhau (i.i.d – independent and identically distributed), thì hàm hợp lý có thể viết lại thành:

\[ L(\mathbf{w}) = \prod_{i=1}^{N} P(\mathbf{x}_i|\mathbf{w}) \]

5.2.2 Nguyên lý tối đa hóa hợp lý

Nguyên lý cốt lõi của MLE là chọn giá trị tham số \(\hat{\mathbf{w}}\) sao cho hàm hợp lý đạt giá trị lớn nhất. Nói cách khác, ta tìm tham số mà dưới nó thì xác suất tạo ra dữ liệu quan sát là cao nhất:

\[ \hat{\mathbf{w}} = \arg\max_{\mathbf{w}} L(\mathbf{w}) \]

Tuy nhiên, việc tối ưu một tích lớn các xác suất như vậy thường không thuận tiện. Vì vậy, ta thường lấy log của hàm hợp lý để thu được hàm log-hợp lý (log-likelihood function):

\[ \ell(\mathbf{w}) = \log L(\mathbf{w}) = \sum_{i=1}^{N} \log P(\mathbf{x}_i|\mathbf{w}) \]

Do hàm log là đơn điệu tăng, nên điểm cực đại của \(\ell(\mathbf{w})\) cũng chính là cực đại của \(L(\mathbf{w})\). Vậy nên:

\[ \hat{\mathbf{w}} = \arg\max_{\mathbf{w}} \ell(\mathbf{w}) \]

5.3 Tóm tắt quy trình MLE

Bước	Mô tả
Bước 1	Xác định mô hình xác suất \(P(x|\mathbf{w})\) và giả định dữ liệu độc lập (i.i.d).
Bước 2	Viết hàm hợp lý \(L(\mathbf{w})\) hoặc log-hợp lý \(\ell(\mathbf{w})\).
Bước 3	Tính đạo hàm của \(\ell(\mathbf{w})\) theo \(\mathbf{w}\), rồi giải phương trình đạo hàm bằng 0.
Bước 4	Kiểm tra điều kiện tối ưu (đạo hàm bậc hai, tính lồi, v.v.).
Bước 5	Kết luận \(\hat{\mathbf{w}}\) là ước lượng MLE.

5.4 Mô hình hóa bài toán dưới dạng toán học

Trong trường hợp dữ liệu phân phối i.i.d, bài toán MLE trở thành:

\[ \hat{\mathbf{w}} = \arg\max_{\mathbf{w}} \sum_{i=1}^{N} \log P(\mathbf{x}_i|\mathbf{w}) \]

Đây là một bài toán tối ưu hóa có thể giải bằng:

• Đạo hàm tường minh (nếu có thể),
• Gradient Descent / Newton-Raphson (trong trường hợp phức tạp),
• Kỹ thuật tối ưu hóa số.

5.5 Tính chất của ước lượng MLE

Tính chất	Mô tả
Nhất quán (Consistency)	Khi \(N \to \infty\), MLE hội tụ về đúng tham số thực sự.
Bất biến (Invariance)	Nếu \(\hat{\theta}\) là MLE của \(\theta\), thì \(g(\hat{\theta})\) là MLE của \(g(\theta)\).
Phân phối tiệm cận chuẩn	MLE tiệm cận theo phân phối chuẩn khi \(N \to \infty\).
Không thiên lệch tiệm cận	Khi mẫu lớn, sai số trung bình của MLE tiến về 0.

5.6 Kết luận

Phương pháp ước lượng hợp lý tối đa (MLE) là một công cụ then chốt trong thống kê và học máy. Bằng cách tối đa hóa xác suất sinh ra dữ liệu đã quan sát, MLE hướng tới mục tiêu chọn mô hình phù hợp nhất với thực tế. Phương pháp này có cơ sở lý thuyết vững chắc và có thể áp dụng để ước lượng tham số của nhiều phân phối phổ biến như Bernoulli, Gaussian, Poisson, v.v.

---

6 Mô hình Hồi quy Logistic và Probit

6.1 Giới thiệu chung

Khi các phương pháp phân tích bảng chéo và kiểm định Chi bình phương đã xác nhận sự tồn tại của mối quan hệ giữa các biến định tính, bước tiếp theo là xây dựng một mô hình dự báo để lượng hóa tác động của các biến độc lập lên một biến phụ thuộc nhị phân. Khi biến phụ thuộc có hai giá trị (ví dụ: 4WD / Không 4WD), mô hình hồi quy tuyến tính thông thường (OLS) sẽ không còn phù hợp vì nó có thể dự đoán các giá trị nằm ngoài khoảng [0, 1], điều vô lý đối với một xác suất.

Để giải quyết vấn đề này, các mô hình xác suất phi tuyến như Hồi quy Logistic và Hồi quy Probit được sử dụng. Cả hai đều thuộc nhóm Mô hình Tuyến tính Tổng quát (Generalized Linear Models - GLM), với mục tiêu mô hình hóa xác suất xảy ra của một sự kiện, \(P(Y=1|X)\), dựa trên các biến độc lập \(X\).

6.2 Hồi quy Logistic (Logistic Regression)

6.2.1 Mục tiêu và nguyên lý

Hồi quy Logistic không mô hình hóa trực tiếp xác suất \(p\), mà mô hình hóa logarit của tỷ số chênh (log-odds), còn gọi là hàm logit.

• Tỷ số chênh (Odds): Là tỷ lệ giữa xác suất xảy ra sự kiện và xác suất không xảy ra sự kiện: \[ \text{Odds} = \frac{p}{1-p} \]
• Hàm Logit: Là logarit tự nhiên của Odds: \[ \text{logit}(p) = \ln\left(\frac{p}{1-p}\right) \]

Hàm logit có một đặc tính tuyệt vời: nó biến một giá trị xác suất \(p\) trong khoảng (0, 1) thành một giá trị trên toàn bộ trục số thực \((-\infty, +\infty)\). Điều này cho phép chúng ta mô hình hóa nó như một tổ hợp tuyến tính của các biến độc lập:

\[ \ln\left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_k X_{ik} \]

6.2.2 Diễn giải hệ số

Việc diễn giải hệ số \(\beta\) trực tiếp khá phức tạp (“sự thay đổi trong log-odds”). Thay vào đó, chúng ta thường lấy lũy thừa cơ số \(e\) của hệ số để có được Tỷ số chênh (Odds Ratio - OR):

\[ \text{OR} = e^\beta \]

• Nếu \(\beta_1 > 0 \implies e^{\beta_1} > 1\): Khi \(X_1\) tăng 1 đơn vị, odds để sự kiện \(Y=1\) xảy ra sẽ tăng lên \(e^{\beta_1}\) lần, giữ các biến khác không đổi.
• Nếu \(\beta_1 < 0 \implies e^{\beta_1} < 1\): Khi \(X_1\) tăng 1 đơn vị, odds để sự kiện \(Y=1\) xảy ra sẽ giảm xuống còn \(e^{\beta_1}\) lần, giữ các biến khác không đổi.
• Nếu \(\beta_1 = 0 \implies e^{\beta_1} = 1\): Biến \(X_1\) không ảnh hưởng đến odds của sự kiện.

6.3 Hồi quy Probit (Probit Regression)

6.3.1 Mục tiêu và nguyên lý

Hồi quy Probit cũng mô hình hóa xác suất \(p\), nhưng thay vì dùng hàm logit, nó sử dụng hàm phân phối tích lũy (CDF) của phân phối chuẩn tắc, ký hiệu là \(\Phi(z)\).

Ý tưởng cơ bản là có một “biến tiềm ẩn” (latent variable) \(Y^*\) tuân theo phân phối chuẩn. Chúng ta không quan sát được \(Y^*\), mà chỉ quan sát được \(Y=1\) nếu \(Y^*\) vượt qua một ngưỡng nào đó.

Mô hình Probit có dạng: \[ P(Y=1|X) = \Phi(\beta_0 + \beta_1 X_1 + \beta_2 X_2 + \ldots + \beta_k X_k) \] Hoặc viết theo dạng hàm ngược (hàm probit): \[ \Phi^{-1}(p_i) = \beta_0 + \beta_1 X_{i1} + \beta_2 X_{i2} + \ldots + \beta_k X_{ik} \]

6.3.2 Diễn giải hệ số

Hệ số \(\beta\) trong mô hình Probit thể hiện sự thay đổi trong z-score (giá trị trên trục hoành của phân phối chuẩn tắc) khi biến độc lập \(X\) thay đổi 1 đơn vị. Việc diễn giải này không trực quan và thường ít được sử dụng trong báo cáo kinh doanh. Thay vào đó, người ta thường tính toán và báo cáo tác động biên (marginal effects), tức là sự thay đổi trong xác suất \(p\) khi \(X\) thay đổi.

6.4 Logistic vs. Probit: Lựa chọn nào?

Tiêu chí	Hồi quy Logistic	Hồi quy Probit
Hàm liên kết	Logit, dựa trên phân phối Logistic	Probit, dựa trên phân phối Chuẩn
Độ dốc	Đường cong S (sigmoid) hơi phẳng hơn ở đuôi	Đường cong S (sigmoid) dốc hơn ở đuôi
Diễn giải	Hệ số có thể diễn giải trực tiếp qua Odds Ratio (OR)	Hệ số khó diễn giải, thường dùng tác động biên
Sử dụng	Rất phổ biến trong y sinh, khoa học xã hội, marketing	Phổ biến trong kinh tế lượng

Trong thực tế, cả hai mô hình thường cho kết quả dự đoán rất giống nhau. Hồi quy Logistic thường được ưa chuộng hơn vì tính dễ diễn giải của Odds Ratio.

6.5 Áp dụng mô hình dự đoán hệ dẫn động hai cầu (4WD)

Câu hỏi nghiên cứu: Những yếu tố kỹ thuật nào (dung tích động cơ, số xi-lanh) ảnh hưởng đến khả năng một chiếc xe được trang bị hệ dẫn động hai cầu (4WD)?

6.5.1 Chuẩn bị dữ liệu

Chúng ta sẽ tạo một biến phụ thuộc mới là is_4wd, nhận giá trị 1 nếu xe là 4WD và 0 trong trường hợp còn lại.

# Tạo biến phụ thuộc nhị phân
mpg_model_data <- mpg %>%
  mutate(is_4wd = ifelse(drv == "4", 1, 0))

# Xem qua dữ liệu mới
head(mpg_model_data %>% select(manufacturer, model, drv, is_4wd, displ, cyl))

6.5.2 Xây dựng và phân tích mô hình Hồi quy Logistic

Chúng ta sẽ xây dựng mô hình dự đoán is_4wd dựa trên displ (dung tích động cơ) và cyl (số xi-lanh).

# Xây dựng mô hình Logistic
logit_model <- glm(is_4wd ~ displ + cyl, 
                   data = mpg_model_data, 
                   family = binomial(link = "logit"))

# Xem kết quả
summary(logit_model)

## 
## Call:
## glm(formula = is_4wd ~ displ + cyl, family = binomial(link = "logit"), 
##     data = mpg_model_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -2.2004     0.6051  -3.637 0.000276 ***
## displ         0.7385     0.3118   2.368 0.017871 *  
## cyl          -0.1056     0.2412  -0.438 0.661626    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 321.03  on 233  degrees of freedom
## Residual deviance: 289.07  on 231  degrees of freedom
## AIC: 295.07
## 
## Number of Fisher Scoring iterations: 4

Diễn giải kết quả mô hình Logistic:

• Hệ số của displ (Dung tích động cơ):
  • Hệ số là 0.7385. Giá trị này dương và có ý nghĩa thống kê (p-value = 0.017 < 0.05).
  • Để diễn giải, ta tính Odds Ratio: exp(0.7385) ≈ 2.09.
  • Diễn giải: Nếu số xi-lanh không đổi thì khi dung tích động cơ tăng 1 lít, tỷ lệ (odds) để chiếc xe đó là xe 4WD tăng lên 2.09 lần.
• Hệ số của cyl (Số xi-lanh):
  • Hệ số là -0.1056. Giá trị âm nhưng không có ý nghĩa thống kê ở mức 5% (p-value = 0.661 > 0.05).

6.5.3 Xây dựng và phân tích mô hình Hồi quy Probit

# Xây dựng mô hình Probit
probit_model <- glm(is_4wd ~ displ + cyl, 
                    data = mpg_model_data, 
                    family = binomial(link = "probit"))

# Xem kết quả
summary(probit_model)

## 
## Call:
## glm(formula = is_4wd ~ displ + cyl, family = binomial(link = "probit"), 
##     data = mpg_model_data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.37030    0.36683  -3.736 0.000187 ***
## displ        0.42691    0.18580   2.298 0.021581 *  
## cyl         -0.04788    0.14594  -0.328 0.742878    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 321.03  on 233  degrees of freedom
## Residual deviance: 289.29  on 231  degrees of freedom
## AIC: 295.29
## 
## Number of Fisher Scoring iterations: 5

Diễn giải kết quả mô hình Probit:

• Hệ số của displ (Dung tích động cơ):
  • Hệ số là 0.4269. Giá trị này dương và có ý nghĩa thống kê (p-value = 0.0001 < 0.05).
  • Diễn giải: Nếu số xi-lanh không đổi thì khi dung tích động cơ tăng 1 lít, z-score dự đoán xác suất xe là 4WD tăng lên 0.4269 đơn vị.
• Hệ số của cyl (Số xi-lanh):
  • Hệ số là -0.0478. Giá trị âm nhưng không có ý nghĩa thống kê ở mức 5% (p-value = 0.7428 > 0.05).

6.5.4 So sánh hai mô hình

Cả hai mô hình đều cho cùng một kết luận: dung tích động cơ (displ) là yếu tố dự báo quan trọng cho việc một chiếc xe có phải là 4WD hay không. Chúng ta có thể so sánh trực tiếp xác suất dự đoán từ hai mô hình.

# Lấy xác suất dự đoán từ hai mô hình
mpg_model_data$pred_logit <- predict(logit_model, type = "response")
mpg_model_data$pred_probit <- predict(probit_model, type = "response")

# Vẽ biểu đồ so sánh
ggplot(mpg_model_data, aes(x = pred_logit, y = pred_probit)) +
  geom_point(alpha = 0.5) +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") +
  labs(
    title = "So sánh xác suất dự đoán giữa mô hình Logit và Probit",
    x = "Xác suất dự đoán từ mô hình Logistic",
    y = "Xác suất dự đoán từ mô hình Probit"
  ) +
  theme_minimal()

Kết luận:

Biểu đồ cho thấy các điểm dữ liệu nằm rất sát đường chéo màu đỏ, chứng tỏ xác suất dự đoán từ hai mô hình gần như tương đương nhau.