1 Mục tiêu: Phân tích các yếu tố ảnh hưởng đến khả năng tăng huyết áp và bị bệnh tiểu đường

2 Nhập dữ liệu

library(tidyverse)
library(xlsx)
data <- read.xlsx(file = "C:/Users/Khanh Hoa/OneDrive/UFM/R/Tuần 4 - PTDLDT/diabetes_data.xlsx", header = T, sheetIndex = 1) #load data
str(data)
## 'data.frame':    1879 obs. of  15 variables:
##  $ Age                          : num  44 51 89 21 27 65 61 74 54 82 ...
##  $ Gender                       : num  0 1 1 1 1 0 1 1 0 1 ...
##  $ Ethnicity                    : num  1 0 0 1 0 0 2 3 0 0 ...
##  $ SocioeconomicStatus          : num  2 1 1 1 1 0 1 0 1 1 ...
##  $ EducationLevel               : num  1 2 3 2 3 0 3 3 2 1 ...
##  $ Smoking                      : num  1 0 0 1 0 1 0 0 0 1 ...
##  $ FamilyHistoryDiabetes        : num  1 0 1 1 0 0 0 0 0 0 ...
##  $ PreviousPreDiabetes          : num  0 0 0 1 0 0 0 0 0 0 ...
##  $ Hypertension                 : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ AntidiabeticMedications      : num  1 0 0 1 0 0 0 0 0 0 ...
##  $ FrequentUrination            : num  0 0 0 0 0 0 1 0 0 0 ...
##  $ ExcessiveThirst              : num  0 0 0 0 0 0 0 0 0 1 ...
##  $ OccupationalExposureChemicals: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ WaterQuality                 : num  0 1 0 1 0 0 0 0 0 0 ...
##  $ Diagnosis                    : num  1 1 0 0 0 0 0 0 1 0 ...

Tập dữ liệu diabetes_data.xlsx chứa thông tin về bệnh nhân và các yếu tố liên quan đến bệnh tiểu đường. Dưới đây là giải thích chi tiết về các biến trong tập dữ liệu:

Biến phụ thuộc:

  • Hypertension (Tăng huyết áp):
    Mã hóa:
    0 = Không,
    1 = Có.

  • Diagnosis (Chẩn đoán tiểu đường - Biến mục tiêu):
    Mã hóa:
    0 = Không mắc bệnh,
    1 = Mắc bệnh tiểu đường.

Biến độc lập:

  • Age (Tuổi):
    Tuổi của bệnh nhân dao động từ 20 đến 90 tuổi.

  • Gender (Giới tính):
    Giới tính của bệnh nhân, mã hóa:
    0 = Nam,
    1 = Nữ.

  • Ethnicity (Dân tộc):
    Dân tộc của bệnh nhân, mã hóa:
    0 = Caucasian,
    1 = African American,
    2 = Asian,
    3 = Khác.

  • SocioeconomicStatus (Tình trạng kinh tế xã hội):
    Tình trạng kinh tế xã hội của bệnh nhân, mã hóa:
    0 = Thấp,
    1 = Trung bình,
    2 = Cao.

  • EducationLevel (Trình độ học vấn):
    Trình độ học vấn của bệnh nhân, mã hóa:
    0 = Không học,
    1 = Trung học phổ thông,
    2 = Đại học,
    3 = Sau đại học hoặc cao hơn.

  • Smoking (Hút thuốc):
    Tình trạng hút thuốc của bệnh nhân, mã hóa:
    0 = Không hút thuốc,
    1 = Có hút thuốc.

  • FamilyHistoryDiabetes (Tiền sử gia đình mắc bệnh tiểu đường):
    Mã hóa:
    0 = Không có,
    1 = Có.

  • PreviousPreDiabetes (Tiền sử tiền tiểu đường):
    Mã hóa:
    0 = Không có,
    1 = Có.

  • AntidiabeticMedications (Sử dụng thuốc tiểu đường):
    Mã hóa:
    0 = Không sử dụng,
    1 = Có sử dụng.

  • FrequentUrination (Đi tiểu thường xuyên):
    Mã hóa:
    0 = Không,
    1 = Có.

  • ExcessiveThirst (Khát nước quá mức):
    Mã hóa:
    0 = Không,
    1 = Có.

  • OccupationalExposureChemicals (Tiếp xúc với hóa chất nghề nghiệp):
    Mã hóa:
    0 = Không,
    1 = Có.

  • WaterQuality (Chất lượng nước):
    Mã hóa:
    0 = Tốt,
    1 = Kém.

3 Thống kê mô tả các biến định tính

#create function to calculate frequency and corresponding rate
thong_ke_dinh_tinh <- function(data, var_name) {
  
  #Convert to factor if needed
  variable <- as.factor(data[[var_name]])
  
  #freq
  freq <- table(variable)
  percent <- prop.table(freq) * 100
  
  #result
  result <- data.frame(
    Value = names(freq),
    Frequency = as.vector(freq),
    Percent = round(as.vector(percent), 2)
  )
  return(result)
}

