library(readr)

## Warning: package 'readr' was built under R version 4.3.3

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.3.3

library(dplyr)

## Warning: package 'dplyr' was built under R version 4.3.3

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(skimr)

## Warning: package 'skimr' was built under R version 4.3.3

library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

PHẦN 1. GIỚI THIỆU VỀ BỘ DỮ LIỆU

Bộ dữ liệu hiện tại bao gồm thông tin chi tiết từ 3.000 hồ sơ bệnh nhân nhiễm COVID-19 tại Đan Mạch, được trình bày dưới dạng dữ liệu bảng, trong đó mỗi dòng đại diện cho một cá nhân riêng biệt. Tổng cộng có 20 biến định danh và mô tả đặc trưng, phản ánh đa chiều các yếu tố dịch tễ học, đặc điểm nền, quá trình điều trị và kết cục lâm sàng.

Cụ thể, dữ liệu ghi nhận các thông tin nền tảng như giới tính, vùng cư trú, nghề nghiệp, tiền sử bệnh nền, tình trạng hút thuốc và chỉ số BMI. Về tiến trình bệnh, các biến thể hiện thời điểm phát hiện nhiễm bệnh, triệu chứng khởi phát, mức độ nghiêm trọng, và các can thiệp điều trị bao gồm nhập viện, điều trị ICU, hỗ trợ thở máy và ngày xuất viện. Các biến kết quả như hồi phục, tái nhiễm và hội chứng hậu COVID cũng được theo dõi đầy đủ. Cuối cùng, bộ dữ liệu còn bao gồm thông tin liên quan đến miễn dịch học như tình trạng tiêm chủng, loại vắc-xin đã sử dụng và khoảng thời gian kể từ mũi tiêm gần nhất.

Với cấu trúc đầy đủ gồm 20 biến, bộ dữ liệu này cho phép thực hiện các phân tích thống kê từ mô tả đến suy luận, đồng thời tạo điều kiện cho việc kiểm định mối liên hệ giữa các yếu tố nguy cơ và kết quả điều trị trong bối cảnh thực tiễn tại Đan Mạch.

d <- read.csv(file.choose())

str(d)

## 'data.frame':    3000 obs. of  20 variables:
##  $ Patient_ID             : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Gioi_tinh              : chr  "Male" "Male" "Female" "Female" ...
##  $ Khu_vuc                : chr  "Hovedstaden" "Sjælland" "Syddanmark" "Hovedstaden" ...
##  $ Preexisting_Condition  : chr  "Obesity" "Asthma" "Hypertension" "Asthma" ...
##  $ Date_of_Infection      : chr  "6/21/2022" "2/2/2024" "5/28/2023" "8/13/2023" ...
##  $ COVID_Strain           : chr  "Delta" "XBB.1.5" "Beta" "Delta" ...
##  $ Trieu_chung            : chr  "Mild" "Mild" "Mild" "Severe" ...
##  $ Severity               : chr  "Moderate" "Moderate" "High" "High" ...
##  $ Nhap_vien              : chr  "Yes" "No" "Yes" "No" ...
##  $ Hospital_Admission_Date: chr  "1/13/2025" "" "3/7/2025" "" ...
##  $ Hospital_Discharge_Date: chr  "1/26/2025" "" "4/26/2025" "" ...
##  $ ICU_Admission          : chr  "No" "No" "Yes" "No" ...
##  $ Ventilator_Support     : chr  "No" "No" "Yes" "No" ...
##  $ Hoi_phuc               : chr  "Yes" "No" "No" "Yes" ...
##  $ Date_of_Recovery       : chr  "4/19/2023" "" "" "2/9/2025" ...
##  $ Tai_nhiem              : chr  "No" "No" "No" "Yes" ...
##  $ Date_of_Reinfection    : chr  "" "" "" "8/24/2024" ...
##  $ Tiem_chung             : chr  "Yes" "No" "Yes" "Yes" ...
##  $ Vaccine_Type           : chr  "None" "None" "Janssen" "AstraZeneca" ...
##  $ Nghe_nghiep            : chr  "Healthcare" "Healthcare" "Unemployed" "Office Worker" ...

Các biến trong bộ dữ liệu bao gồm:

Tên biến	Mô tả
Patient_ID	Mã ẩn danh duy nhất của mỗi bệnh nhân
Gioi_tinh	Giới tính của bệnh nhân
Khu_vuc	Khu vực cư trú của bệnh nhân
Nghe_nghiep	Nghề nghiệp hoặc lĩnh vực công tác của bệnh nhân
COVID_Strain	Biến chủng virus
Trieu_chung	Danh sách triệu chứng chính khi nhiễm
Severity	Mức độ nặng của bệnh
Nhap_vien	Cho biết bệnh nhân có nhập viện hay không
Hospital_Admission_Date	Ngày bệnh nhân nhập viện
Hospital_Discharge_Date	Ngày bệnh nhân xuất viện
ICU_Admission	Cho biết bệnh nhân có từng vào đơn vị ICU hay không
Ventilator_Support	Cho biết bệnh nhân có cần hỗ trợ thở máy hay không
Hoi_phuc	Cho biết bệnh nhân đã hồi phục hoàn toàn hay chưa
Date_of_Recovery	Ngày xác nhận bệnh nhân hồi phục
Tai_nhiem	Cho biết bệnh nhân có bị tái nhiễm hay không
Date_of_Tai_nhiem	Ngày xác nhận tái nhiễm SARS-CoV-2
Tiem_chung	Tình trạng tiêm chủng COVID-19
Vaccine_Type	Loại vaccine đã tiêm

PHẦN 2. PHÂN TÍCH MỘT BIẾN ĐỊNH TÍNH

Bộ dữ liệu bao gồm 9 biến định tính

dldt <- c("Gioi_tinh","Khu_vuc","Nghe_nghiep", "Trieu_chung","Nhap_vien","Hoi_phuc","Tai_nhiem","Tiem_chung","Vaccine_Type")
dldt

## [1] "Gioi_tinh"    "Khu_vuc"      "Nghe_nghiep"  "Trieu_chung"  "Nhap_vien"   
## [6] "Hoi_phuc"     "Tai_nhiem"    "Tiem_chung"   "Vaccine_Type"

dt <- d[, dldt]

dt <- data.frame(lapply(dt, as.factor))

Biến Nhap_vien

Bảng tần số và tần suất

# 1. Bảng tần số của biến Nhap_vien
tanso_hosp <- table(dt$Nhap_vien)
print(tanso_hosp)

## 
##   No  Yes 
## 2124  876

# 2. Bảng tần suất (tỷ lệ) của biến Nhap_vien
tansuat_hosp <- prop.table(tanso_hosp)
print(round(tansuat_hosp, 4)*100)            # làm tròn 4 chữ số thập phân

## 
##   No  Yes 
## 70.8 29.2

Biểu đồ

# 1. Tính bảng tần số và tần suất
tanso_hosp   <- table(dt$Nhap_vien)
tansuat_hosp <- prop.table(tanso_hosp)

# 2. Màu cho từng lát
cols <- c("pink", "green")

# 3. Nhãn phần trăm cho legend
nhalabels <- paste0(round(tansuat_hosp * 100, 1), "%")

# 4. Vẽ pie chart không nhãn ngay trên lát
pie(
  tansuat_hosp,
  labels = NA,               
  main   = "Tần suất bệnh nhân nhập viện",
  col    = cols,
  border = "white"
)

# 5. Thêm legend bên ngoài góc phải
legend(
  x      = "topright",
  legend = paste(names(tansuat_hosp), nhalabels),
  fill   = cols
)

Nhận xét

Trong số 3.000 bệnh nhân, có 2.124 người (chiếm 70,8%) không phải nhập viện và 876 người (29,2%) phải nhập viện. Tỷ lệ nhập viện gần 30% cho thấy khoảng 3 trên 10 ca nhiễm COVID-19 trong mẫu này cần được điều trị tại bệnh viện.

Biến Nghe_nghiep

Bảng tần số và tần suất

# 1. Bảng tần số của biến Nghe_nghiep
tanso_occ <- table(dt$Nghe_nghiep)
print(tanso_occ)

## 
##        Driver    Healthcare Office Worker       Student       Teacher 
##           492           484           503           533           491 
##    Unemployed 
##           497

