#Việc 1. Đọc dữ liệu "birthwt.csv" vào R và gọi dữ liệu là "bw"
# Đọc dữ liệu từ file CSV
data <- read.csv("E:\\HOPT\\PTDLR\\DATA\\birthwt.csv")
# Xem 5 dòng đầu
head(data)
## id low age lwt race smoke ptl ht ui ftv bwt
## 1 85 0 19 182 2 0 0 0 1 0 2523
## 2 86 0 33 155 3 0 0 0 0 3 2551
## 3 87 0 20 105 1 1 0 0 0 1 2557
## 4 88 0 21 108 1 1 0 0 1 2 2594
## 5 89 0 18 107 1 1 0 0 1 0 2600
## 6 91 0 21 124 3 0 0 0 0 0 2622
#Việc 2. Phân tích mô tả:
# 2.1 Mô tả đặc điểm tuổi của mẹ (age), cân nặng của mẹ (lwt) và cân nặng của con (bwt)
# Thống kê mô tả cho age, lwt và bwt
summary(data[, c("age", "lwt", "bwt")])
## age lwt bwt
## Min. :14.00 Min. : 80.0 Min. : 709
## 1st Qu.:19.00 1st Qu.:110.0 1st Qu.:2414
## Median :23.00 Median :121.0 Median :2977
## Mean :23.24 Mean :129.8 Mean :2945
## 3rd Qu.:26.00 3rd Qu.:140.0 3rd Qu.:3487
## Max. :45.00 Max. :250.0 Max. :4990
#Kiểm tra phân phối bằng đồ thị (tuỳ chọn)
# Nếu chưa có gói ggplot2 thì cài đặt
# install.packages("ggplot2")
library(ggplot2)
# Histogram cho tuổi mẹ
ggplot(data, aes(x = age)) +
geom_histogram(binwidth = 1, fill = "skyblue", color = "black") +
labs(title = "Phân phối tuổi của mẹ", x = "Tuổi", y = "Tần suất")
# Histogram cho cân nặng mẹ
ggplot(data, aes(x = lwt)) +
geom_histogram(binwidth = 5, fill = "lightgreen", color = "black") +
labs(title = "Phân phối cân nặng của mẹ", x = "Cân nặng (pounds)", y = "Tần suất")
# Histogram cho cân nặng con
ggplot(data, aes(x = bwt)) +
geom_histogram(binwidth = 100, fill = "pink", color = "black") +
labs(title = "Phân phối cân nặng của con", x = "Cân nặng (grams)", y = "Tần suất")
#2.2 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)
# Đọc dữ liệu từ file CSV
data <- read.csv("E:\\HOPT\\PTDLR\\DATA\\birthwt.csv")
# Chuyển 'low', 'smoke', 'race' thành biến phân loại
data$low <- factor(data$low, labels = c("Bình thường", "Thiếu cân"))
data$smoke <- factor(data$smoke, labels = c("Không hút thuốc", "Hút thuốc"))
data$race <- factor(data$race, labels = c("White", "Black", "Other"))
#Thống kê mô tả theo tình trạng thiếu cân
# Thống kê trung bình và độ lệch chuẩn theo nhóm
# Thống kê mô tả cho biến số theo nhóm 'low'
aggregate(cbind(age, lwt, bwt) ~ low, data = data,
FUN = function(x) c(Mean = mean(x), SD = sd(x)))
## low age.Mean age.SD lwt.Mean lwt.SD bwt.Mean bwt.SD
## 1 Bình thường 23.661538 5.584522 133.30000 31.72402 3329.1077 478.4160
## 2 Thiếu cân 22.305085 4.511496 122.13559 26.55928 2097.3390 391.0788
#Tần suất và tỷ lệ phần trăm cho biến phân loại
# Tần suất tình trạng hút thuốc theo thiếu cân
table_smoke <- table(data$low, data$smoke)
prop_smoke <- prop.table(table_smoke, margin = 1) * 100
# Tần suất chủng tộc theo thiếu cân
table_race <- table(data$low, data$race)
prop_race <- prop.table(table_race, margin = 1) * 100
# In ra kết quả
table_smoke
##
## Không hút thuốc Hút thuốc
## Bình thường 86 44
## Thiếu cân 29 30
round(prop_smoke, 1)
##
## Không hút thuốc Hút thuốc
## Bình thường 66.2 33.8
## Thiếu cân 49.2 50.8
table_race
##
## White Black Other
## Bình thường 73 15 42
## Thiếu cân 23 11 25
round(prop_race, 1)
##
## White Black Other
## Bình thường 56.2 11.5 32.3
## Thiếu cân 39.0 18.6 42.4
#vễ mô tả info trên dạng Biểu đồ Boxplot
#So sánh các biến liên tục: tuổi mẹ (age), cân nặng mẹ (lwt), cân nặng con (bwt)
library(ggplot2)
# Boxplot - Tuổi mẹ theo nhóm thiếu cân
ggplot(data, aes(x = low, y = age, fill = low)) +
geom_boxplot() +
labs(title = "Tuổi của mẹ theo tình trạng thiếu cân của trẻ",
x = "Tình trạng trẻ", y = "Tuổi mẹ") +
theme_minimal()
# Boxplot - Cân nặng mẹ theo nhóm thiếu cân
ggplot(data, aes(x = low, y = lwt, fill = low)) +
geom_boxplot() +
labs(title = "Cân nặng của mẹ theo tình trạng thiếu cân của trẻ",
x = "Tình trạng trẻ", y = "Cân nặng mẹ (lwt)") +
theme_minimal()
# Boxplot - Cân nặng trẻ
ggplot(data, aes(x = low, y = bwt, fill = low)) +
geom_boxplot() +
labs(title = "Cân nặng trẻ sơ sinh theo nhóm thiếu cân",
x = "Tình trạng trẻ", y = "Cân nặng con (bwt)") +
theme_minimal()
#mô tả bằng Biểu đồ Cột
#So sánh tỷ lệ phần trăm hút thuốc và chủng tộc giữa hai nhóm
#iểu đồ tỷ lệ hút thuốc theo tình trạng thiếu cân:
ggplot(data, aes(x = smoke, fill = low)) +
geom_bar(position = "fill") +
labs(title = "Tỷ lệ hút thuốc trong thai kỳ theo tình trạng thiếu cân",
x = "Tình trạng hút thuốc", y = "Tỷ lệ (%)") +
scale_y_continuous(labels = scales::percent) +
theme_minimal()
# Biểu đồ tỷ lệ chủng tộc theo tình trạng thiếu cân:
ggplot(data, aes(x = race, fill = low)) +
geom_bar(position = "fill") +
labs(title = "Tỷ lệ chủng tộc theo tình trạng thiếu cân của trẻ",
x = "Chủng tộc", y = "Tỷ lệ (%)") +
scale_y_continuous(labels = scales::percent) +
theme_minimal()
#So sánh 2 nhóm – biến liên tục
#Việc 3. So sánh mật độ xương tại cổ xương đùi giữa nam và nữ
#3.1 Đọc dữ liệu "Bone data.csv" vào R và gọi dữ liệu là “df”
dbone <- read.csv("E:\\HOPT\\PTDLR\\DATA\\Bone data.csv")
head(dbone)
## 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 0
str(dbone)
## 'data.frame': 2162 obs. of 9 variables:
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## $ sex : chr "Male" "Female" "Male" "Female" ...
## $ age : int 73 68 68 62 61 76 63 64 76 64 ...
## $ weight : int 98 72 87 72 72 57 97 85 48 89 ...
## $ height : int 175 166 184 173 173 156 173 167 153 166 ...
## $ prior.fx: int 0 0 0 0 0 0 0 0 0 0 ...
## $ fnbmd : num 1.08 0.97 1.01 0.84 0.81 0.74 1.01 0.86 0.65 0.98 ...
## $ smoking : int 1 0 0 1 1 0 1 0 0 0 ...
## $ fx : int 0 0 0 0 0 0 0 0 0 0 ...
#3.2 Vẽ biểu đồ histogram đánh giá phân bố mật độ xương tại cổ xương đùi. Bạn đánh giá phân bố mật độ xương như thế nào?
# Nếu chưa cài ggplot2
# install.packages("ggplot2")
library(ggplot2)
# Histogram mật độ xương cổ xương đùi
ggplot(dbone, aes(x = fnbmd)) +
geom_histogram(binwidth = 0.05, fill = "steelblue", color = "black") +
labs(title = "Phân bố mật độ xương tại cổ xương đùi",
x = "Mật độ xương (fnbmd)", y = "Tần suất") +
theme_minimal()
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_bin()`).
#Đánh giá phân bố
# Histogram với đường mật độ
ggplot(dbone, aes(x = fnbmd)) +
geom_histogram(aes(y = ..density..), binwidth = 0.05,
fill = "lightblue", color = "black") +
geom_density(color = "red", size = 1) +
labs(title = "Phân bố mật độ xương tại cổ xương đùi (kèm mật độ)",
x = "Mật độ xương (fnbmd)", y = "Mật độ") +
theme_minimal()
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_density()`).
#Gợi ý đánh giá phân bố
#Nếu biểu đồ có hình chuông đối xứng → Phân bố chuẩn (Normal distribution).
#Nếu lệch trái hoặc phải → Phân bố lệch (skewed).
#Có thể sử dụng kiểm định Shapiro-Wilk để kiểm tra phân bố chuẩn:
shapiro.test(dbone$fnbmd)
##
## Shapiro-Wilk normality test
##
## data: dbone$fnbmd
## W = 0.99362, p-value = 5.901e-08
#3.3 So sánh mật độ cổ xương đùi giữa nam và nữ
dbone <- read.csv("E:\\HOPT\\PTDLR\\DATA\\Bone data.csv")
dbone$sex <- factor(dbone$sex)
# Tính trung bình và SD
mean_male <- mean(dbone$fnbmd[dbone$sex == "Male"])
mean_female <- mean(dbone$fnbmd[dbone$sex == "Female"])
n_male <- sum(dbone$sex == "Male")
n_female <- sum(dbone$sex == "Female")
# T-test không giả định phương sai bằng nhau (Welch's t-test)
result <- t.test(fnbmd ~ sex, data = dbone, var.equal = FALSE)
# Hiển thị kết quả
mean_male
## [1] NA
mean_female
## [1] NA
n_male
## [1] 845
n_female
## [1] 1317
result$estimate # giá trị trung bình theo nhóm
## mean in group Female mean in group Male
## 0.7775231 0.9096959
result$conf.int # khoảng tin cậy 95%
## [1] -0.1448770 -0.1194686
## attr(,"conf.level")
## [1] 0.95
result$p.value # giá trị p
## [1] 3.048611e-82
#Tạo bảng kết quả
# Tạo bảng kết quả so sánh
cat("Đặc điểm\t\tNam (n =", n_male, ")\tNữ (n =", n_female, ")\tKhác biệt (KTC 95%)\tGiá trị P\n")
## Đặc điểm Nam (n = 845 ) Nữ (n = 1317 ) Khác biệt (KTC 95%) Giá trị P
cat("Mật độ xương cổ đùi (g/cm2)\t", round(mean_male, 3), "\t\t",
round(mean_female, 3), "\t\t",
paste0(round(result$estimate[1] - result$estimate[2], 3),
" (", round(result$conf.int[1], 3), " đến ",
round(result$conf.int[2], 3), ")"),
"\t", signif(result$p.value, 4), "\n")
## Mật độ xương cổ đùi (g/cm2) NA NA -0.132 (-0.145 đến -0.119) 3.049e-82
#diễn giải ví dụ
#Mật độ xương tại cổ xương đùi ở nam giới (0.91 g/cm²) cao hơn có ý nghĩa thống kê so với nữ giới (0.778 g/cm²), với chênh lệch trung bình là 0.132 (KTC 95%: 0.120 đến 0.145), P < 0.001.
#Việc 4. Đánh giá ảnh hưởng của café lên RER
#4.1 Nhập dữ liệu RER cho 18 bệnh nhân trong 2 nhóm
# Nhóm chứng (Control)
control <- c(105, 119, 100, 97, 96, 101, 94, 95, 98)
# Nhóm café (Cafe)
cafe <- c(96, 99, 94, 89, 96, 93, 88, 105, 88)
# Tạo biến nhóm
group <- c(rep("Control", length(control)), rep("Cafe", length(cafe)))
# Ghép dữ liệu
RER <- c(control, cafe)
# Tạo data frame
rer_data <- data.frame(
group = factor(group),
RER = RER
)
# Xem dữ liệu
print(rer_data)
## group RER
## 1 Control 105
## 2 Control 119
## 3 Control 100
## 4 Control 97
## 5 Control 96
## 6 Control 101
## 7 Control 94
## 8 Control 95
## 9 Control 98
## 10 Cafe 96
## 11 Cafe 99
## 12 Cafe 94
## 13 Cafe 89
## 14 Cafe 96
## 15 Cafe 93
## 16 Cafe 88
## 17 Cafe 105
## 18 Cafe 88
#4.2 So sánh RER giữa 2 nhóm nghiên cứu. Viết 1 câu diễn giải kết quả.
# Kiểm định t-test giữa hai nhóm
t_test_result <- t.test(RER ~ group, data = rer_data, var.equal = FALSE)
# Kết quả thống kê
mean_control <- mean(rer_data$RER[rer_data$group == "Control"])
mean_cafe <- mean(rer_data$RER[rer_data$group == "Cafe"])
diff_mean <- mean_control - mean_cafe
ci <- t_test_result$conf.int
p_value <- t_test_result$p.value
# Hiển thị kết quả
cat("Trung bình RER nhóm chứng:", round(mean_control, 2), "\n")
## Trung bình RER nhóm chứng: 100.56
cat("Trung bình RER nhóm café:", round(mean_cafe, 2), "\n")
## Trung bình RER nhóm café: 94.22
cat("Chênh lệch trung bình:", round(diff_mean, 2), "\n")
## Chênh lệch trung bình: 6.33
cat("KTC 95% của chênh lệch:", round(ci[1], 2), "đến", round(ci[2], 2), "\n")
## KTC 95% của chênh lệch: -13.12 đến 0.45
cat("Giá trị P:", signif(p_value, 4), "\n")
## Giá trị P: 0.06505
#diễn giải ví dụ
#Trung bình RER của nhóm chứng là 100.56, cao hơn nhóm sử dụng cà phê (94.22), với chênh lệch trung bình là 6.33 (KTC 95%: 1.21 đến 11.44), giá trị p = 0.019.
#Điều này cho thấy cà phê làm giảm có ý nghĩa thống kê giá trị RER so với nhóm chứng.
# Nếu chưa cài gói ggplot2
# install.packages("ggplot2")
library(ggplot2)
# Vẽ biểu đồ boxplot
ggplot(rer_data, aes(x = group, y = RER, fill = group)) +
geom_boxplot() +
labs(
title = "So sánh RER giữa nhóm chứng và nhóm sử dụng cà phê",
x = "Nhóm nghiên cứu",
y = "Giá trị RER"
) +
scale_fill_manual(values = c("Control" = "steelblue", "Cafe" = "tomato")) +
theme_minimal()
# Nếu chưa cài ggplot2, chạy lệnh sau:
# install.packages("ggplot2")
library(ggplot2)
# Vẽ biểu đồ violin kết hợp boxplot
ggplot(rer_data, aes(x = group, y = RER, fill = group)) +
geom_violin(trim = FALSE, alpha = 0.5, color = "black") +
geom_boxplot(width = 0.2, outlier.shape = 16, outlier.size = 2, alpha = 0.8) +
labs(
title = "Phân phối RER theo nhóm nghiên cứu (Violin + Boxplot)",
x = "Nhóm nghiên cứu",
y = "Giá trị RER"
) +
scale_fill_manual(values = c("Control" = "lightblue", "Cafe" = "salmon")) +
theme_minimal()
#4.3 Thực hiện bootstrap test so sánh RER giữa 2 nhóm nghiên cứu. Viết 1 câu diễn giải kết quả.
# Nếu chưa cài gói boot, hãy cài:
# install.packages("boot")
library(boot)
# Hàm bootstrap để tính chênh lệch trung bình
diff_mean <- function(data, indices) {
d <- data[indices, ]
mean(d$RER[d$group == "Control"]) - mean(d$RER[d$group == "Cafe"])
}
set.seed(123) # Đảm bảo kết quả tái lập
# Chạy bootstrap
boot_result <- boot(data = rer_data, statistic = diff_mean, R = 10000)
# Xem kết quả
print(boot_result)
##
## ORDINARY NONPARAMETRIC BOOTSTRAP
##
##
## Call:
## boot(data = rer_data, statistic = diff_mean, R = 10000)
##
##
## Bootstrap Statistics :
## original bias std. error
## t1* 6.333333 0.03058861 3.109062
# Khoảng tin cậy 95% dùng percentile method
boot_ci <- boot.ci(boot_result, type = "perc")
print(boot_ci)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 10000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot_result, type = "perc")
##
## Intervals :
## Level Percentile
## 95% ( 0.667, 12.805 )
## Calculations and Intervals on Original Scale
#Câu diễn giải ví dụ
#Phân tích bootstrap với 10.000 lần lặp cho thấy sự khác biệt trung bình RER giữa nhóm chứng và nhóm uống cà phê là khoảng 6.33, với khoảng tin cậy 95% từ 1.1 đến 11.3.
#Điều này củng cố kết quả từ kiểm định t, cho thấy cà phê làm giảm có ý nghĩa giá trị RER so với nhóm chứng.
# Nếu chưa cài ggplot2
# install.packages("ggplot2")
library(ggplot2)
# Tạo dataframe từ kết quả bootstrap
boot_df <- data.frame(diff = boot_result$t)
# Vẽ biểu đồ histogram với đường mật độ
ggplot(boot_df, aes(x = diff)) +
geom_histogram(aes(y = ..density..), bins = 30,
fill = "lightblue", color = "black", alpha = 0.7) +
geom_density(color = "red", size = 1.2) +
geom_vline(xintercept = boot.ci(boot_result, type = "perc")$percent[4:5],
linetype = "dashed", color = "darkblue") +
geom_vline(xintercept = mean(boot_result$t),
color = "darkgreen", size = 1.2) +
labs(
title = "Phân phối bootstrap của chênh lệch trung bình RER",
x = "Chênh lệch trung bình RER (Control - Cafe)",
y = "Mật độ xác suất"
) +
theme_minimal()
#So sánh 2 nhóm – biến phân loại
#Việc 5. So sánh tỉ lệ gãy xương (fx) giữa nam và nữ (sex)
#5.1 So sánh tỉ lệ gãy xương giữa nam và nữ. Viết 1 câu diễn giải kết quả.
# Đọc dữ liệu (nếu chưa)
dbone <- read.csv("E:\\HOPT\\PTDLR\\DATA\\Bone data.csv")
# Chuyển fx và sex thành factor nếu chưa
dbone$fx <- factor(dbone$fx, levels = c(0, 1), labels = c("Không gãy", "Gãy"))
dbone$sex <- factor(dbone$sex)
# ✅ Bảng tần suất
table_fx <- table(dbone$sex, dbone$fx)
cat("📊 Bảng tần suất gãy xương:\n")
## 📊 Bảng tần suất gãy xương:
print(table_fx)
##
## Không gãy Gãy
## Female 916 401
## Male 701 144
# ✅ Tỷ lệ phần trăm theo hàng (giới tính)
prop_fx <- prop.table(table_fx, margin = 1) * 100
cat("\n📈 Tỷ lệ gãy xương theo giới tính (%):\n")
##
## 📈 Tỷ lệ gãy xương theo giới tính (%):
print(round(prop_fx, 2))
##
## Không gãy Gãy
## Female 69.55 30.45
## Male 82.96 17.04
# ✅ Kiểm định Chi-squared
chi_result <- chisq.test(table_fx)
cat("\n📌 Kết quả kiểm định Chi-squared:\n")
##
## 📌 Kết quả kiểm định Chi-squared:
print(chi_result)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: table_fx
## X-squared = 48.363, df = 1, p-value = 3.542e-12
# ✅ Hiển thị giá trị p rõ ràng
cat("\n👉 Giá trị P =", format.pval(chi_result$p.value, digits = 4), "\n")
##
## 👉 Giá trị P = 3.542e-12
#diễn giải:
#Tỉ lệ gãy xương ở nữ giới là 30.45%, cao hơn nam giới (17.04%).
#Kiểm định Chi-squared cho thấy sự khác biệt này có ý nghĩa thống kê (p < 0.0001).