Bài thực hành ngày 2
Phan Thoai Chieu
2025-10-22
#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"
ob=read.csv("/cloud/project/Obesitydata.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)
ob$hyper=as.factor(ob$hypertension)
ob$diab=as.factor(ob$diabetes)
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~age+gender+weight+height+pcfat+hyper+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] |
| hyper | |
| 0 | 600 (49.3%) |
| 1 | 617 (50.7%) |
| diab | |
| 0 | 1082 (88.9%) |
| 1 | 135 (11.1%) |
#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.
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)
library(table1)
table1(~age+gender+weight+height+pcfat+hypertension+hyper+diabetes+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] |
| gender | |||
| F | 862 (100%) | 0 (0%) | 862 (70.8%) |
| M | 0 (0%) | 355 (100%) | 355 (29.2%) |
| 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] |
| hypertension | |||
| Mean (SD) | 0.501 (0.500) | 0.521 (0.500) | 0.507 (0.500) |
| Median [Min, Max] | 1.00 [0, 1.00] | 1.00 [0, 1.00] | 1.00 [0, 1.00] |
| hyper | |||
| 0 | 430 (49.9%) | 170 (47.9%) | 600 (49.3%) |
| 1 | 432 (50.1%) | 185 (52.1%) | 617 (50.7%) |
| diabetes | |||
| Mean (SD) | 0.118 (0.323) | 0.0930 (0.291) | 0.111 (0.314) |
| Median [Min, Max] | 0 [0, 1.00] | 0 [0, 1.00] | 0 [0, 1.00] |
| 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%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯#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 2#2.2 Kiểm tra phân bổ chuẩn (normal distribution)
shapiro.test(A)##
## Shapiro-Wilk normality test
##
## data: A
## W = 0.92541, p-value = 0.5126shapiro.test(B)##
## Shapiro-Wilk normality test
##
## data: B
## W = 0.89641, p-value = 0.2319#2.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.888889#Nhận xét: Vì P = 0.1345 > 0.05 nên không có sự khác biệt về tài trọng giữa A và B
#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(simpleboot)## Simple Bootstrap Routines (1.1-8)library(boot)
b=two.boot(A,B,mean,R=1000)
boot.ci(b)## Warning in boot.ci(b): bootstrap variances needed for studentized intervals## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b)
##
## Intervals :
## Level Normal Basic
## 95% (-7.153, 0.492 ) (-7.095, 0.476 )
##
## Level Percentile BCa
## 95% (-7.111, 0.460 ) (-7.248, 0.306 )
## Calculations and Intervals on Original Scalehist(b,breaks=50)
#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.
# Cài và nạp gói boot nếu chưa có
install.packages("boot") # chỉ cần cài 1 lần## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.5'
## (as 'lib' is unspecified)library(boot)
# Dữ liệu
A <- c(14, 4, 10, 6, 3, 11, 12)
B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)
# Gộp dữ liệu vào 1 vector và tạo chỉ số nhóm
data <- c(A, B)
group <- c(rep("A", length(A)), rep("B", length(B)))
# Hàm thống kê: chênh lệch trung bình giữa hai nhóm
mean_diff <- function(data, index) {
d <- data[index]
g <- group[index]
mean(d[g == "A"]) - mean(d[g == "B"])
}
# Thực hiện bootstrap (1000 mẫu ngẫu nhiên)
set.seed(123) # để kết quả có thể tái lập
boot_res <- boot(data = data, statistic = mean_diff, R = 1000)
# Kết quả bootstrap
boot_res##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = data, statistic = mean_diff, 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_res, conf = 0.95, type = c("norm", "basic", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_res, conf = 0.95, type = c("norm", "basic",
## "perc", "bca"))
##
## Intervals :
## Level Normal Basic
## 95% (-7.311, 0.630 ) (-7.133, 0.825 )
##
## Level Percentile BCa
## 95% (-7.460, 0.498 ) (-7.590, 0.374 )
## Calculations and Intervals on Original Scale