#Việc 1. Phân tích mô tả

#1.1 Đọc dữ liệu "Obesity data.csv" vào R và gọi dữ liệu là "ob" 
library(readr)
library(gtsummary)
Obesity_data <- read_csv("D:/2025-10-Tap huan_Phan  tich du lieu va cong bo Quoc te danh cho giang vien/Obesity data.csv")
## Rows: 1217 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): gender
## dbl (12): id, height, weight, bmi, age, WBBMC, wbbmd, fat, lean, pcfat, hype...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
ob=Obesity_data

#1.2 Mô tả đặc điểm tuổi (age), giới tính  (gender), cân nặng (weight), chiều cao (height), tỉ trọng mỡ  (pcfat), tiền căn bệnh cao huyết áp (hypertension)và tiền căn bệnh tiểu đường  (diabetes)
#Bước 1: Mô tả biến liên tục (tuổi, cân nặng, chiều cao, pcfat)
summary(ob[, c("age", "weight", "height", "pcfat")])
##       age            weight          height          pcfat     
##  Min.   :13.00   Min.   :34.00   Min.   :136.0   Min.   : 9.2  
##  1st Qu.:35.00   1st Qu.:49.00   1st Qu.:151.0   1st Qu.:27.0  
##  Median :48.00   Median :54.00   Median :155.0   Median :32.4  
##  Mean   :47.15   Mean   :55.14   Mean   :156.7   Mean   :31.6  
##  3rd Qu.:58.00   3rd Qu.:61.00   3rd Qu.:162.0   3rd Qu.:36.8  
##  Max.   :88.00   Max.   :95.00   Max.   :185.0   Max.   :48.4
#Bước 2: Mô tả biến phân loại (giới tính, cao huyết áp, tiểu đường)
table(ob$gender)
## 
##   F   M 
## 862 355
prop.table(table(ob$gender))  # tỉ lệ %
## 
##         F         M 
## 0.7082991 0.2917009
table(ob$hypertension)
## 
##   0   1 
## 600 617
prop.table(table(ob$hypertension))
## 
##         0         1 
## 0.4930156 0.5069844
table(ob$diabetes)
## 
##    0    1 
## 1082  135
prop.table(table(ob$diabetes))
## 
##         0         1 
## 0.8890715 0.1109285
#1.3 Bạn nhận xét như thế nào về kết quả của tiền căn bệnh cao huyết áp và tiểu đường. Làm cách nào để trình bày kết quả tốt hơn?
 #Tỉ lệ người có tiền căn cao huyết áp chiếm hơn một nửa (50.7%), cho thấy đây là một bệnh lý phổ biến trong mẫu nghiên cứu này.
#Chỉ khoảng 11% người tham gia có tiền căn tiểu đường, tức là ít phổ biến hơn so với cao huyết áp trong cùng mẫu dữ liệu.

# Tạo bảng tóm tắt phần trăm
library(ggplot2)

# Dữ liệu hypertension
hypertension_df <- as.data.frame(prop.table(table(ob$hypertension)))
colnames(hypertension_df) <- c("Status", "Proportion")
hypertension_df$Status <- factor(hypertension_df$Status, labels = c("No", "Yes"))

# Dữ liệu diabetes
diabetes_df <- as.data.frame(prop.table(table(ob$diabetes)))
colnames(diabetes_df) <- c("Status", "Proportion")
diabetes_df$Status <- factor(diabetes_df$Status, labels = c("No", "Yes"))

# Vẽ biểu đồ
ggplot(hypertension_df, aes(x = Status, y = Proportion, fill = Status)) +
  geom_bar(stat = "identity") +
  ggtitle("Tỉ lệ tiền căn cao huyết áp") +
  ylab("Tỉ lệ") + xlab("Trạng thái") +
  scale_fill_manual(values = c("#66c2a5", "#fc8d62")) +
  theme_minimal()

#1.4 Bạn muốn trình bày kết quả trung vị (Q1, Q3) thay vi trung vị (min, max) cho biến liên tục.
#Tạo bảng tóm tắt trung vị (Q1–Q3) cho nhiều biến:
vars <- c("age", "weight", "height", "pcfat")