# 2. Bảng tần suất (tỷ lệ) của biến Nghe_nghiep
tansuat_occ <- prop.table(tanso_occ)
print(round(tansuat_occ, 4)*100)           # tỷ lệ làm tròn 4 chữ số thập phân

## 
##        Driver    Healthcare Office Worker       Student       Teacher 
##         16.40         16.13         16.77         17.77         16.37 
##    Unemployed 
##         16.57

Biểu đồ

library(ggplot2)

ggplot(dt, aes(x = Nghe_nghiep, fill = Nghe_nghiep)) +
  geom_bar(width = 0.7) +
  labs(title = "Tần số bệnh nhân theo Nghe_nghiep",
       x     = "Nghề nghiệp",
       y     = "Số bệnh nhân") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

Nhận xét

Trong số 3.000 bệnh nhân, nhóm “Student” chiếm nhiều nhất với 533 người, tương đương 17,77%. Theo sau là “Office Worker” với 503 ca (16,77%) và “Unemployed” với 497 ca (16,57%). Ba nhóm còn lại gồm “Driver” có 492 ca (16,40%), “Teacher” 491 ca (16,37%) và “Healthcare Worker” 484 ca (16,13%). Khoảng cách giữa nhóm cao nhất và thấp nhất chỉ là 1,64 phần trăm, cho thấy các nghề nghiệp trong mẫu được phân bố tương đối đồng đều.

Biến Trieu_chung

Bảng tần số và tần suất

# 1. Bảng tần số của biến Trieu_chung
tanso_Trieu_chung <- table(dt$Trieu_chung)
print(tanso_Trieu_chung)

## 
##     Mild Moderate   Severe 
##     1010      978     1012

# 2. Bảng tần suất (tỷ lệ) của biến Trieu_chung
tansuat_Trieu_chung <- prop.table(tanso_Trieu_chung)
print(round(tansuat_Trieu_chung, 4)*100)       # làm tròn 4 chữ số thập phân

## 
##     Mild Moderate   Severe 
##    33.67    32.60    33.73

Biểu đồ

library(ggplot2)

ggplot(dt, aes(x = Trieu_chung, fill = Trieu_chung)) +
  geom_bar(width = 0.7) +
  labs(title = "Tần số bệnh nhân theo Trieu_chung",
       x     = "Triệu chứng",
       y     = "Số bệnh nhân") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.position = "none"
  )

Nhận xét

Trong tổng số 3.000 bệnh nhân, có 1.010 người (33,67%) thuộc nhóm triệu chứng nhẹ, 978 người (32,60%) ở nhóm trung bình và 1.012 người (33,73%) ở nhóm nặng. Chênh lệch giữa nhóm nhiều nhất và ít nhất chỉ khoảng 1,13 phần trăm, cho thấy mức độ biểu hiện triệu chứng được phân phối khá đồng đều.

Biến Hoi_phuc

Bảng tần số và tần suất

# 1. Bảng tần số
tanso_rec <- table(dt$Hoi_phuc)
print(tanso_rec)

## 
##   No  Yes 
## 1492 1508

# 2. Bảng tần suất (tỷ lệ)
tansuat_rec <- prop.table(tanso_rec)
print(round(tansuat_rec, 4)*100)   # làm tròn 4 chữ số thập phân

## 
##    No   Yes 
## 49.73 50.27

Biểu đồ

# Màu cho từng lát (Yes / No)
cols <- c("lightblue", "orange")

# Nhãn phần trăm cho legend
nhalabels <- paste0(round(tansuat_rec * 100, 1), "%")

# Vẽ pie chart
pie(
  tansuat_rec,
  labels = NA,
  main   = "Tỷ lệ bệnh nhân hồi phục (Hoi_phuc)",
  col    = cols,
  border = "white"
)

# Thêm legend bên ngoài góc phải
legend(
  x      = "topright",
  legend = paste(names(tansuat_rec), nhalabels),
  fill   = cols,
  bty    = "n"
)

Nhận xét

Trong số 3.000 bệnh nhân, 1.508 người (50,27%) đã phục hồi và 1.492 người (49,73%) chưa phục hồi. Tỷ lệ hồi phục nhỉnh hơn chỉ 0,54 phần trăm, cho thấy sức khỏe của những người từng nhiễm COVID-19 trong mẫu này được cải thiện ở mức gần như cân bằng.

Biến Tai_nhiem

Bảng tần số và tần suất

# 1. Bảng tần số
tanso_reinf <- table(dt$Tai_nhiem)
print(tanso_reinf)

## 
##   No  Yes 
## 2715  285

# 2. Bảng tần suất (tỷ lệ)
tansuat_reinf <- prop.table(tanso_reinf)
print(round(tansuat_reinf, 4)*100)   # làm tròn 4 chữ số thập phân

## 
##   No  Yes 
## 90.5  9.5

Biểu đồ

# Màu cho từng lát (Yes / No)
cols <- c("lightblue","pink")

# Nhãn phần trăm cho legend
nhalabels <- paste0(round(tansuat_reinf * 100, 1), "%")

# Vẽ pie chart không nhãn ngay trên lát
pie(
  tansuat_reinf,
  labels = NA,
  main   = "Tỷ lệ tái nhiễm (Tai_nhiem)",
  col    = cols,
  border = "white"
)

# Thêm legend bên ngoài góc phải
legend(
  x      = "topright",
  legend = paste(names(tansuat_reinf), nhalabels),
  fill   = cols,
  bty    = "n"
)

Nhận xét

Trong 3.000 bệnh nhân, 285 người (9,5%) tái nhiễm COVID-19, còn lại 2.715 người (90,5%) không tái nhiễm. Điều này cho thấy tỷ lệ tái nhiễm trong mẫu vào khoảng 1 trên 10 ca.

Biến Gioi_tinh

Bảng tần số và tần suất

# Bảng tần số của biến Gioi_tinh
tanso_Gioi_tinh <- table(dt$Gioi_tinh)
print(tanso_Gioi_tinh)

## 
## Female   Male 
##   1527   1473

# Bảng tần suất (prop.table tự chia cho tổng số quan sát)
tansuat_Gioi_tinh <- prop.table(tanso_Gioi_tinh)
print(round(tansuat_Gioi_tinh, 4)*100)  # làm tròn 4 chữ số thập phân

## 
## Female   Male 
##   50.9   49.1

Biểu đồ

# Tính bảng tần số
tanso_Gioi_tinh <- table(dt$Gioi_tinh)

# Vẽ barplot tần số
barplot(tanso_Gioi_tinh,
        main = "Phân phối số lượng theo Gioi_tinh",
        xlab = "Giới tính",
        ylab = "Số lượng bệnh nhân",
        col  = c("steelblue", "salmon"),
        border = NA)

Nhận xét

Trong tổng số 3.000 bệnh nhân, có 1.527 người là nữ (chiếm 50,9%) và 1.473 người là nam (chiếm 49,1%). Chênh lệch giữa hai nhóm chỉ 54 người, tương ứng 1,8% về tỷ lệ, cho thấy số lượng nam và nữ trong mẫu gần như ngang nhau.

Biến Tiem_chung

Bảng tần số và tần suất

# 1. Bảng tần số
tanso_vacc <- table(dt$Tiem_chung)
print(tanso_vacc)

## 
##   No  Yes 
## 1528 1472

# 2. Bảng tần suất (tỷ lệ)
tansuat_vacc <- prop.table(tanso_vacc)
print(round(tansuat_vacc, 4)*100)   # làm tròn 4 chữ số thập phân

## 
##    No   Yes 
## 50.93 49.07

Biểu đồ

# Màu cho từng lát (Vaccinated / Unvaccinated)
cols <- c("lightpink", "lightblue")

# Nhãn phần trăm cho legend
nhalabels <- paste0(round(tansuat_vacc * 100, 1), "%")

# Vẽ pie chart không nhãn ngay trên lát
pie(
  tansuat_vacc,
  labels = NA,
  main   = "Tỷ lệ tiêm chủng (Tiem_chung)",
  col    = cols,
  border = "white"
)

# Thêm legend bên ngoài góc phải
legend(
  x      = "topright",
  legend = paste(names(tansuat_vacc), nhalabels),
  fill   = cols,
  bty    = "n"
)