library(tidyverse)
library(kableExtra)
thong_ke_dinh_tinh(data, "Gender") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 1: Thống kê mô tả biến Gender</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 1: Thống kê mô tả biến Gender
Value Frequency Percent
0 963 51.25
1 916 48.75

Dựa trên kết quả thống kê được cung cấp, có thể đánh giá rằng sự phân bổ giới tính trong tập dữ liệu này là rất cân bằng và có tính đại diện tốt. Cụ thể, số lượng bệnh nhân nam (mã hóa là 0) là 963, chiếm 51.25%, trong khi số lượng bệnh nhân nữ (mã hóa là 1) là 916, chiếm 48.75%. Với tỷ lệ gần như 50/50 và sự chênh lệch chỉ khoảng 2.5%, tập dữ liệu không cho thấy sự thiên lệch đáng kể về một giới tính nào.

thong_ke_dinh_tinh(data, "Ethnicity") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 2: Thống kê mô tả biến Ethnicity</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 2: Thống kê mô tả biến Ethnicity
Value Frequency Percent
0 1099 58.49
1 357 19.00
2 206 10.96
3 217 11.55

Phân tích thống kê về biến dân tộc cho thấy một sự mất cân đối rõ rệt trong thành phần của tập dữ liệu, với nhóm Caucasian (mã hóa 0) chiếm đa số áp đảo với 1099 người (58.49%). Trong khi đó, các nhóm dân tộc khác có số lượng mẫu ít hơn đáng kể, bao gồm nhóm African American (19.00%), nhóm Khác (11.55%) và nhóm Asian (10.96%). Sự chênh lệch lớn này tạo ra một sự thiên lệch (skew) đáng kể, dẫn đến nguy cơ “sai lệch do mất cân đối”.

thong_ke_dinh_tinh(data, "SocioeconomicStatus") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 3: Thống kê mô tả biến SocioeconomicStatus</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 3: Thống kê mô tả biến SocioeconomicStatus
Value Frequency Percent
0 557 29.64
1 780 41.51
2 542 28.85

Dựa trên kết quả thống kê, có thể đánh giá rằng sự phân bổ về tình trạng kinh tế xã hội (SocioeconomicStatus) trong tập dữ liệu là tương đối cân bằng và hợp lý. Nhóm có tình trạng kinh tế xã hội “Trung bình” (mã hóa 1) chiếm tỷ lệ cao nhất với 41.51% (780 người), trong khi hai nhóm “Thấp” (mã hóa 0) và “Cao” (mã hóa 2) có tỷ lệ gần như nhau, lần lượt là 29.64% (557 người) và 28.85% (542 người). Vì không có nhóm nào bị thiếu hụt mẫu một cách nghiêm trọng, nguy cơ “sai lệch do mất cân đối” đối với biến này là thấp, cho phép các phân tích và mô hình có thể rút ra những kết luận đáng tin cậy về ảnh hưởng của từng mức kinh tế xã hội đến chẩn đoán bệnh tiểu đường.

thong_ke_dinh_tinh(data, "EducationLevel") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 4: Thống kê mô tả biến EducationLevel</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 4: Thống kê mô tả biến EducationLevel
Value Frequency Percent
0 163 8.67
1 615 32.73
2 725 38.58
3 376 20.01

Kết quả phân tích thống kê cho biến trình độ học vấn (EducationLevel) cho thấy sự tập trung rõ rệt của các bệnh nhân ở các bậc học cao. Cụ thể, nhóm có trình độ “Đại học” (mã hóa 2) chiếm tỷ lệ cao nhất với 38.58%, theo sau là nhóm “Trung học phổ thông” (mã hóa 1) với 32.73%. Hai nhóm này chiếm hơn 71% tổng số mẫu, cho thấy phần lớn những người tham gia khảo sát đều có trình độ học vấn từ trung học trở lên. Nhóm “Sau đại học” (mã hóa 3) cũng chiếm một tỷ lệ đáng kể là 20.01%, trong khi nhóm “Không học” (mã hóa 0) là nhóm nhỏ nhất, chỉ chiếm 8.67%.

thong_ke_dinh_tinh(data, "Smoking") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 5: Thống kê mô tả biến Smoking</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 5: Thống kê mô tả biến Smoking
Value Frequency Percent
0 1350 71.85
1 529 28.15

Đánh giá về tình trạng hút thuốc trong tập dữ liệu cho thấy một sự mất cân đối đáng kể, với phần lớn người tham gia không hút thuốc. Cụ thể, có tới 1350 người được ghi nhận là “Không hút thuốc” (mã hóa 0), chiếm 71.85% tổng số mẫu, trong khi chỉ có 529 người “Có hút thuốc” (mã hóa 1), chiếm 28.15%. Sự chênh lệch này, với tỷ lệ người không hút thuốc cao gấp khoảng 2.5 lần người hút thuốc, cho thấy dữ liệu bị thiên lệch về nhóm không hút thuốc

