Phân tích dữ liệu bằng ngôn ngữ R - Viện Phương Nam (10-15/5/2025)
Thach Tran
2025-04-21
Ngày 3: Phân tích so sánh
Phân tích mô tả
Việc 1. Đọc dữ liệu vào R
bw = read.csv("C:\\Thach\\VN trips\\2025_2May\\Phuong Nam Institute\\Datasets\\birthwt.csv")Việc 2. Phân tích mô tả
2.1 Mô tả đặc điểm tuổi mẹ, cân nặng mẹ và cân nặng con
library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-table1(~ age + lwt + bwt, data = bw)| Overall (N=189) |
|
|---|---|
| age | |
| Mean (SD) | 23.2 (5.30) |
| Median [Min, Max] | 23.0 [14.0, 45.0] |
| lwt | |
| Mean (SD) | 130 (30.6) |
| Median [Min, Max] | 121 [80.0, 250] |
| bwt | |
| Mean (SD) | 2940 (729) |
| Median [Min, Max] | 2980 [709, 4990] |
[ChatGPT]
Prompt: “Sử dụng gói lệnh ‘table1’ để mô tả đặc điểm tuổi mẹ (age), cân nặng mẹ (lwt) và cân nặng con (bwt)”
table1(~ age + lwt + bwt, data = bw2)
2.2 Mô tả đặc điểm tuổi mẹ, cân nặng mẹ, tình trạng hút thuốc trong thai kỳ, chủng tộc và cân nặng con theo tình trạng trẻ sanh thiếu cân
bw$race.c = as.character(factor(bw$race, levels = c(1, 2, 3), labels = c("White", "Black", "Other")))
bw$smoke.c = as.character(factor(bw$smoke, levels = c(0, 1), labels = c("No", "Yes")))
bw$low.c = as.character(factor(bw$low, levels = c(0, 1), labels = c("No LBW", "Low birthweight")))
table1(~ age + lwt + smoke.c + race.c + bwt | low.c, data = bw)| Low birthweight (N=59) |
No LBW (N=130) |
Overall (N=189) |
|
|---|---|---|---|
| age | |||
| Mean (SD) | 22.3 (4.51) | 23.7 (5.58) | 23.2 (5.30) |
| Median [Min, Max] | 22.0 [14.0, 34.0] | 23.0 [14.0, 45.0] | 23.0 [14.0, 45.0] |
| lwt | |||
| Mean (SD) | 122 (26.6) | 133 (31.7) | 130 (30.6) |
| Median [Min, Max] | 120 [80.0, 200] | 124 [85.0, 250] | 121 [80.0, 250] |
| smoke.c | |||
| No | 29 (49.2%) | 86 (66.2%) | 115 (60.8%) |
| Yes | 30 (50.8%) | 44 (33.8%) | 74 (39.2%) |
| race.c | |||
| Black | 11 (18.6%) | 15 (11.5%) | 26 (13.8%) |
| Other | 25 (42.4%) | 42 (32.3%) | 67 (35.4%) |
| White | 23 (39.0%) | 73 (56.2%) | 96 (50.8%) |
| bwt | |||
| Mean (SD) | 2100 (391) | 3330 (478) | 2940 (729) |
| Median [Min, Max] | 2210 [709, 2500] | 3270 [2520, 4990] | 2980 [709, 4990] |
[ChatGPT]
PROMPT: “Sử dụng gói lệnh ‘table1’ để mô tả đặc điểm tuổi của mẹ (age), cân nặng của mẹ (lwt), tình trạng hút thuốc trong thai kỳ (smoke), chủng tộc (race), và cân nặng của con (bwt) theo tình trạng trẻ thiếu cân (low)”
ChatGPT:
Gắn nhãn mô tả cho các biến phân loại (giúp bảng dễ hiểu hơn):
bw2\(low <- factor(bw2\)low, levels = c(0, 1), labels = c(“Normal weight”, “Low birth weight”))
bw2\(smoke <- factor(bw2\)smoke, levels = c(0, 1), labels = c(“Non-smoker”, “Smoker”))
bw2\(race <- factor(bw2\)race, levels = c(1, 2, 3), labels = c(“White”, “Black”, “Other”))
Tạo bảng mô tả bằng table1():
table1(~ age + lwt + smoke + race + bwt | low, data = bw2)
So sánh 2 nhóm - biến liên tục
Việc 3. So sánh mật độ xương giữa nam và nữ
3.1 Đọc dữ liệu vào R
df = read.csv("C:\\Thach\\VN trips\\2025_2May\\Phuong Nam Institute\\Datasets\\Bone data.csv")
dim(df)## [1] 2162 9head(df)## id sex age weight height prior.fx fnbmd smoking fx
## 1 1 Male 73 98 175 0 1.08 1 0
## 2 2 Female 68 72 166 0 0.97 0 0
## 3 3 Male 68 87 184 0 1.01 0 0
## 4 4 Female 62 72 173 0 0.84 1 0
## 5 5 Male 61 72 173 0 0.81 1 0
## 6 6 Female 76 57 156 0 0.74 0 03.2 Vẽ histogram đánh giá phân bố mật độ xương
library(lessR)## Warning: package 'lessR' was built under R version 4.3.3##
## lessR 4.3.9 feedback: gerbing@pdx.edu
## --------------------------------------------------------------
## > d <- Read("") Read text, Excel, SPSS, SAS, or R data file
## 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, and descriptive statistics 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':
##
## labelHistogram(fnbmd, fill = "blue", xlab = "Mật độ xương (g/cm2)", ylab = "Frequency", data = df)
## >>> 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(fnbmd, density=TRUE) # smoothed curve + histogram
## Plot(fnbmd) # Violin/Box/Scatterplot (VBS) plot
##
## --- fnbmd ---
##
## n miss mean sd min mdn max
## 2122 40 0.829 0.155 0.280 0.820 1.510
##
##
## --- Outliers --- from the box plot: 33
##
## Small Large
## ----- -----
## 0.3 1.5
## 0.3 1.5
## 0.4 1.4
## 0.4 1.4
## 0.4 1.4
## 0.4 1.4
## 0.4 1.4
## 0.4 1.4
## 0.4 1.3
## 0.4 1.3
## 0.4 1.3
## 1.3
## 1.3
## 1.3
## 1.3
## 1.2
## 1.2
## 1.2
##
## + 15 more outliers
##
##
## Bin Width: 0.1
## Number of Bins: 14
##
## Bin Midpnt Count Prop Cumul.c Cumul.p
## ---------------------------------------------------
## 0.2 > 0.3 0.25 1 0.00 1 0.00
## 0.3 > 0.4 0.35 9 0.00 10 0.00
## 0.4 > 0.5 0.45 15 0.01 25 0.01
## 0.5 > 0.6 0.55 103 0.05 128 0.06
## 0.6 > 0.7 0.65 306 0.14 434 0.20
## 0.7 > 0.8 0.75 522 0.24 956 0.44
## 0.8 > 0.9 0.85 534 0.25 1490 0.69
## 0.9 > 1.0 0.95 371 0.17 1861 0.86
## 1.0 > 1.1 1.05 183 0.08 2044 0.95
## 1.1 > 1.2 1.15 48 0.02 2092 0.97
## 1.2 > 1.3 1.25 21 0.01 2113 0.98
## 1.3 > 1.4 1.35 6 0.00 2119 0.98
## 1.4 > 1.5 1.45 2 0.00 2121 0.98
## 1.5 > 1.6 1.55 1 0.00 2122 0.983.3 So sánh mật độ xương cổ xương đùi giữa nam và nữ
library(table1)
table1(~ fnbmd | sex, data = df)| Female (N=1317) |
Male (N=845) |
Overall (N=2162) |
|
|---|---|---|---|
| fnbmd | |||
| Mean (SD) | 0.778 (0.132) | 0.910 (0.153) | 0.829 (0.155) |
| Median [Min, Max] | 0.770 [0.280, 1.31] | 0.900 [0.340, 1.51] | 0.820 [0.280, 1.51] |
| Missing | 17 (1.3%) | 23 (2.7%) | 40 (1.9%) |
t.test(fnbmd ~ sex, data = df)##
## Welch Two Sample t-test
##
## data: fnbmd by sex
## t = -20.407, df = 1561, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
## 95 percent confidence interval:
## -0.1448770 -0.1194686
## sample estimates:
## mean in group Female mean in group Male
## 0.7775231 0.9096959Dùng gói lessR
library(lessR)
ttest(fnbmd ~ sex, data = df)##
## Compare fnbmd across sex with levels Male and Female
## Grouping Variable: sex
## Response Variable: fnbmd
##
##
## ------ Describe ------
##
## fnbmd for sex Male: n.miss = 23, n = 822, mean = 0.910, sd = 0.153
## fnbmd for sex Female: n.miss = 17, n = 1300, mean = 0.778, sd = 0.132
##
## Mean Difference of fnbmd: 0.132
##
## Weighted Average Standard Deviation: 0.141
##
##
## ------ Assumptions ------
##
## Note: These hypothesis tests can perform poorly, and the
## t-test is typically robust to violations of assumptions.
## Use as heuristic guides instead of interpreting literally.
##
## Null hypothesis, for each group, is a normal distribution of fnbmd.
## Group Male: Sample mean assumed normal because n > 30, so no test needed.
## Group Female: Sample mean assumed normal because n > 30, so no test needed.
##
## Null hypothesis is equal variances of fnbmd, homogeneous.
## Variance Ratio test: F = 0.023/0.018 = 1.336, df = 821;1299, p-value = 0.000
## Levene's test, Brown-Forsythe: t = 3.449, df = 2120, p-value = 0.001
##
##
## ------ Infer ------
##
## --- Assume equal population variances of fnbmd for each sex
##
## t-cutoff for 95% range of variation: tcut = 1.961
## Standard Error of Mean Difference: SE = 0.006
##
## Hypothesis Test of 0 Mean Diff: t-value = 21.080, df = 2120, p-value = 0.000
##
## Margin of Error for 95% Confidence Level: 0.012
## 95% Confidence Interval for Mean Difference: 0.120 to 0.144
##
##
## --- Do not assume equal population variances of fnbmd for each sex
##
## t-cutoff: tcut = 1.961
## Standard Error of Mean Difference: SE = 0.006
##
## Hypothesis Test of 0 Mean Diff: t = 20.407, df = 1560.981, p-value = 0.000
##
## Margin of Error for 95% Confidence Level: 0.013
## 95% Confidence Interval for Mean Difference: 0.119 to 0.145
##
##
## ------ Effect Size ------
##
## --- Assume equal population variances of fnbmd for each sex
##
## Standardized Mean Difference of fnbmd, Cohen's d: 0.939
##
##
## ------ Practical Importance ------
##
## Minimum Mean Difference of practical importance: mmd
## Minimum Standardized Mean Difference of practical importance: msmd
## Neither value specified, so no analysis
##
##
## ------ Graphics Smoothing Parameter ------
##
## Density bandwidth for sex Male: 0.044
## Density bandwidth for sex Female: 0.034
3.4 [ChatGPT]
PROMPT 1: “Dữ liệu ‘Bone data.csv’ đã được đọc vào R và đặt tên là ‘df’. Cung cấp mã R để so sánh mật độ xương (fnbmd) giữa nam và nữ (sex). Trình bày giá trị trung bình (và độ lệch chuẩn) của mật độ xương của 2 nhóm, trung bình khác biệt (và khoảng tin cậy 95%), và giá trị P so sánh khác biệt mật độ xương giữa nam và nữ”
PROMPT 2: “Lệnh trên quá phức tạp. Bạn có thể sử dụng gói lệnh ‘table1’ để mô tả mật độ xương giữa nam và nữ, và lệnh t.test để so sánh được không?”
ChatGPT:
Sử dụng gói table1 để mô tả mật độ xương theo giới tính
table1(~ fnbmd | sex, data = df)
Thực hiện kiểm định t-test để so sánh sự khác biệt mật độ xương giữa nam và nữ
t_test_result <- t.test(fnbmd ~ sex, data = df)
In kết quả kiểm định t-test
print(t_test_result)
PROMPT 3: “Bạn có thể sử dụng lệnh ‘ttest’ trong gói lệnh”lessR” để so sánh mật độ xương (fnbmd) giữa nam và nữ (sex) trong khi có đủ kết quả mô tả không?”
ChatGPT:
ttest(fnbmd ~ sex, data = df, descriptives = TRUE)
PROMPT 4: “Bạn hãy viết 1 câu để diễn giải kết quả trong file ‘fnbmd and sex.jpg’ vừa uploaded”
ChatGPT:
Dưới đây là câu diễn giải kết quả từ hình ảnh bạn đã tải lên:
“Kết quả kiểm định t-test cho thấy sự khác biệt có ý nghĩa thống kê giữa mật độ xương (fnbmd) ở nam và nữ, với giá trị p = 0.000. Trung bình mật độ xương của nam là 0.910 (SD = 0.153) và của nữ là 0.778 (SD = 0.132), với sự khác biệt trung bình giữa hai nhóm là 0.132 (Khoảng tin cậy 95%: 0.119 đến 0.145). Điều này chỉ ra rằng nam có mật độ xương cao hơn nữ một cách có ý nghĩa.”
Việc 4. Đánh giá ảnh hưởng của cafe lên RER
4.1 Nhập nhanh dữ liệu RER
placebo = c(105, 119, 100, 97, 96, 101, 94, 95, 98)
coffee = c(96, 99, 94, 89, 96, 93, 88, 105, 88)4.2 Đánh giá ảnh hưởng của cafe lên RER bằng kiểm định t
t.test(placebo, coffee)##
## Welch Two Sample t-test
##
## data: placebo and coffee
## t = 1.9948, df = 14.624, p-value = 0.06505
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.4490961 13.1157627
## sample estimates:
## mean of x mean of y
## 100.55556 94.222224.3 Đánh giá ảnh hưởng của cafe lên RER bằng bootstrap test
library(simpleboot)## Simple Bootstrap Routines (1.1-7)library(boot)
b = two.boot(placebo, coffee, mean, R = 500)
boot.ci(b)## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 500 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = b)
##
## Intervals :
## Level Normal Basic
## 95% ( 0.557, 12.318 ) ( 0.556, 12.000 )
##
## Level Percentile BCa
## 95% ( 0.667, 12.111 ) ( 0.861, 12.691 )
## Calculations and Intervals on Original Scale
## Some BCa intervals may be unstablehist(b, breaks = 50)
4.4 [ChatGPT]
PROMPT 1: “Tôi có dữ liệu của 2 nhóm như sau: placebo = c(105, 119, 100, 97, 96, 101, 94, 95, 98) và coffee = c(96, 99, 94, 89, 96, 93, 88, 105, 88). Bạn giúp viết lệnh R để phân tích sự khác biệt giữa 2 nhóm bằng phương pháp bootstrap không?”
ChatGPT:
Để phân tích sự khác biệt giữa hai nhóm bằng phương pháp bootstrap trong R, bạn có thể sử dụng một quy trình đơn giản để tính toán trung bình của sự khác biệt giữa các nhóm, và sau đó sử dụng phương pháp bootstrap để ước tính phân phối của sự khác biệt trung bình. Dưới đây là lệnh R để thực hiện phân tích này:
Dữ liệu của 2 nhóm
placebo <- c(105, 119, 100, 97, 96, 101, 94, 95, 98)
coffee <- c(96, 99, 94, 89, 96, 93, 88, 105, 88)
Số lần bootstrap
n_bootstrap <- 10000
Hàm bootstrap để tính sự khác biệt trung bình giữa hai nhóm
bootstrap_diff <- function(placebo, coffee, n) { diff_means <- numeric(n) for (i in 1:n) { placebo_sample <- sample(placebo, length(placebo), replace = TRUE) coffee_sample <- sample(coffee, length(coffee), replace = TRUE) diff_means[i] <- mean(placebo_sample) - mean(coffee_sample) } return(diff_means) }
Thực hiện bootstrap
bootstrap_results <- bootstrap_diff(placebo, coffee, n_bootstrap)
Tính toán khoảng tin cậy 95% từ kết quả bootstrap
ci_lower <- quantile(bootstrap_results, 0.025)
ci_upper <- quantile(bootstrap_results, 0.975)
In kết quả bootstrap
cat(“Khoảng tin cậy 95% của sự khác biệt trung bình là từ”, ci_lower, “đến”, ci_upper, “”)
Vẽ phân phối của sự khác biệt trung bình từ bootstrap
hist(bootstrap_results, main = “Phân phối Bootstrap của Sự Khác Biệt Trung Bình”, xlab = “Sự khác biệt trung bình”, col = “lightblue”, border = “black”)
PROMP 2: “Quá phức tạp. Bạn có thể sử dụng lệnh ‘two.boot’ trong gói lệnh ‘simpleboot’ hoặc ‘boot’ để thực hiện bootstrap không?”
ChatGPT:
Để thực hiện bootstrap đơn giản hơn giữa hai nhóm placebo và coffee bằng cách sử dụng lệnh two.boot từ gói lệnh simpleboot (hoặc boot), bạn có thể làm theo mã dưới đây:
library(simpleboot)
placebo <- c(105, 119, 100, 97, 96, 101, 94, 95, 98)
coffee <- c(96, 99, 94, 89, 96, 93, 88, 105, 88)
Thực hiện bootstrap với lệnh two.boot
bootstrap_result <- two.boot(placebo, coffee, nboot = 10000)
print(bootstrap_result)
Tính khoảng tin cậy 95% từ kết quả bootstrap
ci <- quantile(bootstrap_result$thetastar, c(0.025, 0.975))
So sánh 2 nhóm - biến phân loại
Việc 5. So sánh tỉ lệ gãy xương giữa nam và nữ
5.1 So sánh tỉ lệ gãy xương giữa nam và nữ
table1(~ as.factor(fx) | sex, data = df)| Female (N=1317) |
Male (N=845) |
Overall (N=2162) |
|
|---|---|---|---|
| as.factor(fx) | |||
| 0 | 916 (69.6%) | 701 (83.0%) | 1617 (74.8%) |
| 1 | 401 (30.4%) | 144 (17.0%) | 545 (25.2%) |
chisq.test(df$fx, df$sex)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: df$fx and df$sex
## X-squared = 48.363, df = 1, p-value = 0.0000000000035425.2 [ChatGPT]
PROMPT: “So sánh tỉ lệ gãy xương (fx) giữa nam và nữ (sex)”
ChatGPT:
Dưới đây là mã R để so sánh tỉ lệ gãy xương (fx) giữa nam và nữ (sex) sử dụng gói table1 để mô tả tỉ lệ gãy xương cho từng nhóm và lệnh chisq.test để so sánh tỉ lệ này giữa hai nhóm:
library(table1)
Sử dụng gói table1 để mô tả tỉ lệ gãy xương theo giới tính
table1(~ fx | sex, data = df)
Thực hiện kiểm định Chi-squared để so sánh tỉ lệ gãy xương giữa nam và nữ
chisq_result <- chisq.test(table(df\(fx, df\)sex))
In kết quả kiểm định chi-squared
print(chisq_result)