#viec1
#1.1 Đọc dữ liệu vào R
ob = read.csv ("D:/2025_TAP HUAN_THONG KE DL NCKH/Obesity data.csv")

#1.2 Mô tả dữ liệu
library(table1)
## 
## Attaching package: 'table1'
## The following objects are masked from 'package:base':
## 
##     units, units<-
table1(~ age + gender + weight + height + pcfat + hypertension + diabetes, data = ob)
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Min, Max] 48.0 [13.0, 88.0]
gender
F 862 (70.8%)
M 355 (29.2%)
weight
Mean (SD) 55.1 (9.40)
Median [Min, Max] 54.0 [34.0, 95.0]
height
Mean (SD) 157 (7.98)
Median [Min, Max] 155 [136, 185]
pcfat
Mean (SD) 31.6 (7.18)
Median [Min, Max] 32.4 [9.20, 48.4]
hypertension
Mean (SD) 0.507 (0.500)
Median [Min, Max] 1.00 [0, 1.00]
diabetes
Mean (SD) 0.111 (0.314)
Median [Min, Max] 0 [0, 1.00] 
#1.3. Nhận xét kq
ob$hyper = as.factor(ob$hypertension)
ob$diab = as.factor(ob$diabetes)

table1(~ age + gender + weight + height + pcfat + hypertension + hyper + diabetes + diab, data = ob)
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Min, Max] 48.0 [13.0, 88.0]
gender
F 862 (70.8%)
M 355 (29.2%)
weight
Mean (SD) 55.1 (9.40)
Median [Min, Max] 54.0 [34.0, 95.0]
height
Mean (SD) 157 (7.98)
Median [Min, Max] 155 [136, 185]
pcfat
Mean (SD) 31.6 (7.18)
Median [Min, Max] 32.4 [9.20, 48.4]
hypertension
Mean (SD) 0.507 (0.500)
Median [Min, Max] 1.00 [0, 1.00]
hyper
0 600 (49.3%)
1 617 (50.7%)
diabetes
Mean (SD) 0.111 (0.314)
Median [Min, Max] 0 [0, 1.00]
diab
0 1082 (88.9%)
1 135 (11.1%) 
#1.4. Trung vị (Q1, Q3)
table1(~ age + weight + height + pcfat, data = ob, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]"))
Overall
(N=1217)
age
Mean (SD) 47.2 (17.3)
Median [Q1, Q3] 48.0 [35.0, 58.0]
weight
Mean (SD) 55.1 (9.40)
Median [Q1, Q3] 54.0 [49.0, 61.0]
height
Mean (SD) 157 (7.98)
Median [Q1, Q3] 155 [151, 162]
pcfat
Mean (SD) 31.6 (7.18)
Median [Q1, Q3] 32.4 [27.0, 36.8] 
#1.5. Mô tả tuổi, cân nặng, chiều cao, tỉ trọng mỡ, tiền bệnh cao huyết áp theo giới tính
table1(~ age + weight + height + pcfat + hyper + diab | gender, data = ob)
F
(N=862)		 M
(N=355)		 Overall
(N=1217)
age
Mean (SD) 48.6 (16.4) 43.7 (18.8) 47.2 (17.3)
Median [Min, Max] 49.0 [14.0, 85.0] 44.0 [13.0, 88.0] 48.0 [13.0, 88.0]
weight
Mean (SD) 52.3 (7.72) 62.0 (9.59) 55.1 (9.40)
Median [Min, Max] 51.0 [34.0, 95.0] 62.0 [38.0, 95.0] 54.0 [34.0, 95.0]
height
Mean (SD) 153 (5.55) 165 (6.73) 157 (7.98)
Median [Min, Max] 153 [136, 170] 165 [146, 185] 155 [136, 185]
pcfat
Mean (SD) 34.7 (5.19) 24.2 (5.76) 31.6 (7.18)
Median [Min, Max] 34.7 [14.6, 48.4] 24.6 [9.20, 39.0] 32.4 [9.20, 48.4]
hyper
0 430 (49.9%) 170 (47.9%) 600 (49.3%)
1 432 (50.1%) 185 (52.1%) 617 (50.7%)
diab
0 760 (88.2%) 322 (90.7%) 1082 (88.9%)
1 102 (11.8%) 33 (9.3%) 135 (11.1%) 
#1.6. Đánh giá khác biệt giữa nam và nữ
library(compareGroups)
createTable(compareGroups(gender ~ age + weight + height + pcfat + hyper + diab, data = ob))
## 
## --------Summary descriptives table by 'gender'---------
## 
## ________________________________________ 
##             F           M      p.overall 
##           N=862       N=355              
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## age    48.6 (16.4) 43.7 (18.8)  <0.001   
## weight 52.3 (7.72) 62.0 (9.59)  <0.001   
## height 153 (5.55)  165 (6.73)   <0.001   
## pcfat  34.7 (5.19) 24.2 (5.76)  <0.001   
## hyper:                           0.569   
##     0  430 (49.9%) 170 (47.9%)           
##     1  432 (50.1%) 185 (52.1%)           
## diab:                            0.238   
##     0  760 (88.2%) 322 (90.7%)           
##     1  102 (11.8%) 33 (9.30%)            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯
#Việc 2. 
#2.1 Có dữ liệu về tải trọng của 2 nhóm A và B:
#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 Xem xét tính phân bổ chuẩn (normal distribution)?
#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)

#2.3. Mô tả về tải trọng 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
#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