Bài tập thực hành ngày 2
NguyenThiVuThao
2025-10-22
#Việc 1. PHÂN TÍCH MÔ TẢ #1.1 Đọc dữ liệu “Obesity data.csv”
ob=read.csv("D:\\5_hoctap\\12_nckh\\Gs Nguyen Van Tuan Ts Tran Van Thach\\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) #nạp gói 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] |
#3.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: “tỉ lệ người bị cao huyết áp là hơn 50%, tỉ lệ người mắc bệnh tiểu đường là rất thấp 11.1%”. #Để trình bày kết quả tốt hơn ta sử dụng thêm 2 biến như sau:
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%) |
#4.4 Trình bày kết quả trung vị (Q1, Q3) thay vì 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á đặc điểm nào là khác biệt đáng kể (significant difference) giữa nam và nữ.
library(compareGroups) #nạp gói để có thể dùng các hàm so sánh
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 đọc dữ liệu 2 nhóm A,B vào R
A = c(14, 4, 10, 6, 3, 11, 12)
B = c(16, 17, 13, 12, 7, 16, 11, 8, 7)
#tạo biến wt = 14 4 10 6 3 11 12 16 17 13 12 7 16 11 8 7
wt = c(A, B)
# tạo biến group = "A" "A" "A" "A" "A" "A" "A" "B" "B" "B" "B" "B" "B" "B" "B" "B"
group = c(rep("A", 7), rep("B", 9))
#tạo bảng dữ liệu df có 2 biến là wt, group và 16 quan sát tương ứng
df = data.frame(wt, group)
dim(df) #trả về số dòng số cột của bảng df## [1] 16 2#2.2 Đánh giá phân bố: “Mặc dù sơ đồ không có hình chuông rõ ràng (có thể do số lượng mẫu), nhưng dùng shapiro test cho chỉ số p=0.7>0.05 nên có thể nhận xét wt có thể được coi là phân bố chuẩn”. Sơ đồ và kết quả như sau:
library(lessR) #nạp gói 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':
##
## label#dùng sơ đồ
Histogram(wt, data = df)## >>> Note: wt is not in a data frame (table)
## >>> Note: wt is not in a data frame (table)
## >>> Suggestions
## bin_width: set the width of each bin
## bin_start: set the start of the first bin
## bin_end: set the end of the last bin
## Histogram(wt, density=TRUE) # smoothed curve + histogram
## Plot(wt) # Violin/Box/Scatterplot (VBS) plot
##
## --- wt ---
##
## n miss mean sd min mdn max
## 16 0 10.44 4.29 3.00 11.00 17.00
##
##
## No (Box plot) outliers
##
##
## Bin Width: 2
## Number of Bins: 8
##
## Bin Midpnt Count Prop Cumul.c Cumul.p
## -------------------------------------------------
## 2 > 4 3 2 0.12 2 0.12
## 4 > 6 5 1 0.06 3 0.19
## 6 > 8 7 3 0.19 6 0.38
## 8 > 10 9 1 0.06 7 0.44
## 10 > 12 11 4 0.25 11 0.69
## 12 > 14 13 2 0.12 13 0.81
## 14 > 16 15 2 0.12 15 0.94
## 16 > 18 17 1 0.06 16 1.00
## #Dùng test
shapiro.test(df$wt)##
## Shapiro-Wilk normality test
##
## data: df$wt
## W = 0.96213, p-value = 0.7006#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
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ề kết quả:giá trị trung bình nhóm B cao hơn nhóm A nhưng giá trị chỉ số p = 0.1345 > 0.05 nên không có ý nghĩa thống kê hay không thể kết luận rằng tải trọng giữa hai nhóm là khác biệt rõ rệt.”
#2.5 dùng bootstra Đánh giá khác biệt trung bình tải trọng giữa 2 nhóm
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.187, 0.539 ) (-7.127, 0.556 )
##
## Level Percentile BCa
## 95% (-7.190, 0.492 ) (-7.253, 0.473 )
## 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. #PROMT = “Tôi có dữ liệu
cho 2 nhóm như sau: A = 14, 4, 10, 6, 3, 11, 12 B = 16, 17, 13, 12, 7,
16, 11, 8, 7 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 một cách đơn giản dễ hiểu nhất?
# Dữ liệu
A <- c(14, 4, 10, 6, 3, 11, 12)
B <- c(16, 17, 13, 12, 7, 16, 11, 8, 7)
# Số lần lặp bootstrap
n_boot <- 10000
# Hàm bootstrap cho hiệu trung bình
bootstrap_diff <- function(A, B, n_boot) {
diffs <- numeric(n_boot)
nA <- length(A)
nB <- length(B)
for (i in 1:n_boot) {
sample_A <- sample(A, nA, replace = TRUE)
sample_B <- sample(B, nB, replace = TRUE)
diffs[i] <- mean(sample_A) - mean(sample_B)
}
return(diffs)
}
# Chạy bootstrap
diffs <- bootstrap_diff(A, B, n_boot)
# Tính khoảng tin cậy 95% cho hiệu trung bình
ci <- quantile(diffs, c(0.025, 0.975))
# Tính p-value giả sử H0: hiệu trung bình = 0
# p-value = tỷ lệ mẫu bootstrap có hiệu trung bình <= 0 (hai phía)
p_value <- 2 * min(mean(diffs <= 0), mean(diffs >= 0))
# Kết quả
cat("Khoảng tin cậy 95% cho hiệu trung bình (A - B):", ci, "\n")## Khoảng tin cậy 95% cho hiệu trung bình (A - B): -7.063492 0.3809524cat("p-value kiểm định sự khác biệt trung bình giữa 2 nhóm:", p_value, "\n")## p-value kiểm định sự khác biệt trung bình giữa 2 nhóm: 0.0788