Buoi 3 - Ngay 2

Nap các package

library(readxl)
library(dplyr)

## 
## 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(tidyr)
library(ggplot2)
library(psych)

## 
## Attaching package: 'psych'

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

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:psych':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

library(gtsummary)
library(gt)
library(scales)

## 
## Attaching package: 'scales'

## The following objects are masked from 'package:psych':
## 
##     alpha, rescale

library(binom)
library(purrr)

## 
## Attaching package: 'purrr'

## The following object is masked from 'package:scales':
## 
##     discard

## The following object is masked from 'package:car':
## 
##     some

Nap Data

library(readxl)
data <- read_excel("FETP_Integrated_Surveillance_2026.xlsx", sheet = "Line_list")

Câu 1. Kiểm tra phân phối của biến bmi bằng histogram, Q-Q plot, Shapiro-Wilk, skewness và kurtosis. Ghi kết quả và kết luận.

hist(data$bmi)

qqnorm(data$bmi)
qqline(data$bmi)

qqnorm(data$bmi, main = "Normal Q-Q plot of age")
qqline(data$bmi, col = "red", lwd = 2)

# Shapiro-Wilk
shapiro.test(data$bmi)

## 
##  Shapiro-Wilk normality test
## 
## data:  data$bmi
## W = 0.83776, p-value < 2.2e-16

# Kolmogorov-Smirnov
ks.test(data$bmi, "pnorm", mean(data$bmi), sd(data$bmi))

## Warning in ks.test.default(data$bmi, "pnorm", mean(data$bmi), sd(data$bmi)):
## ties should not be present for the one-sample Kolmogorov-Smirnov test

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  data$bmi
## D = 0.16861, p-value = 9.002e-13
## alternative hypothesis: two-sided

# Skewness và kurtosis

psych::describe(data$bmi)

Kết quả:

Histogram: Biểu đồ không có dạng hình chuông cân đối; dữ liệu tập trung nhiều ở khoảng giá trị nhỏ (< 30) và kéo đuôi dài về phía bên phải (có nhiều giá trị cao).

Q-Q plot: Các điểm dữ liệu lệch xa khỏi đường thẳng tham chiếu ở cả hai đầu (đặc biệt là đuôi phải).

Shapiro - Wilk: W = 0.83776, p-value < 2.2e-16. Dữ liệu bmi hiện không phân phối chuẩn.

Kolmogorov-Smirnov test: p-value = 9.002e-13. Dữ liệu bmi hiện không phân phối chuẩn.

Skewness = 1.66 –> Độ lệch dương tương đối cao, khẳng định phân phối bị lệch phải rõ rệt.

kurtosis =

Kết luận chung: Biến bmi không tuân theo quy luật phân phối chuẩn.

Câu 2. Tính các chỉ số mô tả của bmi và điền vào bảng. Tính thêm khoảng tin cậy 95% của trung bình bmi. Gợi ý: summarise(), quantile(), IQR(), t.test()

# Tính Khoảng tin cậy 95% của Trung bình bằng hàm t.test()
ttest_res <- t.test(data$bmi)
ci_lower <- ttest_res$conf.int[1]
ci_upper <- ttest_res$conf.int[2]

# Tính các chỉ số thống kê mô tả bằng summarise() và các hàm theo gợi ý
bmi_summary <- data %>%
  summarise(
    n = sum(!is.na(bmi)),
    Trung_binh = mean(bmi, na.rm = TRUE),
    Do_lech_chuan = sd(bmi, na.rm = TRUE),
    Trung_vi = median(bmi, na.rm = TRUE),
    Q1 = quantile(bmi, 0.25, na.rm = TRUE),
    Q3 = quantile(bmi, 0.75, na.rm = TRUE),
    Nho_nhat = min(bmi, na.rm = TRUE),
    Lon_nhat = max(bmi, na.rm = TRUE)
  )

# In các chỉ số thống kê mô tả
print(bmi_summary)