thong_ke_dinh_tinh(data, "FamilyHistoryDiabetes") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 6: Thống kê mô tả biến FamilyHistoryDiabetes</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 6: Thống kê mô tả biến FamilyHistoryDiabetes
Value Frequency Percent
0 1431 76.16
1 448 23.84

Dữ liệu về tiền sử gia đình mắc bệnh tiểu đường cho thấy một sự mất cân đối rất lớn, với phần lớn các cá nhân trong mẫu khảo sát không có tiền sử bệnh này trong gia đình. Cụ thể, có tới 1431 người (chiếm 76.16%) được ghi nhận không có tiền sử gia đình mắc bệnh, trong khi chỉ có 448 người (chiếm 23.84%) có tiền sử bệnh này. Với việc nhóm không có tiền sử bệnh chiếm đa số áp đảo, gấp hơn ba lần so với nhóm còn lại, tập dữ liệu này thể hiện sự thiên lệch rõ rệt.

thong_ke_dinh_tinh(data, "PreviousPreDiabetes") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 7: Thống kê mô tả biến PreviousPreDiabetes</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 7: Thống kê mô tả biến PreviousPreDiabetes
Value Frequency Percent
0 1590 84.62
1 289 15.38

Phân tích về tiền sử tiền tiểu đường cho thấy một trong những sự mất cân đối nghiêm trọng nhất trong tập dữ liệu này. Cụ thể, có đến 1590 người, chiếm 84.62% tổng số, không có tiền sử về tiền tiểu đường, trong khi chỉ có 289 người, tương đương 15.38%, có tiền sử này. Sự chênh lệch cực kỳ lớn này, với nhóm không có tiền sử nhiều hơn gấp 5.5 lần, tạo ra một sự thiên vị rất mạnh mẽ. Điều này đặc biệt đáng lo ngại cho việc xây dựng mô hình dự đoán, vì mô hình có thể dễ dàng đạt được độ chính xác cao bằng cách mặc định dự đoán “không có tiền sử tiền tiểu đường”, và do đó có nguy cơ bỏ qua các đặc điểm quan trọng của nhóm thiểu số, vốn là nhóm có nguy cơ cao và cần được quan tâm đặc biệt trong phân tích y tế.

thong_ke_dinh_tinh(data, "Hypertension") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 8: Thống kê mô tả biến Hypertension</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 8: Thống kê mô tả biến Hypertension
Value Frequency Percent
0 1591 84.67
1 288 15.33

Tương tự như biến tiền sử tiền tiểu đường, dữ liệu về tăng huyết áp (Hypertension) cũng cho thấy một sự mất cân đối rất lớn. Phần lớn những người tham gia, với 1591 người (chiếm 84.67%), không bị tăng huyết áp, trong khi chỉ có một nhóm thiểu số nhỏ gồm 288 người (15.33%) mắc bệnh này.

thong_ke_dinh_tinh(data, "AntidiabeticMedications") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 9: Thống kê mô tả biến AntidiabeticMedications</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 9: Thống kê mô tả biến AntidiabeticMedications
Value Frequency Percent
0 1356 72.17
1 523 27.83

Phần lớn người tham gia, với 1356 người (chiếm 72.17%), không sử dụng thuốc điều trị tiểu đường, trong khi chỉ có 523 người (chiếm 27.83%) có sử dụng. Sự chênh lệch này, với nhóm không dùng thuốc lớn hơn gần 2.6 lần, tạo ra một sự thiên vị trong dữ liệu.

thong_ke_dinh_tinh(data, "FrequentUrination") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 10: Thống kê mô tả biến FrequentUrination</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 10: Thống kê mô tả biến FrequentUrination
Value Frequency Percent
0 1508 80.26
1 371 19.74

Dữ liệu về triệu chứng đi tiểu thường xuyên cho thấy một sự mất cân bằng lớn, khi phần lớn người tham gia không gặp phải triệu chứng này. Cụ thể, có 1508 người (chiếm 80.26%) không đi tiểu thường xuyên, trong khi chỉ có 371 người (19.74%) báo cáo có triệu chứng này.

thong_ke_dinh_tinh(data, "ExcessiveThirst") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 11: Thống kê mô tả biến ExcessiveThirst</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 11: Thống kê mô tả biến ExcessiveThirst
Value Frequency Percent
0 1516 80.68
1 363 19.32

Tương tự như triệu chứng đi tiểu thường xuyên, dữ liệu về triệu chứng khát nước quá mức (ExcessiveThirst) cũng cho thấy một sự mất cân bằng lớn. Phần lớn người tham gia, với 1516 người (chiếm 80.68%), không gặp phải tình trạng này, trong khi chỉ có 363 người (19.32%) có triệu chứng

thong_ke_dinh_tinh(data, "OccupationalExposureChemicals") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 12: Thống kê mô tả biến OccupationalExposureChemicals</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 12: Thống kê mô tả biến OccupationalExposureChemicals
Value Frequency Percent
0 1685 89.68
1 194 10.32