summary_stats <- sapply(ob[vars], function(x) {
  q1 <- quantile(x, 0.25, na.rm = TRUE)
  med <- median(x, na.rm = TRUE)
  q3 <- quantile(x, 0.75, na.rm = TRUE)
  sprintf("%.1f (%.1f – %.1f)", med, q1, q3)
})

summary_stats
##                     age                  weight                  height 
##    "48.0 (35.0 – 58.0)"    "54.0 (49.0 – 61.0)" "155.0 (151.0 – 162.0)" 
##                   pcfat 
##    "32.4 (27.0 – 36.8)"
#1.5 Mô tả đặc điểm tuổi (age), cân nặng (weight), chiều cao (height), tỉ trọng mỡ  (pcfat) và tiền căn bệnh cao huyết áp (hypertension) theo giới tính (gender)
#Dữ liệu liên tục (age, weight, height, pcfat) theo giới tính
# Hiển thị trung vị (Q1 – Q3) theo giới tính
vars <- c("age", "weight", "height", "pcfat")

# Hàm mô tả theo gender
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
summary_by_gender <- ob %>%
  group_by(gender) %>%
  summarise(across(all_of(vars), 
                   ~ sprintf("%.1f (%.1f – %.1f)", 
                             median(., na.rm = TRUE), 
                             quantile(., 0.25, na.rm = TRUE), 
                             quantile(., 0.75, na.rm = TRUE))))

summary_by_gender
## # A tibble: 2 × 5
##   gender age                weight             height                pcfat      
##   <chr>  <chr>              <chr>              <chr>                 <chr>      
## 1 F      49.0 (39.0 – 59.0) 51.0 (47.0 – 57.0) 153.0 (150.0 – 157.0) 34.7 (31.5…
## 2 M      44.0 (24.0 – 56.0) 62.0 (55.0 – 68.0) 165.0 (160.0 – 169.0) 24.6 (20.4…
#Biến phân loại (hypertension) theo giới tính
#Dùng table() để đếm số lượng và tỷ lệ:
    # Đếm số lượng
  table(ob$gender, ob$hypertension)
##    
##       0   1
##   F 430 432
##   M 170 185
# Tính tỷ lệ theo hàng (giới tính)
prop.table(table(ob$gender, ob$hypertension), margin = 1)
##    
##             0         1
##   F 0.4988399 0.5011601
##   M 0.4788732 0.5211268
# Tính tỷ lệ theo hàng (giới tính)
prop.table(table(ob$gender, ob$hypertension), margin = 1)
##    
##             0         1
##   F 0.4988399 0.5011601
##   M 0.4788732 0.5211268
#1.6 Đánh giá xem đặc điểm nào là khác biệt đáng kể (significant difference) giữa nam và nữ.
#Kiểm định T-test cho biến liên tục
# Kiểm định t-test giữa nam và nữ cho các biến liên tục
t.test(age ~ gender, data = ob)
## 
##  Welch Two Sample t-test
## 
## data:  age by gender
## t = 4.2608, df = 587.49, p-value = 2.372e-05
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  2.626083 7.117272
## sample estimates:
## mean in group F mean in group M 
##        48.57309        43.70141
t.test(weight ~ gender, data = ob)
## 
##  Welch Two Sample t-test
## 
## data:  weight by gender
## t = -16.952, df = 551.85, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -10.836942  -8.586319
## sample estimates:
## mean in group F mean in group M 
##        52.31090        62.02254
t.test(height ~ gender, data = ob)
## 
##  Welch Two Sample t-test
## 
## data:  height by gender
## t = -29.125, df = 562.31, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -12.56161 -10.97433
## sample estimates:
## mean in group F mean in group M 
##        153.2912        165.0592
t.test(pcfat ~ gender, data = ob)
## 
##  Welch Two Sample t-test
## 
## data:  pcfat by gender
## t = 29.768, df = 602.01, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##   9.822548 11.210140
## sample estimates:
## mean in group F mean in group M 
##        34.67241        24.15607
#Kiểm định Chi-squared cho biến phân loại
# Tạo bảng tần số giữa gender và hypertension
table_gender_htn <- table(ob$gender, ob$hypertension)

# Thực hiện kiểm định chi-squared
chisq.test(table_gender_htn)
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table_gender_htn
## X-squared = 0.32515, df = 1, p-value = 0.5685
#Tạo bảng tổng hợp
library(gtsummary)

