Bài tập ngày 2
Phượng
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”
file.choose()## [1] "F:\\VIEC_KHOA\\SU_KIEN\\2025_2026\\Taphuanphantichdulieu\\Datamophong\\Baitapngay2.html"ob=read.csv("F:\\VIEC_KHOA\\SU_KIEN\\2025_2026\\Taphuanphantichdulieu\\Datamophong\\Obesity data.csv")head(ob)## id gender height weight bmi age WBBMC wbbmd fat lean pcfat hypertension
## 1 1 F 150 49 21.8 53 1312 0.88 17802 28600 37.3 0
## 2 2 M 165 52 19.1 65 1309 0.84 8381 40229 16.8 1
## 3 3 F 157 57 23.1 64 1230 0.84 19221 36057 34.0 1
## 4 4 F 156 53 21.8 56 1171 0.80 17472 33094 33.8 1
## 5 5 M 160 51 19.9 54 1681 0.98 7336 40621 14.8 0
## 6 6 F 153 47 20.1 52 1358 0.91 14904 30068 32.2 1
## diabetes
## 1 1
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0##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] |
##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? Do giá trị mã hóa của tiền căn bệnh cao huyết áp và tiểu đường là 0 và 1 nên hệ thống hiểu nó là biến liên tục. Kết quả thống kê hiển thị Mean. Nhưng điều cần biết là số lượng và tỷ lệ mắc bệnh. => cần dùng hàm as.factor để chuuyển thành biến phân loại
table1(~age + gender + weight + height + pcfat + as.factor(hypertension) + as.factor(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%) |
| as.factor(diabetes) | |
| 0 | 1082 (88.9%) |
| 1 | 135 (11.1%) |
##1.4.Trình bày kết quả trung vị (Q1, Q3) và (min, max) cho biến liên tục.
my.render.cont <- function(x) {
q <- stats::quantile(x, probs = c(0, 0.25, 0.5, 0.75, 1), na.rm = TRUE)
m <- mean(x, na.rm = TRUE)
sdv <- sd(x, na.rm = TRUE)
out <- c(" " = "", # Dòng trống để giữ format
"Mean (SD)" = sprintf("%.1f (%.1f)", m, sdv),
"Median (Q1, Q3)" = sprintf("%.1f (%.1f, %.1f)", q[3], q[2], q[4]),
"Median (Min, Max)" = sprintf("%.1f [%.1f, %.1f]", q[3], q[1], q[5])
)
return(out)
}table1(~ age + gender + height + weight + pcfat + as.factor(hypertension) + as.factor(diabetes),data = ob,render.continuous = my.render.cont)| Overall (N=1217) |
|
|---|---|
| age | |
| Mean (SD) | 47.2 (17.3) |
| Median (Q1, Q3) | 48.0 (35.0, 58.0) |
| Median (Min, Max) | 48.0 [13.0, 88.0] |
| gender | |
| F | 862 (70.8%) |
| M | 355 (29.2%) |
| height | |
| Mean (SD) | 156.7 (8.0) |
| Median (Q1, Q3) | 155.0 (151.0, 162.0) |
| Median (Min, Max) | 155.0 [136.0, 185.0] |
| weight | |
| Mean (SD) | 55.1 (9.4) |
| Median (Q1, Q3) | 54.0 (49.0, 61.0) |
| Median (Min, Max) | 54.0 [34.0, 95.0] |
| pcfat | |
| Mean (SD) | 31.6 (7.2) |
| Median (Q1, Q3) | 32.4 (27.0, 36.8) |
| Median (Min, Max) | 32.4 [9.2, 48.4] |
| as.factor(hypertension) | |
| 0 | 600 (49.3%) |
| 1 | 617 (50.7%) |
| as.factor(diabetes) | |
| 0 | 1082 (88.9%) |
| 1 | 135 (11.1%) |
##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 + height + weight + pcfat + as.factor(hypertension) + as.factor(diabetes)|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] |
| height | |||
| Mean (SD) | 153 (5.55) | 165 (6.73) | 157 (7.98) |
| Median [Min, Max] | 153 [136, 170] | 165 [146, 185] | 155 [136, 185] |
| 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] |
| 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] |
| as.factor(hypertension) | |||
| 0 | 430 (49.9%) | 170 (47.9%) | 600 (49.3%) |
| 1 | 432 (50.1%) | 185 (52.1%) | 617 (50.7%) |
| as.factor(diabetes) | |||
| 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ữ Vì trong compareGroups không được dùng hàm như as.factor() trực tiếp.compareGroups() yêu cầu các biến phải có kiểu dữ liệu phù hợp (numeric hoặc factor) trước khi đưa vào công thức. Do đó, cần tạo biến định danh cho hypertension và diabetes
ob$hyper <- factor(ob$hypertension, labels = c("Không", "Có"))
ob$diab <- factor(ob$diabetes, labels = c("Không", "Có"))library(compareGroups)
createTable(compareGroups(gender~age + height + weight + 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
## height 153 (5.55) 165 (6.73) <0.001
## weight 52.3 (7.72) 62.0 (9.59) <0.001
## pcfat 34.7 (5.19) 24.2 (5.76) <0.001
## hyper: 0.569
## Không 430 (49.9%) 170 (47.9%)
## Có 432 (50.1%) 185 (52.1%)
## diab: 0.238
## Không 760 (88.2%) 322 (90.7%)
## Có 102 (11.8%) 33 (9.30%)
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯Trong bảng này, loại kiểm định tùy theo đặc điểm biến,cụ thể gender là biến định tính, age là biến định lượng => giá trị p là của kiểm định t, hyper là biến định tính => kiểm định chi-square
#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)- c là hàm combine, dùng để ghép các phần tử thành một vector => wt là dữ liệu tải trọng của cả 2 nhóm A và B => tạo wt sẽ ra 1 bảng có 16 quan sát
- group sẽ tạo nhãn cho dữ liệu. Do rep(x, times) lặp lại giá trị x một số lần=> rep(“A”, 7) tạo một vector gồm 7 phần tử “A”. rep(“B”, 9) tạo một vector gồm 9 phần tử “B”.
- data.frame() tạo một bảng dữ liệu (data frame) từ các vector.Ở đây, df sẽ có 2 cột:wt: chứa dữ liệu từ A và B. group: chứa nhãn “A” hoặc “B” tương ứng. Tạo bảng df để gộp dữ liệu của A và B, nhằm có dữ liệu để kiểm định phân phối chuẩn trước khi kiểm định t-test.
Kiểm tra bảng tạo ra:
df## wt group
## 1 14 A
## 2 4 A
## 3 10 A
## 4 6 A
## 5 3 A
## 6 11 A
## 7 12 A
## 8 16 B
## 9 17 B
## 10 13 B
## 11 12 B
## 12 7 B
## 13 16 B
## 14 11 B
## 15 8 B
## 16 7 Bdim(df)## [1] 16 2##2.2 Tải trọng có tuân theo phân bổ chuẩn (normal distribution) không? Tải gói LessR: install.packages(“lessR”) Trong R, khi gọi install.packages(), tên gói phải được đặt trong dấu ngoặc kép. Ví dụ: install.packages(“lessR”) gọi gói về:
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':
##
## labelVẽ biểu đồ: wt là tải trọng cần vẽ, lấy từ data là bảng df đã tạo
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 Shapiro.test để kiểm định phân phối chuẩn của wt.
shapiro.test(df$wt)##
## Shapiro-Wilk normality test
##
## data: df$wt
## W = 0.96213, p-value = 0.7006Kết luận: p=0.7006 > 0.5 => phân phối chuẩn ##2.3 Mô tả đặc điểm về tải trọng giữa 2 nhóm
table1(~wt,data=df)| Overall (N=16) |
|
|---|---|
| wt | |
| Mean (SD) | 10.4 (4.29) |
| Median [Min, Max] | 11.0 [3.00, 17.0] |
##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 có thể dùng 2 lệnh khác nhau, nhưng kết quả như nhau: thứ nhất: t.test(A,B). Lệnh này chỉ cần có dữ liệu của A và B là làm được. thứ hai: t.test(wt~group,data=df). Lệnh này dùng dữ liệu đã gộp thành bảng.
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.888889kết luận: p>0.05 => không có sự khác biệt. Hiểu 1 cách khác, khoảng ước lượng với độ tin cậy 95% cóchứa giá trị0 (từ -7,8 đến 1,18 sẽ đi qua giá trị 0)=> chênh lệch (difference=A-B)=0=> có khả năng A=B=> không khác biệt
##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. Nếu dữ liệu không phân phối chuẩn, sẽ thực hiện bootstrap thay cho t-test
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.042, 0.458 ) (-6.984, 0.634 )
##
## Level Percentile BCa
## 95% (-7.269, 0.349 ) (-7.536, 0.045 )
## Calculations and Intervals on Original ScaleCần nghe lại Thầy giảng để hiểu khi nào dùng Normal, basic, percentile, BCa
vẽ biểu đồ phân bố
hist(b,breaks=50)
b: Đối tượng chứa dữ liệu (ở đây là kết quả bootstrap từ two.boot()).
breaks=50: Số lượng khoảng (bins) để chia dữ liệu → histogram sẽ có 50
cột