Xét về yếu tố phơi nhiễm hóa chất nghề nghiệp, dữ liệu cho thấy đây là một đặc điểm tương đối hiếm gặp trong quần thể được nghiên cứu. Một nhóm nhỏ chỉ gồm 194 cá nhân (chiếm 10.32%) cho biết có tiếp xúc, và nhóm này gần như bị “chìm nghỉm” so với đại đa số 1685 người (89.68%) không phơi nhiễm.

thong_ke_dinh_tinh(data, "WaterQuality") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 13: Thống kê mô tả biến WaterQuality</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 13: Thống kê mô tả biến WaterQuality
Value Frequency Percent
0 1502 79.94
1 377 20.06

Phân tích về chất lượng nước cho thấy một sự chênh lệch đáng kể trong trải nghiệm của những người tham gia, với phần lớn (79.94%) báo cáo sử dụng nước chất lượng tốt. Nhóm thiểu số, chiếm 20.06%, những người sử dụng nước chất lượng kém, bị lấn át hoàn toàn bởi nhóm đa số

thong_ke_dinh_tinh(data, "Diagnosis") %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 14: Thống kê mô tả biến Diagnosis</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 14: Thống kê mô tả biến Diagnosis
Value Frequency Percent
0 1127 59.98
1 752 40.02

Dữ liệu chẩn đoán, là biến mục tiêu của tập dữ liệu, cho thấy một sự phân bổ tương đối khả quan cho việc xây dựng mô hình dự đoán. Cụ thể, có 1127 trường hợp được chẩn đoán “Không mắc bệnh” (chiếm 59.98%) và 752 trường hợp “Mắc bệnh tiểu đường” (chiếm 40.02%). Mặc dù vẫn có sự mất cân bằng với tỷ lệ gần 60/40, mức độ này không quá nghiêm trọng và được xem là khá tốt trong nhiều bài toán thực tế. Quan trọng nhất, với 752 mẫu trong lớp thiểu số (nhóm mắc bệnh), mô hình sẽ có đủ dữ liệu để học hỏi và nhận diện các đặc điểm của nhóm này một cách hiệu quả.

4 Ước lượng khoảng và Kiểm định tỷ lệ

4.1 Ước lượng khoảng cho 2 biến phụ thuộc

Ước lượng khoảng là một kỹ thuật thống kê dùng để ước lượng khoảng giá trị có khả năng chứa tham số thực (như trung bình hoặc tỷ lệ) của tổng thể, dựa trên mẫu quan sát được. Thay vì đưa ra một giá trị duy nhất như ước lượng điểm, ta cung cấp một khoảng giá trị với một mức tin cậy cụ thể.

\[ CI = \hat{p} \pm z_{\alpha/2} \cdot \sqrt{ \frac{ \hat{p}(1 - \hat{p}) }{n} } \]

Trong đó:

  • \(\hat{p}\): Tỷ lệ mẫu

  • \(n\): Kích thước mẫu

  • \(z_{\alpha/2}\): Giá trị tới hạn từ phân phối chuẩn

prop_int_est <- function(data, var_name, value, alpha) {
  #execute
  n <- nrow(data)
  z <- qnorm(1-alpha/2)
  obs_p <- sum(data[[var_name]] == value)/n
  lower_bound <- obs_p - z*sqrt(obs_p*(1-obs_p)/n)
  upper_bound <- obs_p + z*sqrt(obs_p*(1-obs_p)/n)
  result <- data.frame(
    Value = paste0("Giá trị ",value, " trong ", var_name),
    Lower_bound = round(lower_bound,3),
    Upper_bound = round(upper_bound,3),
    Point_est = round(obs_p,3)
  )
  #return
  return(result)
}

prop_int_est(data, var_name = "Hypertension",value = 1,alpha = 0.05) %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 15: Khoảng ước lượng tỷ lệ tăng huyết áp</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 15: Khoảng ước lượng tỷ lệ tăng huyết áp
Value Lower_bound Upper_bound Point_est
Giá trị 1 trong Hypertension 0.137 0.17 0.153

Với độ tin cậy 95% (tương ứng alpha = 0.05), chúng ta có thể phát biểu rằng tỷ lệ thực sự của những người bị tăng huyết áp trong quần thể lớn hơn nằm trong khoảng từ 13.7% (Lower_bound) đến 17.0% (Upper_bound). Nói một cách đơn giản, thay vì chỉ khẳng định tỷ lệ là 15.3%, chúng ta có thể tự tin đến 95% rằng con số thực tế trong dân số nói chung sẽ rơi vào đâu đó trong khoảng này.

prop_int_est(data, var_name = "Diagnosis",value = 1,alpha = 0.05) %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 16: Khoảng ước lượng tỷ lệ bị bệnh tiểu đường</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 16: Khoảng ước lượng tỷ lệ bị bệnh tiểu đường
Value Lower_bound Upper_bound Point_est
Giá trị 1 trong Diagnosis 0.378 0.422 0.4