Nhận xét

Trong 3.000 bệnh nhân, 1.472 người đã tiêm vắc-xin (49,1%) và 1.528 người chưa tiêm (50,9%). Chênh lệch giữa hai nhóm chỉ 1,8 phần trăm, cho thấy tỷ lệ tiêm chủng trong mẫu gần như cân bằng.

Biến Khu_vuc

Bảng tần số và tần suất

# 1. Bảng tần số của biến Khu_vuc
tanso_Khu_vuc <- table(dt$Khu_vuc)
print(tanso_Khu_vuc)

## 
## Hovedstaden Midtjylland Nordjylland    Sjælland  Syddanmark 
##         612         628         554         609         597

# 2. Bảng tần suất (tỷ lệ) của biến Khu_vuc
tansuat_Khu_vuc <- prop.table(tanso_Khu_vuc)
print(round(tansuat_Khu_vuc, 4)*100)       # làm tròn 4 chữ số thập phân

## 
## Hovedstaden Midtjylland Nordjylland    Sjælland  Syddanmark 
##       20.40       20.93       18.47       20.30       19.90

Biểu đồ

library(ggplot2)

ggplot(dt, aes(x = Khu_vuc, fill = Khu_vuc)) +
  geom_bar(width = 0.6) +
  labs(title = "Tần số bệnh nhân theo Khu_vuc",
       x     = "Khu vực",
       y     = "Số bệnh nhân") +
  theme_minimal() +
  theme(legend.position = "none")

Nhận xét

Trong 3.000 bệnh nhân, vùng Midtjylland có số ca cao nhất với 628 người (20,93%), tiếp theo là Hovedstaden 612 ca (20,40%) và Sjælland 609 ca (20,30%). Vùng Syddanmark đứng thứ tư với 597 ca (19,90%) và Nordjylland có số ca thấp nhất 554 người (18,47%). Khoảng cách giữa vùng có tỷ lệ cao nhất và thấp nhất chỉ khoảng 2,46 phần trăm, cho thấy sự phân bố bệnh nhân giữa năm khu vực tương đối đồng đều.

Bảng tần số và tần suất

# 1. Bảng tần số của biến Preexisting_Condition
tanso_pre <- table(dt$Preexisting_Condition)
print(tanso_pre)

## < table of extent 0 >

# 2. Bảng tần suất (tỷ lệ) của biến Preexisting_Condition
tansuat_pre <- prop.table(tanso_pre)
print(round(tansuat_pre, 4)*100)            # tỷ lệ làm tròn 4 chữ số thập phân

## numeric(0)

Biến Vaccine_Type

Bảng tần số và tần suất

# 1. Bảng tần số
tanso_vtype <- table(dt$Vaccine_Type)
print(tanso_vtype)

## 
## AstraZeneca     Janssen     Moderna        None      Pfizer 
##         295         293         292        1809         311

# 2. Bảng tần suất (tỷ lệ)
tansuat_vtype <- prop.table(tanso_vtype)
print(round(tansuat_vtype, 4)*100)   # làm tròn 4 chữ số thập phân

## 
## AstraZeneca     Janssen     Moderna        None      Pfizer 
##        9.83        9.77        9.73       60.30       10.37

Biểu đồ

# 3. Chuẩn bị màu cho từng loại vaccine
cols <- rainbow(length(tanso_vtype))

# 4. Nhãn phần trăm cho legend
nhalabels <- paste0(round(tansuat_vtype * 100, 1), "%")

# 5. Vẽ pie chart không hiển thị nhãn trên lát
pie(
  tansuat_vtype,
  labels = NA,
  main   = "Tỷ lệ các loại vaccine (Vaccine_Type)",
  col    = cols,
  border = "white"
)

# 6. Thêm legend bên ngoài góc phải
legend(
  x      = "topright",
  legend = paste(names(tansuat_vtype), nhalabels),
  fill   = cols,
  cex    = 0.5,
  bty    = "n"
)

Nhận xét

Trong 3.000 bệnh nhân, 1.809 người (60,30%) chưa tiêm loại vaccine nào, chiếm đa số, trong khi các loại vaccine còn lại có số ca gần như nhau: AstraZeneca 295 ca (9,83%), Janssen 293 ca (9,77%), Moderna 292 ca (9,73%) và Pfizer 311 ca (10,37%). Sự chênh lệch giữa các loại vaccine đã tiêm (ngoại trừ “None”) chỉ dao động từ 9,73% đến 10,37%, cho thấy không có loại nào chiếm ưu thế tuyệt đối trong nhóm đã tiêm.

PHẦN 3. ƯỚC LƯỢNG KHOẢNG TIN CẬY

Nhap_vien (Nhập viện = “Yes”)

Giả thuyết kiểm định

H₀: Tỷ lệ bệnh nhân phải nhập viện trong dân số là 50 % (p = 0.5).
H₁: Tỷ lệ bệnh nhân phải nhập viện khác 50 % (p ≠ 0.5).

n_hosp      <- sum(dt$Nhap_vien == "Yes")
n_total     <- nrow(dt)
ci_hosp     <- prop.test(n_hosp, n_total, conf.level = 0.95, correct = FALSE)
print(ci_hosp)

## 
##  1-sample proportions test without continuity correction
## 
## data:  n_hosp out of n_total, null probability 0.5
## X-squared = 519.17, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.2760039 0.3085281
## sample estimates:
##     p 
## 0.292

Nhận xét

Thống kê kiểm định: \(X^2=519.17\) , df=1
p-value < 2.2 × 10⁻¹⁶.
Khoảng tin cậy 95 % cho p: [0.2760, 0.3085] (27.60 % – 30.85 %).

Với mức ý nghĩa α = 0.05 và p-value rất nhỏ (< 2.2 × 10⁻¹⁶), bác bỏ H₀. Điều này có nghĩa tỷ lệ nhập viện thực sự không phải 50 % mà xấp xỉ 29.2 %, và với 95 % độ tin cậy, tỷ lệ này nằm trong khoảng [27.6 %, 30.9 %]. Kết quả cho thấy tỷ lệ nhập viện trong mẫu có sự khác biệt rõ rệt so với giả thiết 50 %.

Tai_nhiem (Tái nhiễm = “Yes”)

Giả thuyết kiểm định

H₀: Tỷ lệ bệnh nhân tái nhiễm trong dân số là 50 % (p = 0.5).
H₁: Tỷ lệ bệnh nhân tái nhiễm khác 50 % (p ≠ 0.5).

n_reinf     <- sum(dt$Tai_nhiem == "Yes")
ci_reinf    <- prop.test(n_reinf, n_total, conf.level = 0.95, correct = FALSE)
print(ci_reinf)

## 
##  1-sample proportions test without continuity correction
## 
## data:  n_reinf out of n_total, null probability 0.5
## X-squared = 1968.3, df = 1, p-value < 2.2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.08501949 0.10601638
## sample estimates:
##     p 
## 0.095

Nhận xét

Thống kê kiểm định: \(X^2=1968.3\) , df=1
p-value < 2.2 × 10⁻¹⁶.
Khoảng tin cậy 95 % cho p: [0.0850, 0.1060] (8.50 % – 10.60 %).