## # A tibble: 1 × 8
##       n Trung_binh Do_lech_chuan Trung_vi    Q1    Q3 Nho_nhat Lon_nhat
##   <int>      <dbl>         <dbl>    <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1   500       37.6          30.8     28.0  16.5  48.5      2.5     172.

# In Khoảng tin cậy 95% của Trung bình
cat(sprintf("KTC 95%% cua trung binh: [%.2f - %.2f]\n", ci_lower, ci_upper))

## KTC 95% cua trung binh: [34.91 - 40.31]

Câu 3.Tạo bảng mô tả các biến age, bmi, temperature, pulse, sbp, dbp, platelet, hematocrit và wbc. Gợi ý: across(), pivot_longer(), sprintf()

# Danh sách các biến cần mô tả
vars_list <- c("age", "bmi", "temperature", "pulse", "sbp", "dbp", 
               "platelet", "hematocrit", "wbc")

# Tạo bảng mô tả các biến bằng across(), pivot_longer() và sprintf()
# 4.5. Mô tả nhiều biến bằng across()
continuous_summary <- data %>%
  summarise(
    across(
      all_of(vars_list),
      list(
        n = ~ sum(!is.na(.x)),
        mean = ~ mean(.x, na.rm = TRUE),
        sd = ~ sd(.x, na.rm = TRUE),
        median = ~ median(.x, na.rm = TRUE),
        q1 = ~ quantile(.x, 0.25, na.rm = TRUE),
        q3 = ~ quantile(.x, 0.75, na.rm = TRUE),
        min = ~ min(.x, na.rm = TRUE),
        max = ~ max(.x, na.rm = TRUE)
      ),
      .names = "{.col}_{.fn}"
    )
  )
  continuous_summary

# Chuyển bảng rộng thành bảng dài, mỗi biến một dòng
continuous_summary_long <- continuous_summary %>%
  pivot_longer(
    cols = everything(),
    names_to = c("variable", ".value"),
    names_pattern = "(.*)_(n|mean|sd|median|q1|q3|min|max)$"
  ) %>%
  mutate(
    `Mean ± SD` = sprintf("%.1f ± %.1f", mean, sd),
    `Median (Q1–Q3)` = sprintf("%.1f (%.1f–%.1f)", median, q1, q3),
    `Min–Max` = sprintf("%.1f–%.1f", min, max)
  ) %>%
  select(variable, n, `Mean ± SD`, `Median (Q1–Q3)`, `Min–Max`)
continuous_summary_long

# Trình bày bảng chuyên nghiệp bằng gt
continuous_summary_gt <- continuous_summary_long %>%
  gt() %>%
  tab_header(title = "Descriptive statistics of continuous variables") %>%
  cols_label(
    variable = "Variable",
    n = "n"
  )
continuous_summary_gt

Descriptive statistics of continuous variables
Variable	n	Mean ± SD	Median (Q1–Q3)	Min–Max
age	500	46.4 ± 26.1	47.0 (25.0–69.0)	0.0–90.0
bmi	500	37.6 ± 30.8	28.0 (16.5–48.5)	2.5–172.5
temperature	500	38.3 ± 1.1	38.3 (37.3–39.2)	36.5–40.0
pulse	500	94.3 ± 20.1	95.0 (76.0–112.0)	60.0–130.0
sbp	500	135.6 ± 25.8	136.0 (113.0–155.0)	90.0–180.0
dbp	500	85.6 ± 14.9	86.0 (73.0–98.0)	60.0–110.0
platelet	500	161792.6 ± 81905.0	164852.5 (93395.5–230052.0)	20180.0–299968.0
hematocrit	500	42.1 ± 7.5	42.5 (34.9–48.9)	30.0–55.0
wbc	500	9.9 ± 4.6	9.8 (5.7–14.0)	2.1–17.9

Câu 4. Lập bảng phân bố disease (diabetes, hypertension, ckd, heart_disease) gồm n và %. Sắp xếp từ nhóm có tần số cao nhất đến thấp nhất.

# Tổng số quan sát (N)
total_n <- nrow(data)