Với độ tin cậy 95% (alpha = 0.05), khoảng tin cậy cho tỷ lệ người thực sự mắc bệnh tiểu đường trong toàn bộ quần thể là từ 37.8% (Lower_bound) đến 42.2% (Upper_bound). Điều này có nghĩa là, chúng ta có thể rất tin tưởng rằng tỷ lệ thực tế của bệnh tiểu đường trong dân số mà mẫu này đại diện không phải là một con số cố định 40%, mà nó dao động trong một khoảng hẹp và hợp lý từ gần 38% đến hơn 42%. Khoảng tin cậy này cung cấp một bức tranh thực tế hơn về mức độ phổ biến của bệnh, thừa nhận sự không chắc chắn tự nhiên khi làm việc với dữ liệu mẫu.

4.2 Kiểm định tỷ lệ

  • H0: Tỷ lệ người không bị tăng huyết áp bằng với số người bị tăng huyết áp
one_spl_pro_test <- function(data, var_name, value, p, alpha) {
  test <- prop.test(x = sum(data[[var_name]] == value), n = nrow(data), p = p, conf.level = 1-alpha)
  result <- data.frame(
    Statistic = round(test$statistic,3),
    P_value = round(test$p.value,3),
    Result = if(test$p.value > alpha) {"Không đủ cơ sở bác bỏ H0"} else {"Bác bỏ H0"}
  )
  return(result)
}

a <- nrow(data[data$Hypertension == 0,])/nrow(data)
one_spl_pro_test(data, var_name = "Hypertension",value = 1 ,p = a,alpha = 0.05) %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 17: Kết quả kiểm định</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 17: Kết quả kiểm định
Statistic P_value Result
X-squared 6956.961 0 Bác bỏ H0

Với mức ý nghĩa 5%, kết quả kiểm định cho thấy tỷ lệ người không bị tăng huyết áp khác với tỷ lệ người bị tăng huyết áp

5 Tác động của hút thuốc đến việc bị bệnh tiểu đường

5.1 Bảng tần số chéo

phan_tich_cap_bien <- function(data, var1, var2) {
  
  tbl <- table(data[[var1]], data[[var2]])
  
  percent_total <- round(prop.table(tbl) * 100, 2)
  percent_row <- round(prop.table(tbl, margin = 1) * 100, 2)
  percent_col <- round(prop.table(tbl, margin = 2) * 100, 2)
  
  chi_result <- chisq.test(tbl)
  
  result <- list(
    Bang_tan_suat = tbl,
    Ti_le_theo_tong = percent_total,
    Ti_le_theo_hang = percent_row,
    Ti_le_theo_cot = percent_col,
    Chi_squared_test = chi_result
  )
  
  return(result)
}
result_hp <- phan_tich_cap_bien(data, "Smoking", "Diagnosis")
#convert to data frame
table1 <- as.data.frame(as.table(result_hp$Bang_tan_suat)) %>%
  rename(Smoking = Var1, Diagnosis = Var2, Frequency = Freq) %>%
  mutate(
    Percent_Total = as.vector(result_hp$Ti_le_theo_tong),
    Percent_of_Smoking   = as.vector(result_hp$Ti_le_theo_hang),
    Percent_Diagnosis   = as.vector(result_hp$Ti_le_theo_cot)
  )
table1 %>% 
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 18: Bảng tần số chéo giữa Smoking x Diagnosis</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 18: Bảng tần số chéo giữa Smoking x Diagnosis
Smoking Diagnosis Frequency Percent_Total Percent_of_Smoking Percent_Diagnosis
0 0 832 44.28 61.63 73.82
1 0 295 15.70 55.77 26.18
0 1 518 27.57 38.37 68.88
1 1 234 12.45 44.23 31.12
data.frame(
  Chi_Squared = round(result_hp$Chi_squared_test$statistic, 3),
  DF = result_hp$Chi_squared_test$parameter,
  P_Value = round(result_hp$Chi_squared_test$p.value, 4)
) %>%
  kable(
    caption = "<b style='display:block; text-align:center;'>Bảng 18: Kiểm định Chi-Squared cho Smoking x Diagnosis</b>",
    escape = FALSE, align = 'c' 
  ) %>%
  kable_styling(
    position = "center", font_size = 13, 
    full_width = FALSE, 
    bootstrap_options = c("striped", "hover") #striped table 
  )
Bảng 18: Kiểm định Chi-Squared cho Smoking x Diagnosis
Chi_Squared DF P_Value
X-squared 5.203 1 0.0225

Giá trị Chi-Squared là 5.203, với bậc tự do (DF) là 1 và p-value là 0.0225. Vì p-value (0.0225) nhỏ hơn mức ý nghĩa thông thường 0.05, có thể kết luận rằng có mối liên quan thống kê giữa việc hút thuốc (smoking) và khả năng bị bệnh tiểu đường (diagnosis).

Từ bảng tần số:

Trong số những người không hút thuốc (Smoking = 0), 61.63% không bị bệnh (Diagnosis = 0) và 38.37% bị bệnh (Diagnosis = 1).