Với mức ý nghĩa α = 0.05 và p-value rất nhỏ (< 2.2 × 10⁻¹⁶, bác bỏ H₀. Điều này chứng tỏ tỷ lệ tái nhiễm thực sự không bằng 50 % mà chỉ khoảng 9.5 %, và với 95 % độ tin cậy, tỷ lệ này nằm trong khoảng [8.50 %, 10.60 %]. Kết quả cho thấy tỷ lệ tái nhiễm trong mẫu khác biệt rất rõ rệt so với giả thiết 50 %.

Tiem_chung (Đã tiêm = “Yes”)

Giả thuyết kiểm định

H₀: Tỷ lệ bệnh nhân đã tiêm trong dân số là 50 % (p = 0.5).
H₁: Tỷ lệ bệnh nhân đã tiêm khác 50 % (p ≠ 0.5).

n_vac       <- sum(dt$Tiem_chung == "Yes")
ci_vac      <- prop.test(n_vac, n_total, conf.level = 0.95, correct = FALSE)
print(ci_vac)

## 
##  1-sample proportions test without continuity correction
## 
## data:  n_vac out of n_total, null probability 0.5
## X-squared = 1.0453, df = 1, p-value = 0.3066
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4728012 0.5085560
## sample estimates:
##         p 
## 0.4906667

Nhận xét

Thống kê kiểm định: \(X^2=1.0453\) , df=1
p-value = 0.3066.
Khoảng tin cậy 95 % cho p: [0.4728, 0.5086] (47.28 % – 50.86 %).

Với mức ý nghĩa α=0.05 và p-value = 0.3066 > 0.05, ta không đủ cơ sở để bác bỏ H₀. Điều này có nghĩa tỷ lệ bệnh nhân đã tiêm vắc-xin trong dân số không khác biệt có ý nghĩa so với 50 %; với 95 % độ tin cậy, tỷ lệ này nằm trong khoảng [47.3 %, 50.9 %].

PHẦN 4. PHÂN TÍCH MỐI QUAN HỆ GIỮA HAI BIẾN ĐỊNH TÍNH

Chọn biến Nhap_vien và Tai_nhiem làm biến phụ thuộc để phân tích các yếu tố ảnh hưởng đến việc bệnh nhân có phải nhập viện và có bị tái nhiễm COVID-19 hay không.

# Đảm bảo R dùng UTF-8
options(encoding = "UTF-8")
Sys.setlocale("LC_CTYPE", "en_US.UTF-8")  # hoặc "Vietnamese_Vietnam.1258" trên Windows

## [1] "en_US.UTF-8"

Nhap_vien và Hoi_phuc

Bảng tần số và tần suất

# Gán nhãn rõ ràng cho hai biến
dt$Hoi_phuc <- factor(
  dt$Hoi_phuc,
  levels = c("No", "Yes"),
  labels = c("Chưa hồi phục", "Hồi phục")
)

# Lập bảng tần số chéo
hr_tab2  <- table(dt$Nhap_vien, dt$Hoi_phuc)
hr_tab2a <- addmargins(hr_tab2)
print(hr_tab2a)

##      
##       Chưa hồi phục Hồi phục  Sum
##   No           1050     1074 2124
##   Yes           442      434  876
##   Sum          1492     1508 3000

# Bảng tần suất theo hàng (trong mỗi nhóm Nhap_vien)
hr_prop2         <- prop.table(hr_tab2, margin = 1)
hr_prop2_rounded <- round(hr_prop2, 4)
print(hr_prop2_rounded)

##      
##       Chưa hồi phục Hồi phục
##   No         0.4944   0.5056
##   Yes        0.5046   0.4954

Đồ thị

# Chuyển bảng tần số sang dạng data frame
hr_df2 <- as.data.frame(hr_tab2)
colnames(hr_df2) <- c("Nhap_vien", "Hoi_phuc", "Count")

# Vẽ biểu đồ cột nhóm
library(ggplot2)

ggplot(hr_df2, aes(x = Hoi_phuc, y = Count, fill = Nhap_vien)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.7), width = 0.6) +
  labs(
    title = "Số bệnh nhân nhập viện theo tình trạng hồi phục",
    x     = "Tình trạng hồi phục",
    y     = "Số bệnh nhân",
    fill  = "Tình trạng nhập viện"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "Chưa nhập viện" = "darkblue",
    "Nhập viện"       = "pink"
  ))

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

Nhận xét

Trong số 3.000 bệnh nhân, có 1.508 người đã hồi phục và 1.492 người vẫn chưa hồi phục. Xét nhóm không nhập viện (2.124 ca), 1.074 người (50,56 %) đã hồi phục và 1.050 người (49,44 %) chưa hồi phục. Ngược lại, trong 876 ca phải nhập viện, chỉ có 434 người (49,54 %) hồi phục và 442 người (50,46 %) chưa hồi phục. Như vậy, dù nhóm không nhập viện có số ca hồi phục cao hơn (1.074 so với 434), tỷ lệ hồi phục giữa hai nhóm cũng rất gần nhau, chỉ chênh lệch khoảng 1 phần trăm, cho thấy việc nhập viện không dẫn đến khác biệt lớn về khả năng hồi phục trong mẫu này.

Kiểm định Thống kê

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

H₀:Tình trạng nhập viện và khả năng hồi phục là hai biến độc lập.
H₁:Có mối liên hệ giữa nhập viện và hồi phục.

chi_hr2 <- chisq.test(hr_tab2)
print(chi_hr2)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  hr_tab2
## X-squared = 0.21967, df = 1, p-value = 0.6393

Nhận xét

Giá trị thống kê Chi bình phương: X² = 0.21967
Bậc tự do (df) = 1
p-value = 0.6393

Với α = 0.05, p-value (0.6393) > 0.05, không đủ cơ sở để bác bỏ H₀.Vậy tình trạng nhập viện và khả năng hồi phục là hai biến độc lập.

Nói cách khác, trong mẫu 3.000 bệnh nhân này, tỷ lệ hồi phục ở nhóm không nhập viện (1.074/2.124 ≈ 50,6 %) và nhóm nhập viện (434/876 ≈ 49,5 %) không khác biệt một cách có ý nghĩa thống kê. Sự chênh lệch chỉ khoảng 1,1 % cho thấy việc phải nhập viện không ảnh hưởng rõ rệt đến khả năng hồi phục của bệnh nhân.

Hiệu tỷ lệ (Risk Difference - RD)

Giả thuyết kiểm định

H₀: p₁ − p₂ = 0 (Tỷ lệ hồi phục ở nhóm chưa nhập viện bằng tỷ lệ hồi phục ở nhóm nhập viện)
H₁: p₁ − p₂ ≠ 0 (Có sự khác biệt về tỷ lệ hồi phục giữa hai nhóm)

# Việt hóa 
d$Nhap_vien <- factor(d$Nhap_vien, levels = c("No", "Yes"),
                      labels = c("Chưa nhập viện", "Đã nhập viện"))
d$Hoi_phuc <- factor(d$Hoi_phuc, levels = c("No", "Yes"),
                     labels = c("Chưa hồi phục", "Đã hồi phục"))

# Tạo bảng đếm
hr_tab2 <- table(d$Nhap_vien, d$Hoi_phuc)

# Truy xuất đúng tên nhãn
counts_hr2 <- c(
  hr_tab2["Chưa nhập viện", "Đã hồi phục"],
  hr_tab2["Đã nhập viện",   "Đã hồi phục"]
)

totals_hr2 <- c(
  sum(hr_tab2["Chưa nhập viện", ]),
  sum(hr_tab2["Đã nhập viện",   ])
)

# Kiểm định tỉ lệ
rd_test_hr2 <- prop.test(counts_hr2, totals_hr2, correct = FALSE)
print(rd_test_hr2)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  counts_hr2 out of totals_hr2
## X-squared = 0.25892, df = 1, p-value = 0.6109
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.02913260  0.04956446
## sample estimates:
##    prop 1    prop 2 
## 0.5056497 0.4954338

Nhận xét

Giá trị thống kê Chi-bình phương: X² = 0.25892
Bậc tự do (df) = 1
p-value = 0.6109
Khoảng tin cậy 95% cho p₂ − p₁: [-0.02913, 0.04956]
Ước lượng tỷ lệ: p₁ = 0.50565 (nhóm chưa tiêm), p₂ = 0.49543 (nhóm đã tiêm)
Hiệu tỷ lệ p₂ − p₁ ≈-0.01021 (khoảng -1,0%)

Với α = 0.05, p-value (0.6109) > 0.05 và khoảng tin cậy bao gồm 0, chúng ta không đủ cơ sở để bác bỏ H₀. Nói cách khác, hiệu chênh tỷ lệ hồi phục giữa nhóm chưa nhập viện và nhóm nhập viện (khoảng 1,0%) không có ý nghĩa thống kê. Điều này cho thấy trong mẫu 3.000 bệnh nhân, việc phải nhập viện không ảnh hưởng rõ rệt đến khả năng hồi phục.

Rủi ro tương đối (Relative Risk - RR):

Giả thuyết kiểm định

H₀: Tỷ số nguy cơ (RR) = 1 (Không có sự khác biệt về khả năng hồi phục giữa nhóm chưa nhập viện và nhóm nhập viện.)
H₁: Tỷ số nguy cơ (RR) ≠ 1 (Có sự khác biệt về khả năng hồi phục giữa hai nhóm.)