# Tạo bảng thống kê và kiểm định
summary_table <- ob %>%
  select(gender, age, weight, height, pcfat, hypertension) %>%
  tbl_summary(by = gender, # so sánh theo giới tính
              statistic = list(all_continuous() ~ "{mean} ± {sd}"),
              digits = all_continuous() ~ 1,
              label = list(
                age ~ "Tuổi",
                weight ~ "Cân nặng (kg)",
                height ~ "Chiều cao (cm)",
                pcfat ~ "Tỉ lệ mỡ (%)",
                hypertension ~ "Cao huyết áp"
              )) %>%
   add_p(test = list(
    all_continuous() ~ "t.test",
    all_categorical() ~ "chisq.test"
  )) %>%
  bold_p(t = 0.05)



# Hiển thị bảng
summary_table
Characteristic F
N = 8621		 M
N = 3551		 p-value2
Tuổi 48.6 ± 16.4 43.7 ± 18.8 <0.001
Cân nặng (kg) 52.3 ± 7.7 62.0 ± 9.6 <0.001
Chiều cao (cm) 153.3 ± 5.6 165.1 ± 6.7 <0.001
Tỉ lệ mỡ (%) 34.7 ± 5.2 24.2 ± 5.8 <0.001
Cao huyết áp 432 (50%) 185 (52%) 0.6
1 Mean ± SD; n (%)
2 Welch Two Sample t-test; Pearson’s Chi-squared test 
#Ketluan: Kết luận về sự khác biệt giữa nam và nữ

  #Khi so sánh đặc điểm giữa nam (n = 355) và nữ (n = 862), kết quả cho thấy:
  
  #Tuổi: Trung bình nữ lớn tuổi hơn nam (48.6 ± 16.4 vs. 43.7 ± 18.8, p < 0.001) → khác biệt có ý nghĩa thống kê.

  #Cân nặng: Nam có cân nặng trung bình cao hơn nữ (62.0 ± 9.6 kg vs. 52.3 ± 7.7 kg, p < 0.001) → khác biệt có ý nghĩa thống kê.

#Chiều cao: Nam cao hơn nữ trung bình khoảng 12 cm (165.1 ± 6.7 cm vs. 153.3 ± 5.6 cm, p < 0.001) → khác biệt có ý nghĩa thống kê.

#Tỉ lệ mỡ cơ thể: Nữ có tỉ lệ mỡ cao hơn nam (34.7 ± 5.2% vs. 24.2 ± 5.8%, p < 0.001) → khác biệt có ý nghĩa thống kê.

#Tiền căn cao huyết áp: Tỷ lệ mắc cao huyết áp ở nam và nữ tương đương (52% vs. 50%, p = 0.6) → không có khác biệt có ý nghĩa thống kê.

#Việc 2. Phân tích khác biệt giữa 2 nhóm: 
  #2.1 Giả sử bạn có dữ liệu về tải trọng của 2 nhóm A và B như sau:
  #Nhóm A (n= 7): 14, 4, 10, 6, 3, 11, 12
  #Nhóm B (n= 9): 16, 17, 13, 12, 7, 16, 11, 8, 7
  #Nhập dữ liệu trên vào R.
#2.1 Nhập dữ liệu nhóm A và B
# Nhập dữ liệu nhóm A và nhóm B
group_A <- c(14, 4, 10, 6, 3, 11, 12)
group_B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)
#2.2 Tải trọng có tuân theo phân bổ chuẩn (normal distribution) không?
#Bước 1: Kiểm định Shapiro-Wilk
# Kiểm tra nhóm A
shapiro.test(group_A)
## 
##  Shapiro-Wilk normality test
## 
## data:  group_A
## W = 0.92541, p-value = 0.5126
# Kiểm tra nhóm B
shapiro.test(group_B)
## 
##  Shapiro-Wilk normality test
## 
## data:  group_B
## W = 0.89641, p-value = 0.2319
#Bước 2: Vẽ biểu đồ trực quan
# Histogram
hist(group_A, main = "Histogram - Group A", xlab = "Tải trọng", col = "lightblue")

hist(group_B, main = "Histogram - Group B", xlab = "Tải trọng", col = "lightgreen")

# Q-Q plot
qqnorm(group_A); qqline(group_A, col = "red", lwd = 2)