Trong số những người hút thuốc (Smoking = 1), 55.77% không bị bệnh (Diagnosis = 0) và 44.23% bị bệnh (Diagnosis = 1).

Như vậy, tỷ lệ bị bệnh tiểu đường ở người hút thuốc (44.23%) cao hơn so với người không hút thuốc (38.37%), cho thấy hút thuốc có thể làm tăng nguy cơ mắc bệnh tiểu đường.

ggplot(table1, aes(x = Smoking, y = Diagnosis, fill = Frequency)) +
    geom_tile(color = "white") + 
    scale_fill_gradient(low = "white", high = "#0099CC") + 
    labs(
      title = "Heatmap tần suất chéo của Smoking x Diagnosis", 
      x = "Smoking", 
      y = "Diagnosis") +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
      axis.text.x = element_text(angle = 0, hjust = 1)
      )

5.2 Kiểm định giả thuyết về sự bằng nhau giữa tỷ lệ người hút thuốc bị bệnh và tỷ lệ người không hút thuốc bị bệnh.

Giả thuyết H0: Tỷ lệ người hút thuốc bị bệnh bằng với tỷ lệ người không hút thuốc bị bệnh.

smok1 <- data[data$Smoking == 1,]
smok0 <- data[data$Smoking == 0,]

smok1_benh <- smok1[smok1$Diagnosis == 1,]
smok0_benh <- smok0[smok0$Diagnosis == 1,]

vt1 <- c(nrow(smok1),nrow(smok0))
vt2 <- c(nrow(smok1_benh),nrow(smok0_benh))

prop.test(vt2, vt1)
## 
##  2-sample test for equality of proportions with continuity correction
## 
## data:  vt2 out of vt1
## X-squared = 5.2031, df = 1, p-value = 0.02255
## alternative hypothesis: two.sided
## 95 percent confidence interval:
##  0.007684176 0.109596507
## sample estimates:
##    prop 1    prop 2 
## 0.4423440 0.3837037

Vì p-value < 0.05, có bằng chứng thống kê để bác bỏ giả thuyết không (H0: hai tỷ lệ bằng nhau), cho thấy tỷ lệ bị bệnh giữa hai nhóm khác nhau ở mức ý nghĩa 5%.

5.3 Relative Risk

Mục đích: Đánh giá tỷ lệ mắc bệnh giữa người hút thuốc và người không hút thuốc.

library(epitools)
a1 <- table(data$Smoking, data$Diagnosis)
a1
##    
##       0   1
##   0 832 518
##   1 295 234
epitab(a1,method = 'riskratio') # Tính cho vị trí 22 và 12
## $tab
##    
##       0        p0   1        p1 riskratio   lower    upper    p.value
##   0 832 0.6162963 518 0.3837037  1.000000      NA       NA         NA
##   1 295 0.5576560 234 0.4423440  1.152827 1.02538 1.296116 0.02125206
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Tỷ số rủi ro (RR) là 1.1528, cho thấy nguy cơ bị bệnh tiểu đường ở người hút thuốc cao hơn 15.28% so với người không hút thuốc. Khoảng tin cậy 95% (1.0254 - 1.2961) không bao gồm 1, và p-value (0.0213) nhỏ hơn 0.05, nên mối quan hệ này có ý nghĩa thống kê.

Mục đích: Đánh giá tỷ lệ không mắc bệnh giữa người hút thuốc và người không hút thuốc.

epitab(a1, method = 'riskratio', rev = 'c')
## $tab
##    
##       1        p0   0        p1 riskratio     lower     upper    p.value
##   0 518 0.3837037 832 0.6162963 1.0000000        NA        NA         NA
##   1 234 0.4423440 295 0.5576560 0.9048504 0.8296334 0.9868869 0.02125206
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Risk Ratio (RR): 0.9049, so sánh tỷ lệ không bị bệnh giữa nhóm hút thuốc và không hút thuốc. Điều này có nghĩa là khả năng không bị bệnh ở nhóm hút thuốc thấp hơn khoảng 9.51% so với nhóm không hút thuốc. Khoảng tin cậy 95%: 0.8296 - 0.9869 (không bao gồm 1), cho thấy sự khác biệt có ý nghĩa thống kê. p-value: 0.0213 (Fisher exact), nhỏ hơn 0.05, xác nhận ý nghĩa thống kê.

5.4 Odds

Mục đích: Xác định cứ 1 người không bị bệnh thì có bao nhiêu người bị bệnh.

table1
##   Smoking Diagnosis Frequency Percent_Total Percent_of_Smoking
## 1       0         0       832         44.28              61.63
## 2       1         0       295         15.70              55.77
## 3       0         1       518         27.57              38.37
## 4       1         1       234         12.45              44.23
##   Percent_Diagnosis
## 1             73.82
## 2             26.18
## 3             68.88
## 4             31.12
odd1 <- table1[table1$Smoking == 1 & table1$Diagnosis == 1,"Frequency"]/
      table1[table1$Smoking == 1 & table1$Diagnosis == 0,"Frequency"]