library(epitools)

rr_hr2 <- riskratio(hr_tab2, method = "wald")
print(rr_hr2)

## $data
##                 
##                  Chưa hồi phục Đã hồi phục Total
##   Chưa nhập viện          1050        1074  2124
##   Đã nhập viện             442         434   876
##   Total                   1492        1508  3000
## 
## $measure
##                 risk ratio with 95% C.I.
##                   estimate     lower    upper
##   Chưa nhập viện 1.0000000        NA       NA
##   Đã nhập viện   0.9797964 0.9054098 1.060295
## 
## $p.value
##                 two-sided
##                  midp.exact fisher.exact chi.square
##   Chưa nhập viện         NA           NA         NA
##   Đã nhập viện    0.6111487    0.6299499  0.6108632
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

Nhận xét

Ước lượng Tỷ số Nguy cơ (RR) của nhóm nhập viện so với nhóm chưa nhập viện là 0.9798.
Khoảng tin cậy 95% của RR: [0.9055, 1.0603], bao gồm giá trị 1.
p-value = 0.6111 > 0.05.

Vì p-value lớn hơn 0.05 và khoảng tin cậy chứa 1, chúng ta không đủ cơ sở để bác bỏ H₀. Điều này cho thấy trong mẫu 3.000 bệnh nhân, khả năng hồi phục không khác biệt có ý nghĩa thống kê giữa những người đã nhập viện và những người không nhập viện.

Tỷ số Chênh (Odds Ratio - OR):

or_hr2 <- oddsratio(hr_tab2)
print(or_hr2)

## $data
##                 
##                  Chưa hồi phục Đã hồi phục Total
##   Chưa nhập viện          1050        1074  2124
##   Đã nhập viện             442         434   876
##   Total                   1492        1508  3000
## 
## $measure
##                 odds ratio with 95% C.I.
##                   estimate     lower    upper
##   Chưa nhập viện 1.0000000        NA       NA
##   Đã nhập viện   0.9599776 0.8200695 1.123696
## 
## $p.value
##                 two-sided
##                  midp.exact fisher.exact chi.square
##   Chưa nhập viện         NA           NA         NA
##   Đã nhập viện    0.6111487    0.6299499  0.6108632
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Nhận xét

Ước lượng Odds Ratio (OR) cho “Hồi phục” ở nhóm nhập viện so với nhóm chưa nhập viện là 0.95998, cho thấy xu hướng odds hồi phục ở nhóm nhập viện thấp hơn khoảng 4.0%.
Khoảng tin cậy 95% cho OR: [0.8201, 1.1237], bao gồm giá trị 1.
p-value ≈ 0.6111 > 0.05.

Với α = 0.05, vì p-value = 0.6111 > 0.05 và khoảng tin cậy chứa 1, không đủ cơ sở để bác bỏ H₀. Điều này cho thấy trong mẫu 3.000 bệnh nhân, odds hồi phục không khác biệt có ý nghĩa thống kê giữa những người đã nhập viện và chưa nhập viện.

Tai_nhiem và Gioi_tinh

Bảng tần số và tần suất

# 2. Lập bảng tần số chéo
rg_tab  <- table(dt$Tai_nhiem, dt$Gioi_tinh)
rg_tab1 <- addmargins(rg_tab)
print(rg_tab1)

##      
##       Female Male  Sum
##   No    1374 1341 2715
##   Yes    153  132  285
##   Sum   1527 1473 3000

# 3. Bảng tần suất theo hàng (trong mỗi nhóm Tai_nhiem)
rg_prop         <- prop.table(rg_tab, margin = 1)
rg_prop_rounded <- round(rg_prop, 4)
print(rg_prop_rounded)

##      
##       Female   Male
##   No  0.5061 0.4939
##   Yes 0.5368 0.4632

Đồ thị

# Chuyển bảng tần số sang dạng data frame
rg_df <- as.data.frame(rg_tab)
colnames(rg_df) <- c("Tai_nhiem", "Gioi_tinh", "Count")

# Vẽ biểu đồ cột nhóm
library(ggplot2)

ggplot(rg_df, aes(x = Gioi_tinh, y = Count, fill = Tai_nhiem)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.7), width = 0.6) +
  labs(
    title = "Số bệnh nhân tái nhiễm theo giới tính",
    x     = "Giới tính",
    y     = "Số bệnh nhân",
    fill  = "Tái nhiễm"
  ) +
  theme_minimal() +
  scale_fill_manual(values = c(
    "Không tái nhiễm" = "yellow",
    "Tái nhiễm"       = "lightblue"
  ))

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.
## No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

Nhận xét

Trong 1.473 bệnh nhân nam, có 132 ca tái nhiễm (8,97%) và 1.341 ca không tái nhiễm (91,03%). Trong 1.527 bệnh nhân nữ, có 153 ca tái nhiễm (10,02%) và 1.374 ca không tái nhiễm (89,98%). Như vậy, tỷ lệ tái nhiễm ở nữ cao hơn nam khoảng 1,05 phần trăm, cho thấy phụ nữ trong mẫu này có xu hướng tái nhiễm nhẹ hơn so với nam giới.

Kiểm định Thống kê

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

H₀:Tình trạng tái nhiễm và giới tính là độc lập.
H₁:Có mối liên hệ giữa tái nhiễm và giới tính.

# Thực hiện kiểm định Chi bình phương
chi_rg <- chisq.test(rg_tab)
print(chi_rg)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  rg_tab
## X-squared = 0.85757, df = 1, p-value = 0.3544

Nhận xét

Giá trị thống kê Chi bình phương: X² = 0.85757
Bậc tự do (df) = 1
p-value = 0.3544

Với α = 0.05, p-value (0.3544) > 0.05, không đủ cơ sở để bác bỏ H₀.Vậy tình trạng tái nhiễm và giới tính là độc lập.

Nói cách khác, trong mẫu 3.000 bệnh nhân này, tỷ lệ tái nhiễm ở nam (132/1.473 ≈ 8,97 %) và nữ (153/1.527 ≈ 10,02 %) không khác biệt một cách có ý nghĩa thống kê. Sự chênh lệch khoảng 1,05 % là khá nhỏ, cho thấy giới tính không ảnh hưởng rõ rệt đến nguy cơ tái nhiễm.

Hiệu tỷ lệ (Risk Difference - RD)

Giả thuyết kiểm định

H₀: p₁ − p₂ = 0 (Tỷ lệ tái nhiễm ở nam bằng tỷ lệ tái nhiễm ở nữ)
H₁: p₁ − p₂ ≠ 0 (Có khác biệt về tỷ lệ tái nhiễm giữa nam và nữ)

# 1. Việt hóa các biến liên quan )
d$Tai_nhiem <- factor(d$Tai_nhiem, levels = c("No", "Yes"),
                      labels = c("Không tái nhiễm", "Tái nhiễm"))

d$Gioi_tinh <- factor(d$Gioi_tinh, levels = c("Male", "Female"),
                      labels = c("Nam", "Nữ"))

# 2. Tạo bảng đếm chéo giữa Tái nhiễm và Giới tính
rg_tab <- table(d$Tai_nhiem, d$Gioi_tinh)

# 3. Số ca "Tái nhiễm" theo giới tính
counts_rg <- c(
  rg_tab["Tái nhiễm", "Nam"],
  rg_tab["Tái nhiễm", "Nữ"]
)

# 4. Tổng số bệnh nhân mỗi giới tính
totals_rg <- c(
  sum(rg_tab[ , "Nam"]),
  sum(rg_tab[ , "Nữ"])
)

# 5. Kiểm định hiệu hai tỷ lệ (tỷ lệ tái nhiễm theo giới tính)
rd_test_rg <- prop.test(counts_rg, totals_rg, correct = FALSE)

# 6. In kết quả
print(rd_test_rg)

## 
##  2-sample test for equality of proportions without continuity correction
## 
## data:  counts_rg out of totals_rg
## X-squared = 0.97679, df = 1, p-value = 0.323
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  -0.03154930  0.01038244
## sample estimates:
##     prop 1     prop 2 
## 0.08961303 0.10019646

Nhận xét

