Thuc hanh R - Ngay 2 - Thanh Phuong
Phuong Nguyen
2025-10-22
Ngày 2: So sánh 2 nhóm - Biến liên tục
Việc 1. Phân tích mô tả
1.1 Đọc dữ liệu vào R
ob = read.csv("E:\\TAP HUAN 21.10.2025\\Obesity data.csv")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)
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] |
table1(~ age + gender + weight + height + pcfat + as.factor(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] |
| as.factor(hypertension) | |
| 0 | 600 (49.3%) |
| 1 | 617 (50.7%) |
| diabetes | |
| Mean (SD) | 0.111 (0.314) |
| Median [Min, Max] | 0 [0, 1.00] |
1.3 Nhận xét 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?
Bệnh huyết áp và tiểu đường đang trình bày kết quả ra giá trị trung bình và độ lệch chuẩn. Để kết quả tốt hơn nên định nghĩa lại bằng as.factor để R hiểu đây là một biến phân loại và kết quả là đếm số lượng và tỷ lệ(%)
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 Trình bày kết quả trung vị (Q1, Q3) thay vi trung vị (min, max) cho biến liên tục.
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ả đặ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)
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á xem đặc điểm nào là khác biệt đáng kể (significant difference) 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%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯=>nam giới có tuổi, cân nặng, chiều cao và tỷ lệ mỡ cơ thể khác biệt đáng kể so với nữ giới (p < 0.001). Trong khi đó, tỷ lệ tăng huyết áp và đái tháo đường không có sự khác biệt có ý nghĩa thống kê giữa hai nhóm (p > 0.05).
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.
A = c(14, 4, 10, 6, 3, 11, 12)
B = c(16, 17, 13, 12, 7, 16, 11, 8, 7)
wt = c(A, B)
group = c(rep("A", 7), rep("B", 9))
df = data.frame(wt, group)
dim(df)## [1] 16 22.2 Tải trọng có tuân theo phân bổ chuẩn (normal distribution) không?
Đánh giá phân bổ
Dùng biểu đồ.
library(lessR)##
## lessR 4.4.5 feedback: gerbing@pdx.edu
## --------------------------------------------------------------
## > d <- Read("") Read data file, many formats available, e.g., Excel
## d is default data frame, data= in analysis routines optional
##
## Many examples of reading, writing, and manipulating data,
## graphics, testing means and proportions, regression, factor analysis,
## customization, forecasting, and aggregation from pivot tables
## Enter: browseVignettes("lessR")
##
## View lessR updates, now including time series forecasting
## Enter: news(package="lessR")
##
## Interactive data analysis
## Enter: interact()##
## Attaching package: 'lessR'## The following object is masked from 'package:table1':
##
## labelDùng test
shapiro.test(df$wt)##
## Shapiro-Wilk normality test
##
## data: df$wt
## W = 0.96213, p-value = 0.70062.3 Mô tả đặc điểm về tải trọng giữa 2 nhóm
library(table1)
table1(~ wt | group, data = df, render.continuous = c(. = "Mean (SD)", . = "Median [Q1, Q3]"))| A (N=7) |
B (N=9) |
Overall (N=16) |
|
|---|---|---|---|
| wt | |||
| Mean (SD) | 8.57 (4.24) | 11.9 (3.95) | 10.4 (4.29) |
| Median [Q1, Q3] | 10.0 [5.00, 11.5] | 12.0 [8.00, 16.0] | 11.0 [7.00, 13.3] |
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(A, B)##
## Welch Two Sample t-test
##
## data: A and 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.888889t.test(wt ~ group, data = df)##
## Welch Two Sample t-test
##
## data: wt by group
## t = -1.6, df = 12.554, p-value = 0.1345
## alternative hypothesis: true difference in means between group A and group B is not equal to 0
## 95 percent confidence interval:
## -7.813114 1.178194
## sample estimates:
## mean in group A mean in group B
## 8.571429 11.888889=> Hai kiểm định cho kết quả giống nhau.
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.
Dùng lệnh:
library(simpleboot)## Simple Bootstrap Routines (1.1-8)library(boot)
b = two.boot(A, B, mean, R = 1000)
boot.ci(b)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b)
##
## Intervals :
## Level Normal Basic
## 95% (-7.092, 0.352 ) (-7.174, 0.000 )
##
## Level Percentile BCa
## 95% (-6.635, 0.539 ) (-6.726, 0.380 )
## Calculations and Intervals on Original Scale=> Tất cả các khoảng tin cậy đều bao trùm giá trị 0. Do đó không thể kết luận trung bình của nhóm A và nhóm B khác nhau một cách có ý nghĩa thống kê (ở mức 95%).
Dùng biểu đồ:
hist(b, breaks = 50)
=> Trong biểu đồ, bootstrap bao trùm cả vùng trước và sau giá trị 0 → phù hợp với kết luận: không có khác biệt có ý nghĩa thống kê giữa hai nhóm
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.
PROMPT: Tôi có dữ liệu cho 2 nhóm như sau. Bạn có thể viết mã R để phân tích sự khác biệt giữa hai nhóm dùng bootstrap? A = c(14, 4, 10, 6, 3, 11, 12) B = c(16, 17, 13, 12, 7, 16, 11, 8, 7
set.seed(123) # để có thể tái lặp kết quả
A <- c(14, 4, 10, 6, 3, 11, 12)
B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)
# 1. Hiệu trung bình thực tế
obs_diff <- mean(A) - mean(B)
# 2. Bootstrap
B_iter <- 10000
boot_diff <- numeric(B_iter)
for (i in 1:B_iter) {
A_boot <- sample(A, replace = TRUE)
B_boot <- sample(B, replace = TRUE)
boot_diff[i] <- mean(A_boot) - mean(B_boot)
}
# 3. Khoảng tin cậy 95%
ci <- quantile(boot_diff, c(0.025, 0.975))
# 4. Kết quả
obs_diff## [1] -3.31746ci## 2.5% 97.5%
## -7.0476190 0.4285714# 5. Vẽ biểu đồ bootstrap
hist(boot_diff, main = "Bootstrap Distribution of Mean(A) - Mean(B)",
xlab = "Bootstrap Mean Differences", col = "lightblue")
abline(v = obs_diff, col = "red", lwd = 2)
PROMPT 2: Bạn có thể sử dụng gói lệnh simpleboot với lệnh two.boot và boot.ci thể thực hiện câu hỏi trên không?
install.packages("simpleboot") # nếu bạn chưa có
install.packages("boot") # boot.ci nằm trong gói boot
library(simpleboot)
library(boot)A <- c(14, 4, 10, 6, 3, 11, 12)
B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)set.seed(123)
boot_res <- two.boot(A, B,
FUN = mean,
R = 10000) # số lần bootstrapboot_res$t0 # Hiệu trung bình quan sát (mean(A) - mean(B))## [1] -3.31746head(boot_res$t) # Một vài giá trị bootstrap## [,1]
## [1,] -1.3650794
## [2,] -3.6507937
## [3,] -2.3174603
## [4,] -0.6666667
## [5,] -5.6666667
## [6,] -5.9523810ci_res <- boot.ci(boot_res, type = c("perc", "bca"))
ci_res## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_res, type = c("perc", "bca"))
##
## Intervals :
## Level Percentile BCa
## 95% (-7.063, 0.413 ) (-7.190, 0.317 )
## Calculations and Intervals on Original Scalehist(boot_res$t,
main = "Bootstrap Distribution of Mean(A) - Mean(B)",
xlab = "Bootstrap Mean Differences",
col = "lightblue")
abline(v = boot_res$t0, col = "red", lwd = 2)