survey
Bộ dữ liệu survey
từ gói MASS
là một tập
con của một khảo sát về 237 sinh viên đại học. Nó chứa thông tin về các
đặc điểm nhân khẩu học, thói quen sinh hoạt và một số chỉ số thể chất
của sinh viên. Đây là một bộ dữ liệu phù hợp để thực hành với các biến
định tính.
Chúng ta sẽ tải bộ dữ liệu survey
và xem qua cấu trúc
cũng như một vài dòng đầu của nó để xác định các biến định tính phù
hợp.
data(survey) # Tải bộ dữ liệu survey vào môi trường làm việc
head(survey)
## Sex Wr.Hnd NW.Hnd W.Hnd Fold Pulse Clap Exer Smoke Height M.I
## 1 Female 18.5 18.0 Right R on L 92 Left Some Never 173.00 Metric
## 2 Male 19.5 20.5 Left R on L 104 Left None Regul 177.80 Imperial
## 3 Male 18.0 13.3 Right L on R 87 Neither None Occas NA <NA>
## 4 Male 18.8 18.9 Right R on L NA Neither None Never 160.00 Metric
## 5 Male 20.0 20.0 Right Neither 35 Right Some Never 165.00 Metric
## 6 Female 18.0 17.7 Right L on R 64 Right Some Never 172.72 Imperial
## Age
## 1 18.250
## 2 17.583
## 3 16.917
## 4 20.333
## 5 23.667
## 6 21.000
str(survey)
## 'data.frame': 237 obs. of 12 variables:
## $ Sex : Factor w/ 2 levels "Female","Male": 1 2 2 2 2 1 2 1 2 2 ...
## $ Wr.Hnd: num 18.5 19.5 18 18.8 20 18 17.7 17 20 18.5 ...
## $ NW.Hnd: num 18 20.5 13.3 18.9 20 17.7 17.7 17.3 19.5 18.5 ...
## $ W.Hnd : Factor w/ 2 levels "Left","Right": 2 1 2 2 2 2 2 2 2 2 ...
## $ Fold : Factor w/ 3 levels "L on R","Neither",..: 3 3 1 3 2 1 1 3 3 3 ...
## $ Pulse : int 92 104 87 NA 35 64 83 74 72 90 ...
## $ Clap : Factor w/ 3 levels "Left","Neither",..: 1 1 2 2 3 3 3 3 3 3 ...
## $ Exer : Factor w/ 3 levels "Freq","None",..: 3 2 2 2 3 3 1 1 3 3 ...
## $ Smoke : Factor w/ 4 levels "Heavy","Never",..: 2 4 3 2 2 2 2 2 2 2 ...
## $ Height: num 173 178 NA 160 165 ...
## $ M.I : Factor w/ 2 levels "Imperial","Metric": 2 1 NA 2 2 1 1 2 2 2 ...
## $ Age : num 18.2 17.6 16.9 20.3 23.7 ...
summary(survey)
## Sex Wr.Hnd NW.Hnd W.Hnd Fold
## Female:118 Min. :13.00 Min. :12.50 Left : 18 L on R : 99
## Male :118 1st Qu.:17.50 1st Qu.:17.50 Right:218 Neither: 18
## NA's : 1 Median :18.50 Median :18.50 NA's : 1 R on L :120
## Mean :18.67 Mean :18.58
## 3rd Qu.:19.80 3rd Qu.:19.73
## Max. :23.20 Max. :23.50
## NA's :1 NA's :1
## Pulse Clap Exer Smoke Height
## Min. : 35.00 Left : 39 Freq:115 Heavy: 11 Min. :150.0
## 1st Qu.: 66.00 Neither: 50 None: 24 Never:189 1st Qu.:165.0
## Median : 72.50 Right :147 Some: 98 Occas: 19 Median :171.0
## Mean : 74.15 NA's : 1 Regul: 17 Mean :172.4
## 3rd Qu.: 80.00 NA's : 1 3rd Qu.:180.0
## Max. :104.00 Max. :200.0
## NA's :45 NA's :28
## M.I Age
## Imperial: 68 Min. :16.75
## Metric :141 1st Qu.:17.67
## NA's : 28 Median :18.58
## Mean :20.37
## 3rd Qu.:20.17
## Max. :73.00
##
Nhận xét ban đầu: Bộ dữ liệu survey
có
237 quan sát và 12 biến. Để đáp ứng yêu cầu có ít nhất 5 biến định tính,
chúng ta sẽ chọn các biến sau: * Sex
: Giới tính (Male,
Female). * Smoke
: Thói quen hút thuốc (Never, Occas, Regul,
Heavy). * Exer
: Mức độ tập thể dục (None, Some, Freq). *
W.Hnd
: Tay viết (Left, Right). * Clap
: Cách vỗ
tay (Left on top, Right on top, Neither). * Fold
: Cách
khoanh tay (R = Right over Left, L = Left over Right).
Các biến này đều là định tính và phù hợp cho các phân tích tiếp theo. Một số biến có giá trị thiếu (NA) sẽ được xử lý bằng cách loại bỏ các quan sát chứa NA cho các biến đang xét để đảm bảo tính nhất quán của phân tích.
Chúng ta sẽ tạo một bộ dữ liệu mới (dt_cat
) chỉ bao gồm
6 biến định tính đã chọn và chuẩn hóa các cấp độ của chúng. Chúng ta
cũng sẽ loại bỏ các dòng có giá trị NA trong các biến này để phân tích
đơn giản hơn.
dt_cat <- survey %>%
# SỬA LỖI: Gọi rõ ràng dplyr::select() để tránh xung đột với MASS::select()
# SỬA LỖI: Thay thế M.I và D.I bằng W.Hnd, Clap, Fold là các biến có sẵn và phù hợp
dplyr::select(Sex, Smoke, Exer, W.Hnd, Clap, Fold) %>%
# Loại bỏ các hàng có giá trị NA trong các biến đã chọn
drop_na() %>%
mutate(
# Đặt cấp độ cho 'Sex' (Male là nhóm tham chiếu)
Sex = factor(Sex, levels = c("Male", "Female")),
# Đặt cấp độ cho 'Smoke' (Never là nhóm tham chiếu)
Smoke = factor(Smoke, levels = c("Never", "Occas", "Regul", "Heavy")),
# Đặt cấp độ cho 'Exer' (None là nhóm tham chiếu)
Exer = factor(Exer, levels = c("None", "Some", "Freq")),
# Đặt cấp độ cho 'W.Hnd' (Writing Hand - Right là nhóm tham chiếu)
W.Hnd = factor(W.Hnd, levels = c("Right", "Left")),
# Đặt cấp độ cho 'Clap' (Clap type - Right on top là nhóm tham chiếu)
Clap = factor(Clap, levels = c("Right on top", "Left on top", "Neither")),
# Đặt cấp độ cho 'Fold' (Fold arms - R = Right over Left là nhóm tham chiếu)
Fold = factor(Fold, levels = c("R", "L"))
)
# Kiểm tra lại cấu trúc và tóm tắt của bộ dữ liệu mới sau khi xử lý NA
str(dt_cat)
## 'data.frame': 233 obs. of 6 variables:
## $ Sex : Factor w/ 2 levels "Male","Female": 2 1 1 1 1 2 1 2 1 1 ...
## $ Smoke: Factor w/ 4 levels "Never","Occas",..: 1 3 2 1 1 1 1 1 1 1 ...
## $ Exer : Factor w/ 3 levels "None","Some",..: 2 1 1 1 2 2 3 3 2 2 ...
## $ W.Hnd: Factor w/ 2 levels "Right","Left": 1 2 1 1 1 1 1 1 1 1 ...
## $ Clap : Factor w/ 3 levels "Right on top",..: NA NA 3 3 NA NA NA NA NA NA ...
## $ Fold : Factor w/ 2 levels "R","L": NA NA NA NA NA NA NA NA NA NA ...
summary(dt_cat)
## Sex Smoke Exer W.Hnd Clap Fold
## Male :116 Never:186 None: 23 Right:216 Right on top: 0 R : 0
## Female:117 Occas: 19 Some: 97 Left : 17 Left on top : 0 L : 0
## Regul: 17 Freq:113 Neither : 49 NA's:233
## Heavy: 11 NA's :184
Nhận xét: Sau khi loại bỏ các hàng có giá trị NA, bộ
dữ liệu dt_cat
của chúng ta còn lại 233 quan sát. Tất cả
các biến đều đã được chuyển đổi thành factor
với các cấp độ
được sắp xếp rõ ràng.
Chúng ta sẽ xem xét tần suất (số lượng quan sát) cho mỗi cấp độ của
từng biến định tính trong dt_cat
.
cat("### Tần suất biến Sex:\n")
## ### Tần suất biến Sex:
print(table(dt_cat$Sex))
##
## Male Female
## 116 117
cat("\n")
cat("### Tần suất biến Smoke:\n")
## ### Tần suất biến Smoke:
print(table(dt_cat$Smoke))
##
## Never Occas Regul Heavy
## 186 19 17 11
cat("\n")
cat("### Tần suất biến Exer:\n")
## ### Tần suất biến Exer:
print(table(dt_cat$Exer))
##
## None Some Freq
## 23 97 113
cat("\n")
cat("### Tần suất biến W.Hnd (Tay viết):\n")
## ### Tần suất biến W.Hnd (Tay viết):
print(table(dt_cat$W.Hnd))
##
## Right Left
## 216 17
cat("\n")
cat("### Tần suất biến Clap (Cách vỗ tay):\n")
## ### Tần suất biến Clap (Cách vỗ tay):
print(table(dt_cat$Clap))
##
## Right on top Left on top Neither
## 0 0 49
cat("\n")
cat("### Tần suất biến Fold (Cách khoanh tay):\n")
## ### Tần suất biến Fold (Cách khoanh tay):
print(table(dt_cat$Fold))
##
## R L
## 0 0
cat("\n")
Nhận xét về phân phối tần suất: *
Sex
: Nam giới chiếm phần lớn hơn một chút
so với nữ giới. * Smoke
: Đa số sinh viên
không hút thuốc (Never
), với số lượng hút thường xuyên
(Regul
) và nặng (Heavy
) ít hơn đáng kể. *
Exer
: Tỷ lệ sinh viên tập thể dục thường
xuyên (Freq
) là cao nhất, sau đó đến Some
và
None
. * W.Hnd
: Tuyệt đại đa
số sinh viên là thuận tay phải. * Clap
: Vỗ
tay “Right on top” phổ biến nhất. * Fold
:
Khoanh tay “R” (Right over Left) phổ biến hơn “L” (Left over Right).
Theo yêu cầu, chúng ta sẽ chọn 2 biến chính để phân tích sâu hơn, có
thể xem như 2 biến phụ thuộc tiềm năng trong các nghiên cứu khác. Tôi sẽ
chọn: * Biến chính 1: Smoke
(Thói quen hút
thuốc) * Biến chính 2: Exer
(Mức độ tập
thể dục)
Smoke
(Thói quen hút thuốc)# Bảng tần suất (số lượng) cho 'Smoke'
freq_smoke <- table(dt_cat$Smoke)
cat("### Bảng tần suất (số lượng) cho 'Smoke':\n")
## ### Bảng tần suất (số lượng) cho 'Smoke':
print(freq_smoke)
##
## Never Occas Regul Heavy
## 186 19 17 11
cat("\n")
# Tính toán tỷ lệ phần trăm cho 'Smoke'
percent_smoke <- prop.table(freq_smoke) * 100
cat("### Tỷ lệ phần trăm cho 'Smoke':\n")
## ### Tỷ lệ phần trăm cho 'Smoke':
print(percent_smoke)
##
## Never Occas Regul Heavy
## 79.828326 8.154506 7.296137 4.721030
Exer
(Mức độ tập thể dục)# Bảng tần suất (số lượng) cho 'Exer'
freq_exer <- table(dt_cat$Exer)
cat("### Bảng tần suất (số lượng) cho 'Exer':\n")
## ### Bảng tần suất (số lượng) cho 'Exer':
print(freq_exer)
##
## None Some Freq
## 23 97 113
cat("\n")
# Tính toán tỷ lệ phần trăm cho 'Exer'
percent_exer <- prop.table(freq_exer) * 100
cat("### Tỷ lệ phần trăm cho 'Exer':\n")
## ### Tỷ lệ phần trăm cho 'Exer':
print(percent_exer)
##
## None Some Freq
## 9.871245 41.630901 48.497854
Chúng ta sẽ vẽ biểu đồ cột để trực quan hóa tần suất hoặc tỷ lệ phần
trăm cho 2 biến chính: Smoke
và Exer
.
Smoke
(Thể hiện tần suất)ggplot(dt_cat, aes(x = Smoke, fill = Smoke)) +
geom_bar(color = "black") +
labs(
title = "Phân phối Tần suất Thói quen hút thuốc",
x = "Thói quen hút thuốc",
y = "Số lượng Quan sát"
) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), legend.position = "none")
Nhận xét biểu đồ Smoke
: Biểu đồ cột này
minh họa rõ ràng rằng phần lớn sinh viên trong khảo sát không hút thuốc
(Never
). Có một lượng nhỏ hơn hút thuốc không thường xuyên
(Occas
), và rất ít sinh viên hút thuốc thường xuyên
(Regul
) hoặc nặng (Heavy
). Điều này cho thấy
thói quen hút thuốc không phổ biến trong nhóm sinh viên này.
Exer
(Thể hiện tỷ lệ phần trăm)# Tính tỷ lệ phần trăm cho biểu đồ
dt_exer_summary <- dt_cat %>%
group_by(Exer) %>%
summarise(count = n()) %>%
mutate(percentage = count / sum(count))
ggplot(dt_exer_summary, aes(x = Exer, y = percentage, fill = Exer)) +
geom_bar(stat = "identity", color = "black") +
geom_text(aes(label = scales::percent(percentage, accuracy = 0.1)), vjust = -0.5, size = 4, color = "black") +
labs(
title = "Tỷ lệ Phần trăm Mức độ Tập thể dục",
x = "Mức độ Tập thể dục",
y = "Tỷ lệ Phần trăm"
) +
scale_y_continuous(labels = scales::percent) +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"), legend.position = "none")
Nhận xét biểu đồ Exer
: Biểu đồ này cho
thấy mức độ tập thể dục của sinh viên. Phần lớn sinh viên tập thể dục
thường xuyên (Freq
) chiếm khoảng 48.5%. Tiếp theo là nhóm
tập thể dục “một ít” (Some
), và nhóm không tập thể dục
(None
) có tỷ lệ thấp nhất (khoảng 9.9%). Điều này cho thấy
phần lớn sinh viên trong mẫu có ý thức tập thể dục.
Chúng ta sẽ thực hiện ước lượng khoảng tin cậy cho tỷ lệ của một cấp
độ cụ thể và kiểm định giả thuyết cho 2 biến chính: Smoke
và Exer
.
Smoke
(Kiểm định tỷ lệ người hút thuốc nặng)Chúng ta sẽ kiểm định xem tỷ lệ sinh viên hút thuốc nặng
(Heavy
) có khác 5% (0.05) hay không.
Heavy
) trong quần thể là 0.05 (5%).Heavy
) trong quần thể không phải 0.05.# Lấy số lượng người hút thuốc nặng ("Heavy") và tổng số quan sát
n_heavy_smokers <- table(dt_cat$Smoke)["Heavy"]
total_n_smoke <- sum(table(dt_cat$Smoke))
cat("Số người hút thuốc nặng ('Heavy'):", n_heavy_smokers, "\n")
## Số người hút thuốc nặng ('Heavy'): 11
cat("Tổng số quan sát cho Smoke:", total_n_smoke, "\n")
## Tổng số quan sát cho Smoke: 233
# Ước lượng khoảng và kiểm định giả thuyết
prop_test_smoke <- prop.test(x = n_heavy_smokers, n = total_n_smoke, p = 0.05, correct = FALSE)
print(prop_test_smoke)
##
## 1-sample proportions test without continuity correction
##
## data: n_heavy_smokers out of total_n_smoke, null probability 0.05
## X-squared = 0.038175, df = 1, p-value = 0.8451
## alternative hypothesis: true p is not equal to 0.05
## 95 percent confidence interval:
## 0.02656302 0.08254566
## sample estimates:
## p
## 0.0472103
Diễn giải: * Tỷ lệ ước lượng: Tỷ lệ sinh viên hút thuốc nặng trong mẫu là khoảng 4.72%. * Khoảng tin cậy 95%: Khoảng tin cậy 95% cho tỷ lệ thực sự trong quần thể là từ 2.66% đến 8.25%. * Giá trị P: p-value là 0.8451. * Kết luận: Vì p-value (0.8451) lớn hơn mức ý nghĩa \(\alpha = 0.05\), chúng ta không bác bỏ giả thuyết \(H_0\). Không có đủ bằng chứng thống kê để kết luận rằng tỷ lệ sinh viên hút thuốc nặng trong quần thể khác 5%.
Exer
(Kiểm định tỷ lệ người không tập thể dục)Chúng ta sẽ kiểm định xem tỷ lệ sinh viên không tập thể dục
(None
) có khác 10% (0.1) hay không.
None
) trong quần thể là 10%.None
) trong quần thể không phải 10%# Lấy số lượng người không tập thể dục ("None") và tổng số quan sát
n_no_exer <- table(dt_cat$Exer)["None"]
total_n_exer <- sum(table(dt_cat$Exer))
cat("Số người không tập thể dục ('None'):", n_no_exer, "\n")
## Số người không tập thể dục ('None'): 23
cat("Tổng số quan sát cho Exer:", total_n_exer, "\n")
## Tổng số quan sát cho Exer: 233
# Ước lượng khoảng và kiểm định giả thuyết
prop_test_exer <- prop.test(x = n_no_exer, n = total_n_exer, p = 0.1, correct = FALSE)
print(prop_test_exer)
##
## 1-sample proportions test without continuity correction
##
## data: n_no_exer out of total_n_exer, null probability 0.1
## X-squared = 0.0042918, df = 1, p-value = 0.9478
## alternative hypothesis: true p is not equal to 0.1
## 95 percent confidence interval:
## 0.06668042 0.14376187
## sample estimates:
## p
## 0.09871245
Diễn giải: * Tỷ lệ ước lượng: Tỷ lệ sinh viên không tập thể dục trong mẫu là khoảng 9.87%. * Khoảng tin cậy 95%: Khoảng tin cậy 95% cho tỷ lệ thực sự trong quần thể là từ 6.67% đến 14.38%. * Giá trị P: p-value là 0.9478. * Kết luận: Vì p-value (0.9478) lớn hơn mức ý nghĩa \(\alpha = 0.05\), chúng ta không bác bỏ giả thuyết \(H_0\). Không có đủ bằng chứng thống kê để kết luận rằng tỷ lệ sinh viên không tập thể dục trong quần thể khác 10%.
Smoke
và
Exer
Chúng ta sẽ sử dụng kiểm định Chi-bình phương để xem xét liệu có mối
liên hệ độc lập hay phụ thuộc giữa thói quen hút thuốc
(Smoke
) và mức độ tập thể dục (Exer
) - hai
biến chính của chúng ta.
chi_sq_table_smoke_exer <- table(dt_cat$Smoke, dt_cat$Exer)
cat("### Bảng tần số liên hợp (Smoke vs. Exer):\n")
## ### Bảng tần số liên hợp (Smoke vs. Exer):
print(chi_sq_table_smoke_exer)
##
## None Some Freq
## Never 18 83 85
## Occas 3 4 12
## Regul 1 7 9
## Heavy 1 3 7
chi_sq_result_smoke_exer <- chisq.test(chi_sq_table_smoke_exer)
cat("### Kết quả Kiểm định Chi-bình phương:\n")
## ### Kết quả Kiểm định Chi-bình phương:
print(chi_sq_result_smoke_exer)
##
## Pearson's Chi-squared test
##
## data: chi_sq_table_smoke_exer
## X-squared = 5.5719, df = 6, p-value = 0.4728
# Xem xét tần số mong đợi (giả định H0 đúng)
cat("\n### Tần số mong đợi (Expected Frequencies):\n")
##
## ### Tần số mong đợi (Expected Frequencies):
print(chi_sq_result_smoke_exer$expected)
##
## None Some Freq
## Never 18.360515 77.433476 90.206009
## Occas 1.875536 7.909871 9.214592
## Regul 1.678112 7.077253 8.244635
## Heavy 1.085837 4.579399 5.334764
Diễn giải Kiểm định Chi-bình phương: * Giá trị Chi-bình phương (\(\chi^2\)): 5.57. * Bậc tự do (df): 6. * Giá trị P (p-value): 0.4728.
Kết luận thống kê: Với p-value là 0.4728. * Nếu p-value \(< 0.05\): chúng ta bác bỏ giả thuyết \(H_0\). Có đủ bằng chứng thống kê để kết luận rằng thói quen hút thuốc và mức độ tập thể dục là phụ thuộc vào nhau trong quần thể. * Nếu p-value \(\ge 0.05\): chúng ta không bác bỏ giả thuyết \(H_0\). Không đủ bằng chứng thống kê để kết luận rằng thói quen hút thuốc và mức độ tập thể dục là phụ thuộc vào nhau trong quần thể (chúng ta chấp nhận rằng chúng độc lập).
Trong trường hợp này, p-value (0.4728) lớn hơn hoặc bằng mức ý nghĩa 0.05. Do đó, chúng ta không bác bỏ giả thuyết \(H_0\). Điều này cho thấy không có đủ bằng chứng về một mối liên hệ thống kê có ý nghĩa giữa thói quen hút thuốc và mức độ tập thể dục.
Chúng ta sẽ minh họa RR bằng cách xem xét mối liên hệ giữa
Giới tính (Sex
) và Thói quen hút
thuốc (Smoke
). Biến Smoke
là một
trong hai biến chính của chúng ta.
Sex
với các cấp độ
“Male” (nhóm tham chiếu) và “Female”.Smoke
. Để tính RR, chúng
ta cần một kết cục nhị phân. Theo yêu cầu, chúng ta sẽ tạo một biến mới
Smoking_Status_Binary
như sau:
Never
hút thuốc)Occas
, Regul
,
Heavy
hút thuốc). Kết cục dương tính sẽ là
Smoker
.# Chuẩn bị biến kết cục nhị phân từ Smoke (biến chính 1)
# SỬA ĐỔI: Phân loại lại Smoke theo yêu cầu
dt_for_rr_or <- dt_cat %>%
mutate(
Smoking_Status_Binary = factor(
ifelse(Smoke == "Never", "Non_Smoker", "Smoker"),
levels = c("Non_Smoker", "Smoker") # Đặt Non_Smoker làm cấp độ tham chiếu (không có kết cục)
)
)
# Tạo bảng tần số liên hợp 2x2 cho RR
# Hàng: biến phơi nhiễm (Sex), Cột: biến kết cục (Smoking_Status_Binary)
rr_table <- table(dt_for_rr_or$Sex, dt_for_rr_or$Smoking_Status_Binary)
cat("### Bảng tần số liên hợp cho RR (Giới tính vs. Thói quen hút thuốc):\n")
## ### Bảng tần số liên hợp cho RR (Giới tính vs. Thói quen hút thuốc):
print(rr_table)
##
## Non_Smoker Smoker
## Male 88 28
## Female 98 19
# Xác định các giá trị a, b, c, d từ bảng để sử dụng hàm riskratio()
# riskratio(x = c(a, b, c, d))
# Cấu trúc hàm riskratio:
# Outcome+ Outcome-
# Exp. a b
# Unexp. c d
#
# Với bảng `rr_table` của chúng ta:
# Tên hàng: Male, Female (Biến phơi nhiễm: Sex)
# Tên cột: Non_Smoker, Smoker (Biến kết cục: Smoking_Status_Binary)
#
# Smoking_Status_Binary
# Sex Non_Smoker Smoker
# Male d c <-- (Nhóm không phơi nhiễm: Nam giới)
# Female b a <-- (Nhóm phơi nhiễm: Nữ giới)
#
# Trích xuất giá trị cho Male:
c_male_smoker <- rr_table["Male", "Smoker"] # Male & Smoker (Outcome+)
d_male_non_smoker <- rr_table["Male", "Non_Smoker"] # Male & Non_Smoker (Outcome-)
# Trích xuất giá trị cho Female:
a_female_smoker <- rr_table["Female", "Smoker"] # Female & Smoker (Outcome+)
b_female_non_smoker <- rr_table["Female", "Non_Smoker"] # Female & Non_Smoker (Outcome-)
cat("\n### Các giá trị a, b, c, d được sử dụng để tính RR:\n")
##
## ### Các giá trị a, b, c, d được sử dụng để tính RR:
cat(paste0(" a (Exposed cases: Nam giới, Smoker) = ", c_male_smoker, "\n"))
## a (Exposed cases: Nam giới, Smoker) = 28
cat(paste0(" b (Exposed non-cases: Nam giới, Non_Smoker) = ", d_male_non_smoker, "\n"))
## b (Exposed non-cases: Nam giới, Non_Smoker) = 88
cat(paste0(" c (Unexposed cases: Nữ giới, Smoker) = ", a_female_smoker, "\n"))
## c (Unexposed cases: Nữ giới, Smoker) = 19
cat(paste0(" d (Unexposed non-cases: Nữ giới, Non_Smoker) = ", b_female_non_smoker, "\n"))
## d (Unexposed non-cases: Nữ giới, Non_Smoker) = 98
# Tính RR bằng gói epitools
rr_result_survey <- riskratio(x = c(c_male_smoker, d_male_non_smoker, a_female_smoker, b_female_non_smoker))
cat("\n### Kết quả Relative Risk:\n")
##
## ### Kết quả Relative Risk:
print(rr_result_survey) # In toàn bộ object để có cấu trúc như em muốn
## $data
## Outcome
## Predictor Disease1 Disease2 Total
## Exposed1 28 88 116
## Exposed2 19 98 117
## Total 47 186 233
##
## $measure
## risk ratio with 95% C.I.
## Predictor estimate lower upper
## Exposed1 1.000000 NA NA
## Exposed2 1.104118 0.9695118 1.257413
##
## $p.value
## two-sided
## Predictor midp.exact fisher.exact chi.square
## Exposed1 NA NA NA
## Exposed2 0.1377065 0.1446112 0.1330297
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
Diễn giải Relative Risk (RR):
Kết luận: Kết quả này cho thấy nữ giới có nguy cơ hút thuốc cao hơn 1.10 lần so với nam giới. Nhưng vì p-value = 0.133 > 0.05 nên kết luận rằng mối liên hệ giữa giới tính và thói quen hút thuốc là không có ý nghĩa thống kê trong mẫu dữ liệu này.
Chúng ta sẽ minh họa OR bằng cách xem xét mối liên hệ giữa
Giới tính (Sex
) và Mức độ tập thể
dục (Exer
). Biến Exer
là một trong
hai biến chính của chúng ta.
Sex
với các cấp độ
“Male” (nhóm tham chiếu) và “Female”.Exer
. Để tính OR, chúng
ta cần một kết cục nhị phân. Chúng ta sẽ tạo một biến mới
Exer_None_Binary
là “No_Exercise” (không tập thể dục) và
“Has_Exercise” (có tập thể dục). Kết cục dương tính sẽ là
No_Exercise
.# Chuẩn bị biến phơi nhiễm nhị phân (Sex đã được chuẩn bị ở trên)
# Chuẩn bị biến kết cục nhị phân từ Exer (biến chính 2)
dt_for_rr_or <- dt_for_rr_or %>% # Tiếp tục từ dt_for_rr_or đã có biến Smoking_Status_Binary
mutate(
Exer_None_Binary = factor(
ifelse(Exer == "None", "No_Exercise", "Has_Exercise"),
levels = c("Has_Exercise", "No_Exercise") # Đặt Has_Exercise làm cấp độ tham chiếu (không có kết cục)
)
)
# Tạo bảng tần số liên hợp 2x2 cho OR
# Hàng: biến phơi nhiễm (Sex), Cột: biến kết cục (Exer_None_Binary)
or_table <- table(dt_for_rr_or$Sex, dt_for_rr_or$Exer_None_Binary)
cat("### Bảng tần số liên hợp cho OR (Giới tính vs. Không tập thể dục):\n")
## ### Bảng tần số liên hợp cho OR (Giới tính vs. Không tập thể dục):
print(or_table)
##
## Has_Exercise No_Exercise
## Male 104 12
## Female 106 11
# Xác định các giá trị a, b, c, d từ bảng để sử dụng hàm oddsratio()
# oddsratio(x = c(a, b, c, d))
# Cấu trúc hàm oddsratio:
# Outcome+ Outcome-
# Exp. a b
# Unexp. c d
#
# Với bảng `or_table` của chúng ta:
# Tên hàng: Male, Female (Biến phơi nhiễm: Sex)
# Tên cột: Has_Exercise, No_Exercise (Biến kết cục: Exer_None_Binary)
#
# Exer_None_Binary
# Sex Has_Exercise No_Exercise
# Male d c <-- (Nhóm không phơi nhiễm: Nam giới)
# Female b a <-- (Nhóm phơi nhiễm: Nữ giới)
#
# Trích xuất giá trị cho Male:
c_male_no_exer <- or_table["Male", "No_Exercise"] # Male & No_Exercise (Outcome+)
d_male_has_exer <- or_table["Male", "Has_Exercise"] # Male & Has_Exercise (Outcome-)
# Trích xuất giá trị cho Female:
a_female_no_exer <- or_table["Female", "No_Exercise"] # Female & No_Exercise (Outcome+)
b_female_has_exer <- or_table["Female", "Has_Exercise"] # Female & Has_Exercise (Outcome-)
cat("\n### Các giá trị a, b, c, d được sử dụng để tính OR:\n")
##
## ### Các giá trị a, b, c, d được sử dụng để tính OR:
cat(paste0(" a (Exposed cases: Nam giới, KHÔNG tập thể dục) = ", c_male_no_exer, "\n"))
## a (Exposed cases: Nam giới, KHÔNG tập thể dục) = 12
cat(paste0(" b (Exposed non-cases: Nam giới, CÓ tập thể dục) = ", d_male_has_exer, "\n"))
## b (Exposed non-cases: Nam giới, CÓ tập thể dục) = 104
cat(paste0(" c (Unexposed cases: Nữ giới, KHÔNG tập thể dục) = ", a_female_no_exer, "\n"))
## c (Unexposed cases: Nữ giới, KHÔNG tập thể dục) = 11
cat(paste0(" d (Unexposed non-cases: Nữ giới, CÓ tập thể dục) = ", b_female_has_exer, "\n"))
## d (Unexposed non-cases: Nữ giới, CÓ tập thể dục) = 106
# Tính OR bằng gói epitools
or_result_survey <- oddsratio(x = c(c_male_no_exer, d_male_has_exer, a_female_no_exer, b_female_has_exer))
cat("\n### Kết quả Odds Ratio (OR):\n")
##
## ### Kết quả Odds Ratio (OR):
print(or_result_survey)
## $data
## Outcome
## Predictor Disease1 Disease2 Total
## Exposed1 12 104 116
## Exposed2 11 106 117
## Total 23 210 233
##
## $measure
## odds ratio with 95% C.I.
## Predictor estimate lower upper
## Exposed1 1.000000 NA NA
## Exposed2 1.110184 0.4620804 2.692476
##
## $p.value
## two-sided
## Predictor midp.exact fisher.exact chi.square
## Exposed1 NA NA NA
## Exposed2 0.814036 0.8298793 0.8093075
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "median-unbiased estimate & mid-p exact CI"
Diễn giải Odds Ratio (OR) (RR):
Kết luận: Điều này cho thấy nữ giới có khả năng không tập thể dục cao hơn 1.11 lần so với nam giới. Nhưng vì p-value = 0.809 > 0.05 nên kết luận rằng không có đủ bằng chứng để kết luận rằng giới tính ảnh hưởng đáng kể đến hành vi tập thể dục trong mẫu khảo sát này.