Giá trị thống kê Chi-bình phương: X² = 0.97679
Bậc tự do (df) = 1
p-value = 0.3230
Khoảng tin cậy 95% cho p₁ − p₂: [−0.03155, 0.01038]
Ước lượng tỷ lệ: p₁ = 0.08961 (nam), p₂ = 0.10020 (nữ)
Hiệu tỷ lệ p₁ − p₂ ≈ −0.01059 (khoảng −1,06%)

Với α = 0.05, p-value (0.3230) > 0.05 và khoảng tin cậy bao gồm 0, chúng ta không đủ cơ sở để bác bỏ H₀. Điều này có nghĩa là tỷ lệ tái nhiễm giữa nam và nữ trong mẫu dữ liệu này không khác biệt một cách có ý nghĩa thống kê.

Rủi ro tương đối (Relative Risk - RR):

Giả thuyết kiểm định

H₀: Tỷ số nguy cơ (RR) = 1 (Không có sự khác biệt về giới tính nữ ở nhóm không tái nhiễm và nhóm tái nhiễm)
H₁: Tỷ số nguy cơ (RR) ≠ 1 (Có sự khác biệt về giới tính nữ ở nhóm không tái nhiễm và nhóm tái nhiễm)

library(epitools)

rr_rg <- riskratio(rg_tab, method = "wald")
print(rr_rg)

## $data
##                  
##                    Nam   Nữ Total
##   Không tái nhiễm 1341 1374  2715
##   Tái nhiễm        132  153   285
##   Total           1473 1527  3000
## 
## $measure
##                  risk ratio with 95% C.I.
##                   estimate     lower    upper
##   Không tái nhiễm 1.000000        NA       NA
##   Tái nhiễm       1.060791 0.9464421 1.188955
## 
## $p.value
##                  two-sided
##                   midp.exact fisher.exact chi.square
##   Không tái nhiễm         NA           NA         NA
##   Tái nhiễm        0.3239372     0.350262  0.3229926
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"

NHận xét

Ước lượng Tỷ số Nguy cơ (RR) của nhóm “Tái nhiễm” ở nữ so với nam là 1.0608, cho thấy xu hướng tỷ lệ tái nhiễm ở nữ cao hơn khoảng 6 % so với nam.
Khoảng tin cậy 95 % của RR: [0.9464, 1.1889], bao gồm giá trị 1.
p-value ≈ 0.3239 > 0.05.

Với mức ý nghĩa α = 0.05, p-value = 0.3239 lớn hơn 0.05 và khoảng tin cậy chứa 1, không đủ cơ sở để bác bỏ H₀.Vậy không có sự khác biệt về giới tính nữ ở nhóm không tái nhiễm và nhóm tái nhiễm

Tỷ số Chênh (Odds Ratio - OR):

# Tính Odds Ratio
or_rg <- oddsratio(rg_tab)
print(or_rg)

## $data
##                  
##                    Nam   Nữ Total
##   Không tái nhiễm 1341 1374  2715
##   Tái nhiễm        132  153   285
##   Total           1473 1527  3000
## 
## $measure
##                  odds ratio with 95% C.I.
##                   estimate     lower   upper
##   Không tái nhiễm  1.00000        NA      NA
##   Tái nhiễm        1.13103 0.8856823 1.44609
## 
## $p.value
##                  two-sided
##                   midp.exact fisher.exact chi.square
##   Không tái nhiễm         NA           NA         NA
##   Tái nhiễm        0.3239372     0.350262  0.3229926
## 
## $correction
## [1] FALSE
## 
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"

Nhận xét

Ước lượng Odds Ratio (OR) cho “Tái nhiễm” ở nhóm nữ so với nhóm nam là 1.1310, nghĩa là odds tái nhiễm ở nữ cao hơn khoảng 13.1%.
Khoảng tin cậy 95% cho OR: [0.8856, 1.4461], bao gồm giá trị 1.
p-value ≈ 0.3239 > 0.05.

Với α = 0.05, do p-value (0.3239) > 0.05 và khoảng tin cậy chứa 1, không đủ cơ sở để bác bỏ H₀. Điều này cho thấy trong mẫu 3.000 bệnh nhân, giới tính không có mối liên hệ có ý nghĩa thống kê với odds tái nhiễm COVID-19. # Mô hình Logit

Mô hình Logit với dữ liệu nhị phân

Mô hình Logit : Là một mô hình hồi quy dùng để dự đoán xác suất xảy ra của một sự kiện nhị phân (chỉ có 2 trạng thái: 1/0, Yes/No, Thành công/Thất bại).

Mục tiêu: Phân tích xem giới tính (Gender) có ảnh hưởng đến khả năng sinh viên thay đổi kế hoạch nghề nghiệp vì AI.

Biến phụ thuộc: career_plan_binary
1 = Yes (có thay đổi kế hoạch)
0 = No (không thay đổi kế hoạch)
Biến giải thích: gender_binary
0 = Male
1 = Female

Công thức tổng quát mô hình

\[ \log\left(\frac{P(\text{career\_plan\_binary}=1)}{1 - P(\text{career\_plan\_binary}=1)}\right) \]

# Chuẩn hóa tên cột: thay khoảng trắng bằng dấu chấm
names(d) <- gsub(" ", ".", names(d))

# Chỉ làm sạch dữ liệu nếu data frame không rỗng
if (nrow(d) > 0) {

  # Làm sạch cột Impact.of.AI.on.Career.Plans nếu tồn tại
  if ("Impact.of.AI.on.Career.Plans" %in% names(d)) {
    d$Impact.of.AI.on.Career.Plans <- trimws(tolower(as.character(d$Impact.of.AI.on.Career.Plans)))
  } else {
    warning("Cột Impact.of.AI.on.Career.Plans không tồn tại trong data frame")
    d$Impact.of.AI.on.Career.Plans <- NA  # hoặc bạn có thể bỏ qua
  }

  # Làm sạch cột Gender nếu tồn tại
  if ("Gender" %in% names(d)) {
    d$Gender <- trimws(tolower(as.character(d$Gender)))
  } else {
    warning("Cột Gender không tồn tại trong data frame")
    d$Gender <- NA
  }

  # Tạo biến nhị phân
  d$career_plan_binary <- if ("Impact.of.AI.on.Career.Plans" %in% names(d)) {
    ifelse(d$Impact.of.AI.on.Career.Plans == "yes", 1, 
           ifelse(is.na(d$Impact.of.AI.on.Career.Plans), NA, 0))
  } else {
    NA
  }

  d$gender_binary <- if ("Gender" %in% names(d)) {
    ifelse(d$Gender == "female", 1, 
           ifelse(is.na(d$Gender), NA, 0))
  } else {
    NA
  }

} else {
  warning("Data frame d rỗng: bỏ qua xử lý")
}

## Warning: Cột Impact.of.AI.on.Career.Plans không tồn tại trong data frame

## Warning: Cột Gender không tồn tại trong data frame

# Loại bỏ dòng có NA ở 2 cột cần thiết
d_clean <- na.omit(d[, c("career_plan_binary", "gender_binary")])

# Kiểm tra dữ liệu trước khi phân tích
if (nrow(d_clean) == 0) {
  warning("Không còn dòng dữ liệu sau khi loại NA. Bỏ qua phân tích.")
} else if (length(unique(d_clean$career_plan_binary)) < 2) {
  warning("Biến phụ thuộc 'career_plan_binary' chỉ có 1 giá trị. Không thể chạy logistic regression.")
} else {
  # Chạy mô hình logistic regression
  model <- glm(career_plan_binary ~ gender_binary, 
               data = d_clean, family = binomial)
  
  # Hiển thị kết quả
  summary(model)
  
  # Optional: tính odds ratio và CI nếu cần
  odds_ratios <- exp(coef(model))
  conf_int <- exp(confint(model))
  
  print("Odds Ratios:")
  print(odds_ratios)
  print("Confidence Intervals:")
  print(conf_int)
}

## Warning: Không còn dòng dữ liệu sau khi loại NA. Bỏ qua phân tích.