cat("Trong nhóm người hút thuốc, cứ 1 người không bị bệnh tiểu đường thì có khoảng", 
    odd1, "người bị bệnh tiểu đường")
## Trong nhóm người hút thuốc, cứ 1 người không bị bệnh tiểu đường thì có khoảng 0.7932203 người bị bệnh tiểu đường
table1
##   Smoking Diagnosis Frequency Percent_Total Percent_of_Smoking
## 1       0         0       832         44.28              61.63
## 2       1         0       295         15.70              55.77
## 3       0         1       518         27.57              38.37
## 4       1         1       234         12.45              44.23
##   Percent_Diagnosis
## 1             73.82
## 2             26.18
## 3             68.88
## 4             31.12
odd2 <- table1[table1$Smoking == 0 & table1$Diagnosis == 1,"Frequency"]/
      table1[table1$Smoking == 0 & table1$Diagnosis == 0,"Frequency"]
cat("Trong nhóm người không hút thuốc, cứ 1 người không bị bệnh tiểu đường thì có khoảng", 
    odd2, "người bị bệnh tiểu đường")
## Trong nhóm người không hút thuốc, cứ 1 người không bị bệnh tiểu đường thì có khoảng 0.6225962 người bị bệnh tiểu đường

5.5 Odds Ratio

Mục đích: Đánh giá odds mắc bệnh trong nhóm người hút thuốc so với nhóm người không hút thuốc

cat("Odds bị bệnh (tỷ lệ người bị bệnh so với người không bị bệnh) trong nhóm người hút thuốc cao hơn", odd1/odd2, "lần so với nhóm người không hút thuốc.")
## Odds bị bệnh (tỷ lệ người bị bệnh so với người không bị bệnh) trong nhóm người hút thuốc cao hơn 1.274053 lần so với nhóm người không hút thuốc.
epitab(a1, method = 'oddsratio')
## $tab
##    
##       0        p0   1        p1 oddsratio    lower    upper    p.value
##   0 832 0.7382431 518 0.6888298  1.000000       NA       NA         NA
##   1 295 0.2617569 234 0.3111702  1.274053 1.039311 1.561815 0.02125206
## 
## $measure
## [1] "wald"
## 
## $conf.level
## [1] 0.95
## 
## $pvalue
## [1] "fisher.exact"

Điều này thể hiện rằng việc hút thuốc có thể làm tăng nguy cơ mắc bệnh tiểu đường.

5.6 Hồi quy Logistic

logit_model <- glm(Diagnosis ~ Smoking, data = data, family = "binomial")
summary(logit_model)
## 
## Call:
## glm(formula = Diagnosis ~ Smoking, family = "binomial", data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.47386    0.05597  -8.467   <2e-16 ***
## Smoking      0.24220    0.10390   2.331   0.0198 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2529.5  on 1878  degrees of freedom
## Residual deviance: 2524.1  on 1877  degrees of freedom
## AIC: 2528.1
## 
## Number of Fisher Scoring iterations: 4

Dựa trên kết quả hồi quy, ta thấy rằng log-odds(Diagnosis = 1|Smoking = 0) = -0.47386 => P(Diagnosis = 1|Smoking = 0) = 0.383. Vậy người không hút thuốc có xác suất mắc tiểu đường khoảng 38.3%.

Hệ số hồi quy của Smoking = 0.242 => Khi chuyển từ không hút sang hút thuốc, log-odds mắc tiểu đường tăng 0.242, khi đó Odds-ratio của nó sẽ bằng 1.274, thể hiện rằng người hút thuốc có nguy cơ mắc tiểu đường cao hơn khoảng 27.4% so với người không hút thuốc.

Để dễ dàng diễn giải, ta có thể tính toán hiệu ứng biên bằng cách lấy đạo hàm riêng của \(p\) theo biến độc lập \(X_i\) để đo lường sự thay đổi của xác suất \(p\) khi \(X_i\) thay đổi một đơn vị. Công thức tính hiệu ứng biên như sau:

\[ \frac{\partial p}{\partial X_i} = p \times (1 - p) \times \beta_i, \]