qqnorm(group_B); qqline(group_B, col = "blue", lwd = 2)

#Kết luận:  Vì cả p-value của nhóm A và nhóm B đều > 0.05, ta không có đủ bằng chứng để bác bỏ giả thuyết rằng dữ liệu có phân phối chuẩn.
#Kết luận: Dữ liệu tải trọng của cả nhóm A và nhóm B tuân theo phân phối chuẩn.
#2.3 Mô tả đặc điểm về tải trọng giữa 2 nhóm
summary(group_A)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   3.000   5.000  10.000   8.571  11.500  14.000
summary(group_B)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    7.00    8.00   12.00   11.89   16.00   17.00
mean(group_A); sd(group_A)
## [1] 8.571429
## [1] 4.237025
mean(group_B); sd(group_B)
## [1] 11.88889
## [1] 3.95109
# Tạo bảng mô tả 
data.frame(
  Group = c("A", "B"),
  N = c(length(group_A), length(group_B)),
  Mean = c(mean(group_A), mean(group_B)),
  SD = c(sd(group_A), sd(group_B)),
  Median = c(median(group_A), median(group_B)),
  Min = c(min(group_A), min(group_B)),
  Max = c(max(group_A), max(group_B))
)
##   Group N      Mean       SD Median Min Max
## 1     A 7  8.571429 4.237025     10   3  14
## 2     B 9 11.888889 3.951090     12   7  17
#2.4 Thực hiện phép kiểm t để đánh giá khác biệt về tải trọng của 2 nhóm. Bạn nhận xét gì về kết quả này.
t.test(group_A, group_B)
## 
##  Welch Two Sample t-test
## 
## data:  group_A and group_B
## t = -1.6, df = 12.554, p-value = 0.1345
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.813114  1.178194
## sample estimates:
## mean of x mean of y 
##  8.571429 11.888889
#Ketluan: Kết quả kiểm định t cho thấy không có sự khác biệt có ý nghĩa thống kê về tải trọng giữa nhóm A và nhóm B (p = 0.1345). Mặc dù nhóm B có tải trọng trung bình cao hơn (11.89 so với 8.57), nhưng sự khác biệt này không đủ mạnh để được xem là có ý nghĩa về mặt thống kê.
#2.5 Thực hiện bootstrap để đánh giá khác biệt trung bình tải trọng giữa 2 nhóm. Bạn nhận xét gì về kết quả này.
library(boot)
# Gộp dữ liệu và tạo nhãn nhóm
values <- c(group_A, group_B)
group <- c(rep("A", length(group_A)), rep("B", length(group_B)))
data <- data.frame(group, values)

# Hàm bootstrap: lấy hiệu số trung bình giữa nhóm B và A
boot_diff_mean <- function(data, indices) {
  d <- data[indices, ]
  mean(d$values[d$group == "B"]) - mean(d$values[d$group == "A"])
}

# Thực hiện bootstrap
set.seed(123)  # để kết quả tái lập được
boot_result <- boot(data = data, statistic = boot_diff_mean, R = 1000)

# Xem kết quả
boot_result
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = data, statistic = boot_diff_mean, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1*  3.31746 -0.02297878    2.025754
# Tính khoảng tin cậy 95%
boot.ci(boot_result, type = "perc")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_result, type = "perc")
## 
## Intervals : 
## Level     Percentile     
## 95%   (-0.498,  7.460 )  
## Calculations and Intervals on Original Scale
plot(boot_result)
abline(v = 0, col = "red", lty = 2)  # vẽ đường tại 0 để trực quan khoảng CI

#Ketluan:Phân tích bootstrap với 1000 mẫu lặp lại cho thấy trung bình tải trọng của nhóm B cao hơn nhóm A khoảng 3.32 đơn vị. Tuy nhiên, khoảng tin cậy 95% cho hiệu số trung bình là (-0.50 ; 7.46), cho thấy không có sự khác biệt có ý nghĩa thống kê giữa hai nhóm (vì khoảng tin cậy chứa 0). Điều này đồng nhất với kết quả kiểm định t-test truyền thống.

#2.6 Sử dụng ChatGPT viết code để phân tích sự khác biệt trung vị tải trọng giữa 2 nhóm dùng phương pháp bootstrap.