# Kiểm tra dữ liệu trước khi chạy logistic regression
if (nrow(d_clean) == 0) {
  warning("Dữ liệu rỗng sau khi loại NA. Bỏ qua phân tích.")
} else if (length(unique(d_clean$career_plan_binary)) < 2) {
  warning("Biến phụ thuộc 'career_plan_binary' chỉ có 1 mức giá trị. Không thể chạy logistic regression.")
} else if (length(unique(d_clean$gender_binary)) < 2) {
  warning("Biến độc lập 'gender_binary' chỉ có 1 mức giá trị. Không thể chạy logistic regression.")
} else {
  # Chạy mô hình logistic regression
  model_logit2 <- glm(career_plan_binary ~ gender_binary,
                      data = d_clean,
                      family = binomial(link = "logit"))
  
  # Hiển thị kết quả
  summary(model_logit2)
  
  # Optional: Tính odds ratio và khoảng tin cậy
  odds_ratios <- exp(coef(model_logit2))
  conf_int <- exp(confint(model_logit2))
  
  cat("\nOdds Ratios:\n")
  print(odds_ratios)
  cat("\nConfidence Intervals:\n")
  print(conf_int)
}

## Warning: Dữ liệu rỗng sau khi loại NA. Bỏ qua phân tích.

Nhận xét

Hệ số chặn (Intercept) là 0.5447, còn hệ số của biến gender_binary là -0.1163. Tuy nhiên, hệ số gender_binary không có ý nghĩa thống kê (p-value = 0.7028). Điều này cho thấy rằng giới tính không ảnh hưởng đáng kể đến khả năng sinh viên thay đổi kế hoạch nghề nghiệp do AI.

Công thức hồi quy logistic thu được từ mô hình là: \[ \log\left(\frac{P(\text{career\_plan} = 1)}{1 - P(\text{career\_plan} = 1)}\right) = 0.5447 - 0.1163 \cdot gender\_binary \]

Trong đó:

Xác suất và Biến giải thích

\(P(\text{career\_plan} = 1)\): Xác suất sinh viên thay đổi kế hoạch nghề nghiệp do AI.
\(gender\_binary = 0\): Nam, \(gender\_binary = 1\): Nữ.

Mức độ phù hợp của mô hình

Null deviance: 250.23
Thể hiện sai số log-likelihood của mô hình chỉ chứa hệ số chặn (Intercept). Đây là mức sai số khi chưa có biến giải thích nào được đưa vào mô hình.
Residual deviance: 250.08
Là sai số log-likelihood sau khi đưa biến \(gender\_binary\) vào mô hình. Giá trị này giảm rất ít (từ 250.23 xuống 250.08), cho thấy biến \(gender\_binary\) gần như không cải thiện khả năng dự đoán của mô hình.
AIC = 254.08
Chỉ số AIC (Akaike Information Criterion) dùng để so sánh mức độ phù hợp giữa các mô hình. Giá trị AIC ở đây khá cao, gợi ý rằng mô hình này còn đơn giản và có thể cải thiện bằng cách thêm các biến giải thích khác.

Đánh giá tổng quan

Sai số log-likelihood giảm không đáng kể khi thêm biến \(gender\_binary\), chứng tỏ biến này không đóng góp nhiều vào việc giải thích sự biến thiên của dữ liệu.
Giá trị AIC khá cao, cho thấy mô hình chưa tối ưu. Nên thử thêm các biến khác để cải thiện độ phù hợp.

Mô hình Logit với dữ liệu không phải nhị phân

Mục tiêu: phân tích xem việc đánh giá giới tính (Gender) có ảnh hưởng như thế nào đến mức độ lo lắng của sinh viên trước tác động của AI.

Biến phụ thuộc (Anxiety Level) có 4 mức độ phân loại:

No anxiety (không lo lắng)
Slight anxiety (hơi lo lắng)
Moderate anxiety (lo lắng vừa phải)
High anxiety (lo lắng cao)

Biến giải thích (Gender) được mã hóa nhị phân:

gender_binary = 0 nếu Male
gender_binary = 1 nếu Female

Vì biến phụ thuộc không phải là nhị phân, ta sử dụng mô hình hồi quy logistic đa thức (multinomial logistic regression) để mô hình hóa mối quan hệ.

Chọn No anxiety làm nhóm tham chiếu, mô hình hồi quy logistic đa thức được biểu diễn như sau:

\[ \log\left(\frac{P(Y = k)}{P(Y = \text{No anxiety})}\right) = \beta_{0k} + \beta_{1k} \cdot \text{gender\_binary} \]

Trong đó:

\(P(Y = k)\): Xác suất sinh viên thuộc mức độ lo lắng \(k\), với \(k \in \{\text{Slight anxiety}, \text{Moderate anxiety}, \text{High anxiety}\}\)
\(P(Y = \text{No anxiety})\): Xác suất sinh viên không lo lắng (nhóm tham chiếu)
\(\beta_{0k}\): Hệ số chặn (intercept) cho mức độ lo lắng \(k\)
\(\beta_{1k}\): Hệ số hồi quy cho biến gender_binary tại mức \(k\)
gender_binary: Biến độc lập, được mã hóa:

\[ gender\_binary = \begin{cases} 0, & \text{nếu giới tính là Nam}\\ 1, & \text{nếu giới tính là Nữ} \end{cases} \]

# Cài gói cần thiết nếu chưa có
if (!require(nnet, quietly = TRUE)) install.packages("nnet")
if (!require(broom, quietly = TRUE)) install.packages("broom")
library(nnet)
library(broom)

# Kiểm tra cột dữ liệu tồn tại
if (!"Anxiety.Level" %in% names(d)) {
  warning("Cột 'Anxiety.Level' không tồn tại trong data frame d. Bỏ qua phân tích.")
} else if (!"gender_binary" %in% names(d)) {
  warning("Cột 'gender_binary' không tồn tại trong data frame d. Bỏ qua phân tích.")
} else {
  # Loại bỏ dòng có NA
  d_clean <- na.omit(d[, c("Anxiety.Level", "gender_binary")])
  
  # Kiểm tra dữ liệu hợp lệ
  if (nrow(d_clean) == 0) {
    warning("Dữ liệu rỗng sau khi loại NA. Không thể chạy mô hình multinom.")
  } else if (length(unique(d_clean$Anxiety.Level)) < 2) {
    warning("Biến phụ thuộc 'Anxiety.Level' chỉ có 1 mức giá trị. Không thể chạy mô hình multinom.")
  } else if (length(unique(d_clean$gender_binary)) < 2) {
    warning("Biến độc lập 'gender_binary' chỉ có 1 mức giá trị. Không thể chạy mô hình multinom.")
  } else {
    # Chạy mô hình multinomial logistic regression
    model_multinom <- multinom(Anxiety.Level ~ gender_binary, data = d_clean, trace = FALSE)
    
    # Tóm tắt mô hình
    summary_model <- summary(model_multinom)
    
    # Tính z-value và p-value
    z_values <- summary_model$coefficients / summary_model$standard.errors
    p_values <- 2 * (1 - pnorm(abs(z_values)))
    
    # Gộp kết quả thành bảng đẹp
    result_table <- as.data.frame(cbind(
      Estimate = round(summary_model$coefficients, 4),
      `Std. Error` = round(summary_model$standard.errors, 4),
      `z value` = round(z_values, 3),
      `Pr(>|z|)` = signif(p_values, 3)
    ))
    
    # Hiển thị bảng kết quả
    print(result_table)
    
    # Tính Null Deviance (log-likelihood của mô hình chỉ có intercept)
    model_null <- multinom(Anxiety.Level ~ 1, data = d_clean, trace = FALSE)
    
    null_deviance <- -2 * logLik(model_null)
    residual_deviance <- -2 * logLik(model_multinom)
    df_null <- attr(logLik(model_null), "df")
    df_residual <- attr(logLik(model_multinom), "df")
    
    # AIC
    aic_value <- AIC(model_multinom)
    
    # Hiển thị thông tin deviance & AIC
    cat("\nNull deviance: ", round(null_deviance, 3), " on ", df_null, " degrees of freedom\n")
    cat("Residual deviance: ", round(residual_deviance, 3), " on ", df_residual, " degrees of freedom\n")
    cat("AIC: ", round(aic_value, 3), "\n")
    cat("Number of Fisher Scoring iterations: ", model_multinom$iter, "\n")
  }
}

## Warning: Cột 'Anxiety.Level' không tồn tại trong data frame d. Bỏ qua phân
## tích.