trong đó, \(p \times (1 - p)\) là xác suất biên, phụ thuộc vào giá trị cụ thể của \(X_i\). Ở nghiên cứu này, hiệu ứng biên sẽ được tính cho từng quan sát của \(X_i\) và rồi lấy trung bình các giá trị đó ((Average Marginal Effect - AME).

library(margins)
## Warning: package 'margins' was built under R version 4.4.3
summary(margins(logit_model))
##   factor    AME     SE      z      p  lower  upper
##  Smoking 0.0580 0.0247 2.3442 0.0191 0.0095 0.1064

Tiếp tục, hiệu ứng biên trung bình là 0.058, thể hiện rằng khi các yếu tố khác như nhau, người hút thuốc sẽ có trung bình xác suất mắc bệnh cao hơn 5.8 điểm phần trăm so với người không hút thuốc ở mức ý nghĩa thống kê 5%.

5.7 Hồi quy Probit

probit_model <- glm(Diagnosis ~ Smoking, data = data, family = binomial(link = "probit"))
summary(probit_model)
## 
## Call:
## glm(formula = Diagnosis ~ Smoking, family = binomial(link = "probit"), 
##     data = data)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -0.29577    0.03466  -8.534   <2e-16 ***
## Smoking      0.15074    0.06476   2.328   0.0199 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2529.5  on 1878  degrees of freedom
## Residual deviance: 2524.1  on 1877  degrees of freedom
## AIC: 2528.1
## 
## Number of Fisher Scoring iterations: 4

Kết quả ước lượng mô hình Probit cho thấy biến Smoking (hút thuốc) có ảnh hưởng có ý nghĩa thống kê đến xác suất được chẩn đoán mắc bệnh (Diagnosis). Hệ số hồi quy của biến Smoking là 0.15074, với giá trị p tương ứng 0.0199, nhỏ hơn ngưỡng 0.05, cho thấy mối quan hệ này có ý nghĩa thống kê ở mức 5%. Do hồi quy Probit sử dụng thang đo z (chuẩn hóa), hệ số này cho biết khi một người hút thuốc (biến Smoking = 1), thì điểm z (tức logit-probit score) tăng trung bình khoảng 0.15 đơn vị, từ đó làm tăng xác suất mắc bệnh. Ngược lại, hệ số chặn (intercept) có giá trị -0.29577 và có ý nghĩa thống kê rất cao (p < 2e-16), phản ánh xác suất cơ bản mắc bệnh ở nhóm không hút thuốc

Tương tự, hiệu ứng biên cho mô hình Probit được tính toán như sau:

summary(margins(probit_model))
##   factor    AME     SE      z      p  lower  upper
##  Smoking 0.0581 0.0249 2.3382 0.0194 0.0094 0.1068

Kết quả tính hiệu ứng biên trung bình (AME) từ mô hình Probit cho thấy biến Smoking (hút thuốc) có ảnh hưởng dương và có ý nghĩa thống kê đến xác suất được chẩn đoán mắc bệnh (Diagnosis). Cụ thể, AME của biến Smoking là 0.0581, nghĩa là khi một cá nhân hút thuốc, thì xác suất được chẩn đoán mắc bệnh tăng trung bình khoảng 5.81 điểm phần trăm, giữ các yếu tố khác không đổi.

6 Tác động của tình trạng kinh tế xã hội đến việc bị bệnh tiểu đường

Tình trạng kinh tế xã hội là một biến độc lập có 3 nhãn với các giá trị:

0 = Thấp,
1 = Trung bình,
2 = Cao.

Khi thực hiện hồi quy với biến độc lập có nhiều hơn 2 nhãn, R sẽ tự động chuyển nó thành các biến giả. Và kết quả hồi quy sẽ so sánh từng nhóm với nhóm tham chiếu (reference group). Nhóm tham chiếu sẽ là nhóm đầu tiên theo thứ tự alphabet.

6.1 Hồi quy logistic

# Chuyển thành định dạng factor
data$SocioeconomicStatus <- factor(data$SocioeconomicStatus)
logit_model <- glm(Diagnosis ~ SocioeconomicStatus, data = data, family = "binomial")
summary(logit_model)
## 
## Call:
## glm(formula = Diagnosis ~ SocioeconomicStatus, family = "binomial", 
##     data = data)
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -0.50222    0.08743  -5.744 9.23e-09 ***
## SocioeconomicStatus1  0.13931    0.11377   1.225    0.221    
## SocioeconomicStatus2  0.13657    0.12359   1.105    0.269    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2529.5  on 1878  degrees of freedom
## Residual deviance: 2527.7  on 1876  degrees of freedom
## AIC: 2533.7
## 
## Number of Fisher Scoring iterations: 4

Mô hình hồi quy logistic được xây dựng để phân tích mối quan hệ giữa tình trạng chẩn đoán bệnh (Diagnosis) và mức độ kinh tế xã hội (SocioeconomicStatus). Trong đó, SocioeconomicStatus là một biến định tính với ba nhóm, và nhóm tham chiếu là mức 0. Kết quả cho thấy hệ số chặn (intercept) là -0.50222 với giá trị p rất nhỏ (p < 0.001), cho thấy xác suất cơ bản mắc bệnh ở nhóm SocioeconomicStatus = 0 là có ý nghĩa thống kê. Hai hệ số của SocioeconomicStatus1 (0.13931) và SocioeconomicStatus2 (0.13657) lần lượt thể hiện sự thay đổi log-odds mắc bệnh so với nhóm có tình trạng kinh tế thấp nhất, tuy nhiên cả hai đều có giá trị p lớn hơn 0.05 (lần lượt là 0.221 và 0.269), cho thấy không có sự khác biệt có ý nghĩa thống kê giữa các mức độ kinh tế xã hội đối với xác suất mắc bệnh trong mô hình này. Nói cách khác, theo mô hình này, mức độ kinh tế xã hội không phải là yếu tố dự báo đáng tin cậy cho khả năng mắc bệnh.