# Thống kê từng bệnh riêng lẻ (sắp xếp giảm dần theo tần số n)
disease_cols <- c("diabetes", "hypertension", "ckd", "heart_disease")

individual_diseases <- data %>%
  summarise(across(all_of(disease_cols), ~ sum(. == "Yes", na.rm = TRUE))) %>%
  pivot_longer(cols = everything(), names_to = "Bệnh", values_to = "n") %>%
  mutate(`%` = (n / total_n) * 100) %>%
  arrange(desc(n))  # Sắp xếp từ cao nhất đến thấp nhất

print("--- Bảng từng bệnh (Sắp xếp giảm dần) ---")

## [1] "--- Bảng từng bệnh (Sắp xếp giảm dần) ---"

print(individual_diseases)

## # A tibble: 4 × 3
##   Bệnh              n   `%`
##   <chr>         <int> <dbl>
## 1 hypertension    273  54.6
## 2 heart_disease   271  54.2
## 3 ckd             248  49.6
## 4 diabetes        243  48.6

# Tính số lượng bệnh (# of diseases) của mỗi người
data <- data %>%
  mutate(
    num_diseases = rowSums(select(., all_of(disease_cols)) == "Yes", na.rm = TRUE),
    disease_cat = case_when(
      num_diseases == 0 ~ "0",
      num_diseases == 1 ~ "1",
      num_diseases >= 2 ~ ">= 2"
    )
  )

num_diseases_summary <- data %>%
  count(disease_cat, name = "n") %>%
  mutate(`%` = (n / total_n) * 100)

print("--- Bảng số lượng bệnh (# of diseases) ---")

## [1] "--- Bảng số lượng bệnh (# of diseases) ---"

print(num_diseases_summary)

## # A tibble: 3 × 3
##   disease_cat     n   `%`
##   <chr>       <int> <dbl>
## 1 0              25   5  
## 2 1             109  21.8
## 3 >= 2          366  73.2

Câu 5. Mô tả BMI theo từng bệnh.

data <- read_excel("FETP_Integrated_Surveillance_2026.xlsx", sheet = "Line_list")
# Danh sách 4 bệnh cần phân tích
disease_cols <- c("diabetes", "hypertension", "ckd", "heart_disease")

# Tính chỉ số BMI cho từng bệnh đơn lẻ (nhóm mắc bệnh = "Yes")
bmi_by_disease <- map_dfr(disease_cols, function(col) {
  data %>%
    filter(.data[[col]] == "Yes") %>%
    summarise(
      Bệnh = col,
      n = sum(!is.na(bmi)),
      `Mean ± SD` = sprintf("%.2f ± %.2f", mean(bmi, na.rm = TRUE), sd(bmi, na.rm = TRUE)),
      `Median (Q1–Q3)` = sprintf("%.2f (%.2f–%.2f)", 
                                 median(bmi, na.rm = TRUE), 
                                 quantile(bmi, 0.25, na.rm = TRUE), 
                                 quantile(bmi, 0.75, na.rm = TRUE))
    )
})

# 4. Tính chỉ số BMI theo số lượng bệnh mắc phải (# of diseases)
bmi_by_num_diseases <- data %>%
  mutate(
    # Tính số bệnh mắc phải của từng người
    num_diseases_raw = rowSums(select(., all_of(disease_cols)) == "Yes", na.rm = TRUE),
    disease_group = case_when(
      num_diseases_raw == 0 ~ "0",
      num_diseases_raw == 1 ~ "1",
      num_diseases_raw >= 2 ~ ">= 2"
    )
  ) %>%
  group_by(disease_group) %>%
  summarise(
    n = sum(!is.na(bmi)),
    `Mean ± SD` = sprintf("%.2f ± %.2f", mean(bmi, na.rm = TRUE), sd(bmi, na.rm = TRUE)),
    `Median (Q1–Q3)` = sprintf("%.2f (%.2f–%.2f)", 
                               median(bmi, na.rm = TRUE), 
                               quantile(bmi, 0.25, na.rm = TRUE), 
                               quantile(bmi, 0.75, na.rm = TRUE))
  ) %>%
  rename(Bệnh = disease_group)