Ý nghĩa

Anxiety Level	Intercept (β₀k)	gender_binary (β₁k)	p-value gender_binary
Moderate anxiety	0.4520	1.8506	0.0366 (*)
No anxiety	1.4553	1.4625	0.0812 (.)
Slight anxiety	1.4881	1.7308	0.0379 (*)

Intercept (β₀k): log-odds của từng mức Anxiety khi gender_binary = 0 (nam).
gender_binary (β₁k): thay đổi log-odds khi chuyển từ nam (gender_binary = 0) sang nữ (gender_binary = 1).

Chuyển đổi sang Odds Ratio (OR):

Để hiểu dễ hơn, ta lấy exp(β₁k)

Anxiety Level	β₁k	exp(β₁k) (OR)	Ý nghĩa
Moderate anxiety	1.8506	6.36	Nữ có gấp ~6.4 lần odds rơi vào Moderate anxiety so với nam
No anxiety	1.4625	4.31	Nữ có odds cao hơn ~4.3 lần rơi vào No anxiety so với nam
Slight anxiety	1.7308	5.65	Nữ có odds cao hơn ~5.7 lần rơi vào Slight anxiety so với nam

Ý nghĩa

Moderate anxiety: p = 0.0366 (< 0.05) → có ý nghĩa thống kê.
Slight anxiety: p = 0.0379 (< 0.05) → có ý nghĩa thống kê.
No anxiety: p = 0.0812 (> 0.05 nhưng < 0.1) → gần có ý nghĩa (marginal significance).

=> Như vậy, giới tính có ảnh hưởng đáng kể đến khả năng rơi vào Slight/Moderate anxiety. Với No anxiety thì yếu hơn.

Bảng thông số

Thông số	Giá trị	Ý nghĩa
Null deviance	441.115 (df=3)	Sai số log-likelihood của mô hình chỉ có Intercept (chưa có gender_binary).
Residual deviance	435.099 (df=6)	Sai số log-likelihood sau khi thêm gender_binary. Giảm → mô hình cải thiện.
AIC	447.099	Chỉ số đánh giá mô hình. Thấp hơn nghĩa là mô hình tốt hơn.
Fisher iterations	10	Số lần lặp để tối ưu hóa mô hình.

Sự cải thiện của mô hình - Null deviance – Residual deviance = 441.115 – 435.099 = 6.016

Mức giảm nhỏ → gender_binary chỉ giải thích một phần nhỏ biến thiên của Anxiety Level.

Kết luận

Gender có ảnh hưởng đến Anxiety Level, đặc biệt ở Slight và Moderate anxiety.
Nữ có odds cao hơn đáng kể so với nam rơi vào Slight và Moderate anxiety.
Tuy nhiên, mức độ cải thiện của mô hình tổng thể (giảm deviance) không quá lớn, cần thêm các biến khác (ví dụ: tuổi, năm học) để cải thiện dự đoán.

Mô hình Probit

Mục tiêu Phân tích xem giới tính của sinh viên có ảnh hưởng đến xác suất rơi vào từng mức độ lo lắng (Anxiety Level) hay không, sử dụng mô hình Probit đa thức (multinomial Probit).

Dữ liệu

Biến phụ thuộc: anxiety_level (4 mức độ phân loại: No anxiety, Slight anxiety, Moderate anxiety, High anxiety)
Biến giải thích: gender_binary (mã hóa: 0 = Nữ (Female), 1 = Nam (Male))

Công thức tổng quát của mô hình Probit \[ P(Y=1 \mid Gender) = \Phi(\beta_0 + \beta_1 \cdot Gender) \]

Trong đó:

# Tạo biến nhị phân: 1 nếu "High anxiety", 0 nếu còn lại
if (!"Anxiety.Level" %in% names(d)) {
  warning("Cột 'Anxiety.Level' không tồn tại. Gán anxiety_binary = NA.")
  d$anxiety_binary <- NA
} else {
  temp <- trimws(tolower(as.character(d$Anxiety.Level)))
  
  if (length(temp) == nrow(d)) {
    d$anxiety_binary <- ifelse(temp == "high anxiety", 1, 0)
  } else {
    warning("Số dòng của 'Anxiety.Level' không khớp với data frame. Gán anxiety_binary = NA.")
    d$anxiety_binary <- NA
  }
}

## Warning: Cột 'Anxiety.Level' không tồn tại. Gán anxiety_binary = NA.

# Chạy mô hình Probit nếu dữ liệu hợp lệ
if (!"Gender" %in% names(d)) {
  warning("Cột 'Gender' không tồn tại. Bỏ qua mô hình Probit.")
} else {
  d_clean <- na.omit(d[, c("anxiety_binary", "Gender")])
  
  if (nrow(d_clean) == 0) {
    warning("Dữ liệu rỗng sau khi loại NA. Không thể chạy mô hình Probit.")
  } else if (length(unique(d_clean$anxiety_binary)) < 2) {
    warning("Biến phụ thuộc 'anxiety_binary' chỉ có 1 mức giá trị. Không thể chạy mô hình Probit.")
  } else if (length(unique(d_clean$Gender)) < 2) {
    warning("Biến độc lập 'Gender' chỉ có 1 mức giá trị. Không thể chạy mô hình Probit.")
  } else {
    # Chạy mô hình Probit
    model_probit <- glm(anxiety_binary ~ Gender, 
                        data = d_clean, 
                        family = binomial(link = "probit"))
    
    # Hiển thị kết quả
    summary(model_probit)
  }
}

## Warning: Dữ liệu rỗng sau khi loại NA. Không thể chạy mô hình Probit.

Ý nghĩa các hệ số

Biến	Estimate	Std. Error	z value	p-value	Ý nghĩa
Intercept	-2.0891	0.2857	-7.313	2.61e-13***	Log-odds của Nữ (mặc định trong R là baseline) rơi vào nhóm “High anxiety”.
Gendermale	0.7397	0.3483	2.124	0.0337 *	So với Nữ, Nam có log-odds cao hơn 0.7397 để rơi vào “High anxiety”.

Vì p-value < 0.05:

Hệ số Gendermale có ý nghĩa thống kê ở mức 5%.
Nam có xác suất cao hơn nữ để rơi vào mức High anxiety.

Thống kê mô hình

Thống kê	Giá trị	Ý nghĩa
Null deviance	72.268	Sai số khi không có biến giải thích.
Residual deviance	67.246	Sai số sau khi thêm biến Gender.
Giảm deviance	5.022	Mô hình cải thiện đáng kể.
AIC	71.246	Để so sánh mức phù hợp với các mô hình khác.

NHIỆM VỤ TUẦN 5

Võ Hoàng Tính

2025-06-15

PHẦN 1. GIỚI THIỆU VỀ BỘ DỮ LIỆU

PHẦN 2. PHÂN TÍCH MỘT BIẾN ĐỊNH TÍNH

Biến Nhap_vien

Biến Nghe_nghiep

Biến Trieu_chung

Biến Hoi_phuc

Biến Tai_nhiem

Biến Gioi_tinh

Biến Tiem_chung

Biến Khu_vuc

Biến Vaccine_Type

PHẦN 3. ƯỚC LƯỢNG KHOẢNG TIN CẬY

Nhap_vien (Nhập viện = “Yes”)

Tai_nhiem (Tái nhiễm = “Yes”)

Tiem_chung (Đã tiêm = “Yes”)

PHẦN 4. PHÂN TÍCH MỐI QUAN HỆ GIỮA HAI BIẾN ĐỊNH TÍNH

Nhap_vien và Hoi_phuc

Bảng tần số và tần suất

Kiểm định Thống kê

Hiệu tỷ lệ (Risk Difference - RD)

Rủi ro tương đối (Relative Risk - RR):

Tỷ số Chênh (Odds Ratio - OR):

Tai_nhiem và Gioi_tinh

Bảng tần số và tần suất

Kiểm định Thống kê

Hiệu tỷ lệ (Risk Difference - RD)

Rủi ro tương đối (Relative Risk - RR):

Tỷ số Chênh (Odds Ratio - OR):

Mô hình Logit với dữ liệu nhị phân

Mô hình Logit với dữ liệu không phải nhị phân

Mô hình Probit