# In kết quả ra màn hình
print("--- BMI THEO TỪNG BỆNH ---")

## [1] "--- BMI THEO TỪNG BỆNH ---"

print(bmi_by_disease)

## # A tibble: 4 × 4
##   Bệnh              n `Mean ± SD`   `Median (Q1–Q3)`   
##   <chr>         <int> <chr>         <chr>              
## 1 diabetes        243 37.96 ± 29.08 28.90 (17.70–51.30)
## 2 hypertension    273 36.43 ± 30.68 27.90 (16.20–46.70)
## 3 ckd             248 38.28 ± 32.37 28.20 (15.67–53.45)
## 4 heart_disease   271 37.24 ± 31.97 27.70 (16.10–47.15)

print("--- BMI THEO SỐ LƯỢNG BỆNH ---")

## [1] "--- BMI THEO SỐ LƯỢNG BỆNH ---"

print(bmi_by_num_diseases)

## # A tibble: 3 × 4
##   Bệnh      n `Mean ± SD`   `Median (Q1–Q3)`   
##   <chr> <int> <chr>         <chr>              
## 1 0        25 44.78 ± 32.16 38.40 (16.20–61.30)
## 2 1       109 35.98 ± 27.24 27.40 (18.80–42.60)
## 3 >= 2    366 37.61 ± 31.64 27.90 (16.22–48.58)

Câu 6. Lập bảng tỷ lệ nhập viện trong từng nhóm bệnh. Tỷ lệ phần trăm phải được tính theo hàng.

# Danh sách 4 bệnh cần phân tích
disease_cols <- c("diabetes", "hypertension", "ckd", "heart_disease")

# Tính tỷ lệ nhập viện cho từng bệnh đơn lẻ (nhóm mắc bệnh = "Yes")
hosp_by_disease <- map_dfr(disease_cols, function(col) {
  data %>%
    filter(.data[[col]] == "Yes") %>%
    summarise(
      Bệnh = col,
      no_n = sum(hospitalized == "No", na.rm = TRUE),
      yes_n = sum(hospitalized == "Yes", na.rm = TRUE),
      Tổng = n()
    ) %>%
    mutate(
      `Không nhập viện n (%)` = sprintf("%d (%.1f%%)", no_n, (no_n / Tổng) * 100),
      `Có nhập viện n (%)` = sprintf("%d (%.1f%%)", yes_n, (yes_n / Tổng) * 100)
    ) %>%
    select(Bệnh, `Không nhập viện n (%)`, `Có nhập viện n (%)`, Tổng)
})

# Tính tỷ lệ nhập viện theo số lượng bệnh mắc phải (# of diseases)
hosp_by_num_diseases <- data %>%
  mutate(
    num_diseases_raw = rowSums(select(., all_of(disease_cols)) == "Yes", na.rm = TRUE),
    disease_group = case_when(
      num_diseases_raw == 0 ~ "0",
      num_diseases_raw == 1 ~ "1",
      num_diseases_raw >= 2 ~ ">= 2"
    )
  ) %>%
  group_by(disease_group) %>%
  summarise(
    no_n = sum(hospitalized == "No", na.rm = TRUE),
    yes_n = sum(hospitalized == "Yes", na.rm = TRUE),
    Tổng = n()
  ) %>%
  mutate(
    `Không nhập viện n (%)` = sprintf("%d (%.1f%%)", no_n, (no_n / Tổng) * 100),
    `Có nhập viện n (%)` = sprintf("%d (%.1f%%)", yes_n, (yes_n / Tổng) * 100)
  ) %>%
  rename(Bệnh = disease_group) %>%
  select(Bệnh, `Không nhập viện n (%)`, `Có nhập viện n (%)`, Tổng)

# In kết quả ra màn hình Console
print("--- TỶ LỆ NHẬP VIỆN THEO TỪNG BỆNH ---")

## [1] "--- TỶ LỆ NHẬP VIỆN THEO TỪNG BỆNH ---"

print(hosp_by_disease)

## # A tibble: 4 × 4
##   Bệnh          `Không nhập viện n (%)` `Có nhập viện n (%)`  Tổng
##   <chr>         <chr>                   <chr>                <int>
## 1 diabetes      127 (52.3%)             116 (47.7%)            243
## 2 hypertension  133 (48.7%)             140 (51.3%)            273
## 3 ckd           127 (51.2%)             121 (48.8%)            248
## 4 heart_disease 126 (46.5%)             145 (53.5%)            271

print("--- TỶ LỆ NHẬP VIỆN THEO SỐ LƯỢNG BỆNH ---")

## [1] "--- TỶ LỆ NHẬP VIỆN THEO SỐ LƯỢNG BỆNH ---"

print(hosp_by_num_diseases)

## # A tibble: 3 × 4
##   Bệnh  `Không nhập viện n (%)` `Có nhập viện n (%)`  Tổng
##   <chr> <chr>                   <chr>                <int>
## 1 0     14 (56.0%)              11 (44.0%)              25
## 2 1     64 (58.7%)              45 (41.3%)             109
## 3 >= 2  178 (48.6%)             188 (51.4%)            366

Câu 7.Vẽ biểu đồ cột phân bố giới tính bằng ggplot2, có nhãn số lượng trên mỗi cột.

data <- data %>%
  mutate(
    sex_vn = case_when(
      sex == "Female" ~ "Nữ",
      sex == "Male"   ~ "Nam",
      TRUE            ~ sex
    ),
  
    sex_vn = factor(sex_vn, levels = c("Nữ", "Nam"))
  )

# 4. Vẽ biểu đồ cột bằng ggplot2
ggplot(data, aes(x = sex_vn, fill = sex_vn)) +
  # Tạo cột
  geom_bar(width = 0.5, color = "black") +
  # Thêm nhãn số lượng trên đỉnh mỗi cột
  geom_text(stat = "count", aes(label = after_stat(count)), 
            vjust = -0.5, size = 5, fontface = "bold") +
  # Tùy chỉnh thang đo trục Y (giúp nhãn số lượng không bị che mất)
  scale_y_continuous(limits = c(0, 300)) +
  # Đổi màu cột phù hợp cho "Nữ" và "Nam"
  scale_fill_manual(values = c("Nữ" = "#E75480", "Nam" = "#4682B4")) +
  # Đặt tiêu đề và tên các trục bằng tiếng Việt
  labs(
    title = "Biểu đồ phân bố giới tính",
    x = "Giới tính",
    y = "Số lượng (n)"
  ) +
  # Tùy chỉnh giao diện (Theme)
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), # Căn giữa tiêu đề
    legend.position = "none" # Ẩn phần chú thích vì tên trục X đã rõ ràng
  )

Câu 8. Vẽ boxplot của age theo disease.

# Vẽ Boxplot bằng ggplot2
ggplot(data, aes(x = disease, y = age, fill = disease)) +
  # Vẽ biểu đồ hộp
  geom_boxplot(alpha = 0.7, outlier.color = "red", outlier.shape = 16, outlier.size = 2) +
  # Bổ sung điểm trung bình (Mean) dạng điểm đỏ nổi bật
  stat_summary(fun = mean, geom = "point", shape = 18, size = 3, color = "darkred") +
  # Tùy chỉnh màu sắc sinh động cho từng loại bệnh
  scale_fill_brewer(palette = "Set2") +
  # Đặt tiêu đề và tên các trục bằng tiếng Việt
  labs(
    title = "Biểu đồ Boxplot phân bố Tuổi theo Loại bệnh (Disease)",
    x = "Loại bệnh",
    y = "Tuổi (Age)"
  ) +
  # Tùy chỉnh giao diện (Theme) đẹp mắt
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"), # Căn giữa tiêu đề
    axis.text.x = element_text(face = "bold"),
    legend.position = "none" # Ẩn chú thích vì màu đã khớp với trục X
  )