Phân tích dữ liệu sống còn và mô hình dự báo tái phạm bằng R
Thien Duong
2026-08-18
1 1. Phân tích dữ liệu sống còn – So sánh tuổi thọ bóng đèn
1.1 1.1 Khởi tạo dữ liệu
dat <- data.frame(
id = 1:30,
group = c(
rep("Thường", 15),
rep("LED", 15)
),
voltage = c(
213.3, 231.1, 210.0, 220.6, 233.7,
215.2, 216.2, 214.9, 217.7, 221.1,
229.8, 213.6, 211.4, 218.7, 211.4,
218.9, 215.2, 202.5, 221.9, 217.9,
227.2, 227.5, 231.7, 225.7, 226.6,
217.7, 231.3, 232.0, 214.7, 213.2
),
time = c(
1060, 1615, 380, 17, 11,
233, 230, 2000, 130, 289,
184, 735, 2000, 308, 2000,
352, 2000, 2000, 2000, 2000,
172, 527, 2000, 2000, 2000,
2000, 66, 530, 2000, 2000
),
status = c(
1, 1, 1, 1, 1,
1, 1, 0, 1, 1,
1, 1, 0, 1, 0,
1, 0, 0, 0, 0,
1, 1, 0, 0, 0,
0, 1, 1, 0, 0
)
)
dim(dat)## [1] 30 5head(dat)## id group voltage time status
## 1 1 Thường 213.3 1060 1
## 2 2 Thường 231.1 1615 1
## 3 3 Thường 210.0 380 1
## 4 4 Thường 220.6 17 1
## 5 5 Thường 233.7 11 1
## 6 6 Thường 215.2 233 1Dữ liệu gồm 30 bóng đèn thuộc hai nhóm là bóng đèn thường và bóng đèn
LED.
Biến time thể hiện thời gian hoạt động và
status cho biết bóng đèn đã hỏng hay bị kiểm duyệt.
1.2 1.2 Kaplan-Meier
library(survival)
library(survminer)
km <- survfit(
Surv(time, status) ~ group,
data = dat
)
summary(km)## Call: survfit(formula = Surv(time, status) ~ group, data = dat)
##
## group=LED
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 66 15 1 0.933 0.0644 0.815 1.000
## 172 14 1 0.867 0.0878 0.711 1.000
## 352 13 1 0.800 0.1033 0.621 1.000
## 527 12 1 0.733 0.1142 0.540 0.995
## 530 11 1 0.667 0.1217 0.466 0.953
##
## group=Thường
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 11 15 1 0.933 0.0644 0.8153 1.000
## 17 14 1 0.867 0.0878 0.7106 1.000
## 130 13 1 0.800 0.1033 0.6212 1.000
## 184 12 1 0.733 0.1142 0.5405 0.995
## 230 11 1 0.667 0.1217 0.4661 0.953
## 233 10 1 0.600 0.1265 0.3969 0.907
## 289 9 1 0.533 0.1288 0.3322 0.856
## 308 8 1 0.467 0.1288 0.2717 0.802
## 380 7 1 0.400 0.1265 0.2152 0.743
## 735 6 1 0.333 0.1217 0.1630 0.682
## 1060 5 1 0.267 0.1142 0.1152 0.617
## 1615 4 1 0.200 0.1033 0.0727 0.550ggsurvplot(
km,
data = dat,
risk.table = TRUE,
pval = TRUE,
xlab = "Thời gian",
ylab = "Xác suất sống còn",
title = "Kaplan-Meier theo loại bóng đèn",
ggtheme = theme_minimal()
)
1.2.1 Kiểm định Log-rank
survdiff(
Surv(time, status) ~ group,
data = dat
)## Call:
## survdiff(formula = Surv(time, status) ~ group, data = dat)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## group=LED 15 5 10.05 2.54 6.31
## group=Thường 15 12 6.95 3.67 6.31
##
## Chisq= 6.3 on 1 degrees of freedom, p= 0.01Kiểm định Log-rank được sử dụng để đánh giá sự khác biệt giữa hai đường sống còn.
1.3 1.3 Hồi quy Cox
modelCox_lightBulb <- coxph(
Surv(time, status) ~ group,
data = dat
)
summary(modelCox_lightBulb)## Call:
## coxph(formula = Surv(time, status) ~ group, data = dat)
##
## n= 30, number of events= 17
##
## coef exp(coef) se(coef) z Pr(>|z|)
## groupThường 1.2655 3.5448 0.5364 2.359 0.0183 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## groupThường 3.545 0.2821 1.239 10.14
##
## Concordance= 0.65 (se = 0.059 )
## Likelihood ratio test= 6.23 on 1 df, p=0.01
## Wald test = 5.57 on 1 df, p=0.02
## Score (logrank) test = 6.31 on 1 df, p=0.01Mô hình Cox được sử dụng để đánh giá mối liên quan giữa loại bóng đèn và nguy cơ bóng đèn bị hỏng theo thời gian.
2 2. Hồi quy Cox – Đánh giá mối liên quan đến tái phạm
2.1 2.1 Đọc dữ liệu Rossi
library(carData)
data("Rossi")
names(Rossi)## [1] "week" "arrest" "fin" "age" "race" "wexp" "mar" "paro"
## [9] "prio" "educ" "emp1" "emp2" "emp3" "emp4" "emp5" "emp6"
## [17] "emp7" "emp8" "emp9" "emp10" "emp11" "emp12" "emp13" "emp14"
## [25] "emp15" "emp16" "emp17" "emp18" "emp19" "emp20" "emp21" "emp22"
## [33] "emp23" "emp24" "emp25" "emp26" "emp27" "emp28" "emp29" "emp30"
## [41] "emp31" "emp32" "emp33" "emp34" "emp35" "emp36" "emp37" "emp38"
## [49] "emp39" "emp40" "emp41" "emp42" "emp43" "emp44" "emp45" "emp46"
## [57] "emp47" "emp48" "emp49" "emp50" "emp51" "emp52"dat_Rossi <- Rossi[
c(
"week",
"arrest",
"fin",
"age",
"race",
"wexp",
"mar",
"paro",
"prio",
"educ"
)
]
head(dat_Rossi)## week arrest fin age race wexp mar paro prio educ
## 1 20 1 no 27 black no not married yes 3 3
## 2 17 1 no 18 black no not married yes 8 4
## 3 25 1 no 19 other yes not married yes 13 3
## 4 52 0 yes 23 black yes married yes 1 5
## 5 52 0 no 19 other yes not married yes 3 3
## 6 52 0 no 24 black yes not married no 2 42.2 2.2 Kích thước dữ liệu
nrow(dat_Rossi)## [1] 432ncol(dat_Rossi)## [1] 102.3 2.3 Mô tả dữ liệu
library(table1)
table1(
~ arrest + week + fin + age + race +
wexp + mar + paro + prio + educ,
data = dat_Rossi
)| Overall (N=432) |
|
|---|---|
| arrest | |
| Mean (SD) | 0.264 (0.441) |
| Median [Min, Max] | 0 [0, 1.00] |
| week | |
| Mean (SD) | 45.9 (12.7) |
| Median [Min, Max] | 52.0 [1.00, 52.0] |
| fin | |
| no | 216 (50.0%) |
| yes | 216 (50.0%) |
| age | |
| Mean (SD) | 24.6 (6.11) |
| Median [Min, Max] | 23.0 [17.0, 44.0] |
| race | |
| black | 379 (87.7%) |
| other | 53 (12.3%) |
| wexp | |
| no | 185 (42.8%) |
| yes | 247 (57.2%) |
| mar | |
| married | 53 (12.3%) |
| not married | 379 (87.7%) |
| paro | |
| no | 165 (38.2%) |
| yes | 267 (61.8%) |
| prio | |
| Mean (SD) | 2.98 (2.90) |
| Median [Min, Max] | 2.00 [0, 18.0] |
| educ | |
| Mean (SD) | 3.48 (0.834) |
| Median [Min, Max] | 3.00 [2.00, 6.00] |
2.4 2.4 Kaplan-Meier của toàn bộ mẫu nghiên cứu
dat_Rossi <- na.omit(dat_Rossi)
km_Rossi <- survfit(
Surv(week, arrest) ~ 1,
data = dat_Rossi
)
summary(km_Rossi)## Call: survfit(formula = Surv(week, arrest) ~ 1, data = dat_Rossi)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 432 1 0.998 0.00231 0.993 1.000
## 2 431 1 0.995 0.00327 0.989 1.000
## 3 430 1 0.993 0.00400 0.985 1.000
## 4 429 1 0.991 0.00461 0.982 1.000
## 5 428 1 0.988 0.00515 0.978 0.999
## 6 427 1 0.986 0.00563 0.975 0.997
## 7 426 1 0.984 0.00607 0.972 0.996
## 8 425 5 0.972 0.00791 0.957 0.988
## 9 420 2 0.968 0.00852 0.951 0.984
## 10 418 1 0.965 0.00881 0.948 0.983
## 11 417 2 0.961 0.00935 0.942 0.979
## 12 415 2 0.956 0.00987 0.937 0.976
## 13 413 1 0.954 0.01011 0.934 0.974
## 14 412 3 0.947 0.01080 0.926 0.968
## 15 409 2 0.942 0.01123 0.920 0.964
## 16 407 2 0.937 0.01165 0.915 0.961
## 17 405 3 0.931 0.01223 0.907 0.955
## 18 402 3 0.924 0.01278 0.899 0.949
## 19 399 2 0.919 0.01313 0.894 0.945
## 20 397 5 0.907 0.01395 0.880 0.935
## 21 392 2 0.903 0.01425 0.875 0.931
## 22 390 1 0.900 0.01440 0.873 0.929
## 23 389 1 0.898 0.01455 0.870 0.927
## 24 388 4 0.889 0.01512 0.860 0.919
## 25 384 3 0.882 0.01552 0.852 0.913
## 26 381 3 0.875 0.01591 0.844 0.907
## 27 378 2 0.870 0.01616 0.839 0.903
## 28 376 2 0.866 0.01640 0.834 0.898
## 30 374 2 0.861 0.01664 0.829 0.894
## 31 372 1 0.859 0.01675 0.827 0.892
## 32 371 2 0.854 0.01698 0.822 0.888
## 33 369 2 0.850 0.01720 0.816 0.884
## 34 367 2 0.845 0.01742 0.811 0.880
## 35 365 4 0.836 0.01783 0.801 0.871
## 36 361 3 0.829 0.01813 0.794 0.865
## 37 358 4 0.819 0.01851 0.784 0.857
## 38 354 1 0.817 0.01860 0.781 0.854
## 39 353 2 0.812 0.01878 0.777 0.850
## 40 351 4 0.803 0.01913 0.767 0.842
## 42 347 2 0.799 0.01929 0.762 0.837
## 43 345 4 0.789 0.01962 0.752 0.829
## 44 341 2 0.785 0.01977 0.747 0.824
## 45 339 2 0.780 0.01993 0.742 0.820
## 46 337 4 0.771 0.02022 0.732 0.812
## 47 333 1 0.769 0.02029 0.730 0.809
## 48 332 2 0.764 0.02043 0.725 0.805
## 49 330 5 0.752 0.02077 0.713 0.794
## 50 325 3 0.745 0.02096 0.705 0.788
## 52 322 4 0.736 0.02121 0.696 0.779ggsurvplot(
km_Rossi,
data = dat_Rossi,
risk.table = TRUE,
risk.table.height = 0.22,
risk.table.y.text.col = TRUE,
risk.table.y.text = FALSE,
conf.int = TRUE,
censor = TRUE,
xlab = "Thời gian theo dõi (tuần)",
ylab = "Xác suất chưa bị bắt lại",
title = "Kaplan-Meier của mẫu nghiên cứu",
break.time.by = 10,
xlim = c(0, 52),
ylim = c(0, 1),
ggtheme = theme_minimal(),
font.title = c(16, "bold"),
font.x = 13,
font.y = 13,
font.tickslab = 11
)
2.5 2.5 Kaplan-Meier theo hỗ trợ tài chính
km_Rossi_finance <- survfit(
Surv(week, arrest) ~ fin,
data = dat_Rossi
)
summary(km_Rossi_finance)## Call: survfit(formula = Surv(week, arrest) ~ fin, data = dat_Rossi)
##
## fin=no
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 216 1 0.995 0.00462 0.986 1.000
## 2 215 1 0.991 0.00652 0.978 1.000
## 3 214 1 0.986 0.00796 0.971 1.000
## 4 213 1 0.981 0.00917 0.964 1.000
## 5 212 1 0.977 0.01023 0.957 0.997
## 6 211 1 0.972 0.01118 0.951 0.994
## 8 210 2 0.963 0.01285 0.938 0.988
## 10 208 1 0.958 0.01360 0.932 0.985
## 11 207 1 0.954 0.01430 0.926 0.982
## 13 206 1 0.949 0.01496 0.920 0.979
## 14 205 2 0.940 0.01618 0.909 0.972
## 16 203 1 0.935 0.01675 0.903 0.969
## 17 202 3 0.921 0.01832 0.886 0.958
## 18 199 2 0.912 0.01927 0.875 0.951
## 19 197 1 0.907 0.01972 0.870 0.947
## 20 196 1 0.903 0.02016 0.864 0.943
## 21 195 2 0.894 0.02099 0.853 0.936
## 23 193 1 0.889 0.02138 0.848 0.932
## 24 192 2 0.880 0.02214 0.837 0.924
## 25 190 3 0.866 0.02320 0.821 0.912
## 26 187 3 0.852 0.02417 0.806 0.901
## 27 184 2 0.843 0.02478 0.795 0.893
## 28 182 2 0.833 0.02536 0.785 0.885
## 30 180 1 0.829 0.02564 0.780 0.881
## 31 179 1 0.824 0.02591 0.775 0.876
## 32 178 1 0.819 0.02617 0.770 0.872
## 33 177 2 0.810 0.02668 0.760 0.864
## 35 175 1 0.806 0.02693 0.754 0.860
## 37 174 2 0.796 0.02740 0.744 0.852
## 38 172 1 0.792 0.02763 0.739 0.848
## 39 171 1 0.787 0.02786 0.734 0.844
## 40 170 3 0.773 0.02850 0.719 0.831
## 43 167 4 0.755 0.02928 0.699 0.814
## 44 163 2 0.745 0.02964 0.689 0.806
## 45 161 1 0.741 0.02982 0.685 0.802
## 46 160 1 0.736 0.02999 0.680 0.797
## 47 159 1 0.731 0.03016 0.675 0.793
## 48 158 1 0.727 0.03032 0.670 0.789
## 49 157 2 0.718 0.03063 0.660 0.780
## 50 155 1 0.713 0.03078 0.655 0.776
## 52 154 4 0.694 0.03134 0.636 0.759
##
## fin=yes
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 7 216 1 0.995 0.00462 0.986 1.000
## 8 215 3 0.981 0.00917 0.964 1.000
## 9 212 2 0.972 0.01118 0.951 0.994
## 11 210 1 0.968 0.01205 0.944 0.991
## 12 209 2 0.958 0.01360 0.932 0.985
## 14 207 1 0.954 0.01430 0.926 0.982
## 15 206 2 0.944 0.01559 0.914 0.975
## 16 204 1 0.940 0.01618 0.909 0.972
## 18 203 1 0.935 0.01675 0.903 0.969
## 19 202 1 0.931 0.01730 0.897 0.965
## 20 201 4 0.912 0.01927 0.875 0.951
## 22 197 1 0.907 0.01972 0.870 0.947
## 24 196 2 0.898 0.02058 0.859 0.939
## 30 194 1 0.894 0.02099 0.853 0.936
## 32 193 1 0.889 0.02138 0.848 0.932
## 34 192 2 0.880 0.02214 0.837 0.924
## 35 190 3 0.866 0.02320 0.821 0.912
## 36 187 3 0.852 0.02417 0.806 0.901
## 37 184 2 0.843 0.02478 0.795 0.893
## 39 182 1 0.838 0.02507 0.790 0.889
## 40 181 1 0.833 0.02536 0.785 0.885
## 42 180 2 0.824 0.02591 0.775 0.876
## 45 178 1 0.819 0.02617 0.770 0.872
## 46 177 3 0.806 0.02693 0.754 0.860
## 48 174 1 0.801 0.02717 0.749 0.856
## 49 173 3 0.787 0.02786 0.734 0.844
## 50 170 2 0.778 0.02829 0.724 0.835ggsurvplot(
km_Rossi_finance,
data = dat_Rossi,
risk.table = TRUE,
risk.table.height = 0.25,
risk.table.y.text.col = TRUE,
risk.table.y.text = TRUE,
pval = TRUE,
pval.coord = c(3, 0.18),
conf.int = TRUE,
censor = TRUE,
xlab = "Thời gian theo dõi (tuần)",
ylab = "Xác suất chưa bị bắt lại",
title = "Kaplan-Meier theo tình trạng hỗ trợ tài chính",
legend.title = "Hỗ trợ tài chính",
legend.labs = c("Không hỗ trợ", "Có hỗ trợ"),
break.time.by = 10,
xlim = c(0, 52),
ylim = c(0, 1),
ggtheme = theme_minimal(),
font.title = c(16, "bold"),
font.x = 13,
font.y = 13,
font.tickslab = 11
)
2.6 2.6 Xác suất bị bắt lại trong 10 tuần đầu
km_Rossi_fin_10 <- survfit(
Surv(week, arrest) ~ fin,
data = dat_Rossi
)
summary_10 <- summary(
km_Rossi_fin_10,
times = 10
)
km_rossi_prob <- data.frame(
group = summary_10$strata,
survival_prob = summary_10$surv,
arrest_prob = 1 - summary_10$surv
)
km_rossi_prob$arrest_percent <-
km_rossi_prob$arrest_prob * 100
km_rossi_prob## group survival_prob arrest_prob arrest_percent
## 1 fin=no 0.9583333 0.04166667 4.166667
## 2 fin=yes 0.9722222 0.02777778 2.777778Xác suất bị bắt lại được tính theo:
\[ P(T \leq 10) = 1 - S(10) \]
2.7 2.7 Mô hình Cox với hỗ trợ tài chính
modelCox_RossiFin <- coxph(
Surv(week, arrest) ~ fin,
data = dat_Rossi
)
summary(modelCox_RossiFin)## Call:
## coxph(formula = Surv(week, arrest) ~ fin, data = dat_Rossi)
##
## n= 432, number of events= 114
##
## coef exp(coef) se(coef) z Pr(>|z|)
## finyes -0.3691 0.6914 0.1897 -1.945 0.0517 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## finyes 0.6914 1.446 0.4767 1.003
##
## Concordance= 0.546 (se = 0.023 )
## Likelihood ratio test= 3.84 on 1 df, p=0.05
## Wald test = 3.78 on 1 df, p=0.05
## Score (logrank) test = 3.83 on 1 df, p=0.05exp(coef(modelCox_RossiFin))## finyes
## 0.6913776exp(confint(modelCox_RossiFin))## 2.5 % 97.5 %
## finyes 0.476676 1.002784exp(coef) là Hazard Ratio (HR). HR nhỏ hơn 1 cho thấy
yếu tố liên quan đến giảm hazard bị bắt lại, trong khi HR lớn hơn 1 cho
thấy hazard tăng.
2.8 2.8 Cox hiệu chỉnh yếu tố nhiễu
modelCox_adjusted <- coxph(
Surv(week, arrest) ~
fin + age + race + prio + educ,
data = dat_Rossi
)
summary(modelCox_adjusted)## Call:
## coxph(formula = Surv(week, arrest) ~ fin + age + race + prio +
## educ, data = dat_Rossi)
##
## n= 432, number of events= 114
##
## coef exp(coef) se(coef) z Pr(>|z|)
## finyes -0.34523 0.70806 0.19107 -1.807 0.07079 .
## age -0.06625 0.93590 0.02072 -3.198 0.00138 **
## raceother -0.39697 0.67235 0.30763 -1.290 0.19690
## prio 0.09119 1.09548 0.02794 3.263 0.00110 **
## educ -0.20090 0.81800 0.13106 -1.533 0.12530
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## finyes 0.7081 1.4123 0.4869 1.0297
## age 0.9359 1.0685 0.8987 0.9747
## raceother 0.6724 1.4873 0.3679 1.2287
## prio 1.0955 0.9128 1.0371 1.1571
## educ 0.8180 1.2225 0.6327 1.0576
##
## Concordance= 0.648 (se = 0.027 )
## Likelihood ratio test= 33.04 on 5 df, p=4e-06
## Wald test = 31.44 on 5 df, p=8e-06
## Score (logrank) test = 32.42 on 5 df, p=5e-06coef(modelCox_adjusted)## finyes age raceother prio educ
## -0.34522629 -0.06624708 -0.39697417 0.09118848 -0.20089618exp(coef(modelCox_adjusted))## finyes age raceother prio educ
## 0.7080601 0.9358996 0.6723514 1.0954755 0.8179974exp(confint(modelCox_adjusted))## 2.5 % 97.5 %
## finyes 0.4868945 1.0296874
## age 0.8986600 0.9746823
## raceother 0.3679100 1.2287145
## prio 1.0370926 1.1571449
## educ 0.6326966 1.0575680Mô hình đánh giá mối liên quan giữa hỗ trợ tài chính và thời gian bị bắt lại sau khi hiệu chỉnh cho tuổi, chủng tộc, số lần kết án trước và trình độ học vấn.
3 3. Xây dựng mô hình dự báo xác suất bị bắt lại
df_Rossi <- Rossi[
c(
"arrest",
"fin",
"age",
"race",
"wexp",
"mar",
"paro",
"prio",
"educ"
)
]
df_Rossi <- na.omit(df_Rossi)
summary(df_Rossi)## arrest fin age race wexp
## Min. :0.0000 no :216 Min. :17.0 black:379 no :185
## 1st Qu.:0.0000 yes:216 1st Qu.:20.0 other: 53 yes:247
## Median :0.0000 Median :23.0
## Mean :0.2639 Mean :24.6
## 3rd Qu.:1.0000 3rd Qu.:27.0
## Max. :1.0000 Max. :44.0
## mar paro prio educ
## married : 53 no :165 Min. : 0.000 Min. :2.000
## not married:379 yes:267 1st Qu.: 1.000 1st Qu.:3.000
## Median : 2.000 Median :3.000
## Mean : 2.984 Mean :3.477
## 3rd Qu.: 4.000 3rd Qu.:4.000
## Max. :18.000 Max. :6.0003.1 3.1 Bayesian Model Averaging – BMA
library(BMA)
model_BMA_Rossi <- bic.glm(
arrest ~ .,
data = df_Rossi,
glm.family = "binomial"
)
summary(model_BMA_Rossi)##
## Call:
## bic.glm.formula(f = arrest ~ ., data = df_Rossi, glm.family = "binomial")
##
##
## 12 models were selected
## Best 5 models (cumulative posterior probability = 0.7985 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 0.817350 0.81303 5.022e-01 9.561e-01 1.466e+00
## fin 14.6
## .yes -0.058569 0.16599 . . .
## age 100.0 -0.078012 0.02328 -7.796e-02 -8.305e-02 -7.692e-02
## race 3.1
## .other -0.010614 0.08693 . . .
## wexp 4.9
## .yes -0.014130 0.08348 . . .
## mar 5.0
## .not married 0.027326 0.15084 . . .
## paro 2.1
## .yes -0.002227 0.03709 . . .
## prio 75.0 0.077046 0.05457 1.046e-01 . 9.433e-02
## educ 20.8 -0.062332 0.14083 . . -2.802e-01
##
## nVar 2 1 3
## BIC -2.129e+03 -2.127e+03 -2.126e+03
## post prob 0.395 0.122 0.107
## model 4 model 5
## Intercept 6.565e-01 2.054e+00
## fin
## .yes -4.057e-01 .
## age -7.648e-02 -8.108e-02
## race
## .other . .
## wexp
## .yes . .
## mar
## .not married . .
## paro
## .yes . .
## prio 1.056e-01 .
## educ . -3.351e-01
##
## nVar 3 2
## BIC -2.126e+03 -2.126e+03
## post prob 0.095 0.079
##
## 1 observations deleted due to missingness.3.1.1 Biểu đồ BMA
par(
mfrow = c(1, 1),
mar = c(5, 5, 4, 2)
)
imageplot.bma(model_BMA_Rossi)
BMA đánh giá nhiều mô hình ứng viên thay vì phụ thuộc hoàn toàn vào một mô hình duy nhất. Các biến có xác suất hậu nghiệm cao hơn có bằng chứng mạnh hơn về vai trò dự báo.
3.2 3.2 Dự báo xác suất bị bắt lại đến tuần thứ 30
Xét hai người:
- Người 1: 30 tuổi và chưa từng bị kết án trước đó
(
prio = 0). - Người 2: 35 tuổi và có một lần bị kết án trước đó
(
prio = 1).
modelCox_predict <- coxph(
Surv(week, arrest) ~ age + prio,
data = Rossi
)
summary(modelCox_predict)## Call:
## coxph(formula = Surv(week, arrest) ~ age + prio, data = Rossi)
##
## n= 432, number of events= 114
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.06934 0.93301 0.02079 -3.335 0.000854 ***
## prio 0.09466 1.09928 0.02711 3.492 0.000479 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.933 1.0718 0.8958 0.9718
## prio 1.099 0.9097 1.0424 1.1593
##
## Concordance= 0.633 (se = 0.028 )
## Likelihood ratio test= 25.68 on 2 df, p=3e-06
## Wald test = 24.51 on 2 df, p=5e-06
## Score (logrank) test = 25.63 on 2 df, p=3e-06people_con <- data.frame(
age = c(30, 35),
prio = c(0, 1)
)
people_con## age prio
## 1 30 0
## 2 35 13.2.1 Dự báo tại tuần 30
pred_cox <- survfit(
modelCox_predict,
newdata = people_con
)
pred_30 <- summary(
pred_cox,
times = 30
)
surv_30 <- as.vector(
pred_30$surv
)
pre_arrested <- 1 - surv_30
result_30 <- data.frame(
Nguoi = c("Người 1", "Người 2"),
age = c(30, 35),
prior = c(0, 1),
Xac_suat_chua_bi_bat_lai = surv_30,
Xac_suat_bi_bat_lai = pre_arrested
)
result_30## Nguoi age prior Xac_suat_chua_bi_bat_lai Xac_suat_bi_bat_lai
## 1 Người 1 30 0 0.9329806 0.06701942
## 2 Người 2 35 1 0.9475105 0.05248945Xác suất bị bắt lại đến tuần 30 được tính:
\[ P(T \leq 30) = 1 - S(30) \]
Theo kết quả dự báo:
## Người 1 (30 tuổi, prior = 0): 6.7%## Người 2 (35 tuổi, prior = 1): 5.25%3.3 3.3 Bootstrap validation và Calibration
library(rms)
dd <- datadist(Rossi)
options(datadist = "dd")
model_boot <- lrm(
arrest ~ fin + age + race + wexp +
mar + paro + prio + educ,
data = Rossi,
x = TRUE,
y = TRUE
)
validate_boot <- validate(
model_boot,
method = "boot",
B = 100
)
validate_boot## index.orig training test optimism index.corrected Lower Upper
## Dxy 0.3671 0.4043 0.3316 0.0727 0.2944 0.1765 0.3952
## R2 0.1113 0.1413 0.0922 0.0491 0.0623 -0.0186 0.1281
## Intercept 0.0000 0.0000 -0.1956 0.1956 -0.1956 -0.5346 0.2311
## Slope 1.0000 1.0000 0.7876 0.2124 0.7876 0.4847 1.1701
## Emax 0.0000 0.0000 0.1019 -0.1019 0.1019 -0.0248 0.2545
## D 0.0770 0.1000 0.0629 0.0370 0.0400 -0.0225 0.0911
## U -0.0046 -0.0046 0.0047 -0.0093 0.0047 -0.0093 0.0306
## Q 0.0816 0.1046 0.0583 0.0463 0.0353 -0.0529 0.0965
## B 0.1786 0.1742 0.1825 -0.0083 0.1869 0.1672 0.2074
## g 0.7922 0.9216 0.7090 0.2125 0.5796 0.2427 0.8731
## gp 0.1376 0.1536 0.1244 0.0292 0.1083 0.0631 0.1529
## n
## Dxy 100
## R2 100
## Intercept 100
## Slope 100
## Emax 100
## D 100
## U 100
## Q 100
## B 100
## g 100
## gp 100Bootstrap được sử dụng để đánh giá nội bộ mức độ ổn định và khả năng dự báo của mô hình thông qua việc lấy mẫu lặp lại từ dữ liệu ban đầu.
3.3.1 Calibration plot
cal_boot <- calibrate(
model_boot,
method = "boot",
B = 100
)
plot(
cal_boot,
xlab = "Xác suất dự báo",
ylab = "Xác suất quan sát",
main = "Calibration Plot - Bootstrap"
)##
## n=432 Mean absolute error=0.011 Mean squared error=0.00028
## 0.9 Quantile of absolute error=0.032abline(
0,
1,
lty = 2,
lwd = 2
)
Đường calibration càng gần đường chéo 45 độ thì xác suất mô hình dự báo càng gần xác suất thực tế quan sát được.
4 4. Đánh giá tầm quan trọng của biến số
db_Rossi <- Rossi[
c(
"arrest",
"fin",
"age",
"race",
"wexp",
"mar",
"paro",
"prio",
"educ"
)
]
df_model <- na.omit(db_Rossi)
summary(df_model)## arrest fin age race wexp
## Min. :0.0000 no :216 Min. :17.0 black:379 no :185
## 1st Qu.:0.0000 yes:216 1st Qu.:20.0 other: 53 yes:247
## Median :0.0000 Median :23.0
## Mean :0.2639 Mean :24.6
## 3rd Qu.:1.0000 3rd Qu.:27.0
## Max. :1.0000 Max. :44.0
## mar paro prio educ
## married : 53 no :165 Min. : 0.000 Min. :2.000
## not married:379 yes:267 1st Qu.: 1.000 1st Qu.:3.000
## Median : 2.000 Median :3.000
## Mean : 2.984 Mean :3.477
## 3rd Qu.: 4.000 3rd Qu.:4.000
## Max. :18.000 Max. :6.0004.1 4.1 Đánh giá bằng BMA
modelBMA_testing <- bic.glm(
arrest ~
fin + age + race + wexp +
mar + paro + prio + educ,
data = df_model,
glm.family = "binomial"
)
summary(modelBMA_testing)##
## Call:
## bic.glm.formula(f = arrest ~ fin + age + race + wexp + mar + paro + prio + educ, data = df_model, glm.family = "binomial")
##
##
## 12 models were selected
## Best 5 models (cumulative posterior probability = 0.7985 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100 0.817350 0.81303 5.022e-01 9.561e-01 1.466e+00
## fin 14.6
## .yes -0.058569 0.16599 . . .
## age 100.0 -0.078012 0.02328 -7.796e-02 -8.305e-02 -7.692e-02
## race 3.1
## .other -0.010614 0.08693 . . .
## wexp 4.9
## .yes -0.014130 0.08348 . . .
## mar 5.0
## .not married 0.027326 0.15084 . . .
## paro 2.1
## .yes -0.002227 0.03709 . . .
## prio 75.0 0.077046 0.05457 1.046e-01 . 9.433e-02
## educ 20.8 -0.062332 0.14083 . . -2.802e-01
##
## nVar 2 1 3
## BIC -2.129e+03 -2.127e+03 -2.126e+03
## post prob 0.395 0.122 0.107
## model 4 model 5
## Intercept 6.565e-01 2.054e+00
## fin
## .yes -4.057e-01 .
## age -7.648e-02 -8.108e-02
## race
## .other . .
## wexp
## .yes . .
## mar
## .not married . .
## paro
## .yes . .
## prio 1.056e-01 .
## educ . -3.351e-01
##
## nVar 3 2
## BIC -2.126e+03 -2.126e+03
## post prob 0.095 0.079
##
## 1 observations deleted due to missingness.BMA đánh giá tầm quan trọng dựa trên mức độ hỗ trợ hậu nghiệm của các biến trong tập hợp các mô hình ứng viên.
4.2 4.2 Đánh giá bằng LASSO
library(glmnet)
X <- model.matrix(
arrest ~
fin + age + race + wexp +
mar + paro + prio + educ,
data = df_model
)[, -1]
Y <- df_model$arrest
modelLASSO_testing <- cv.glmnet(
x = X,
y = Y,
family = "binomial",
alpha = 1
)
modelLASSO_testing$lambda.min## [1] 0.009327493coef_LASSO <- coef(
modelLASSO_testing,
s = "lambda.min"
)
coef_LASSO## 9 x 1 sparse Matrix of class "dgCMatrix"
## lambda.min
## (Intercept) 0.69196727
## finyes -0.30560446
## age -0.05763655
## raceother -0.20255291
## wexpyes -0.08735927
## marnot married 0.33353467
## paroyes .
## prio 0.08227426
## educ -0.196956144.2.1 Bảng tầm quan trọng của biến theo LASSO
coef_table <- data.frame(
Variable = rownames(as.matrix(coef_LASSO)),
Coefficient = as.vector(coef_LASSO)
)
coef_table <- coef_table[
coef_table$Variable != "(Intercept)",
]
coef_table$Importance <-
abs(coef_table$Coefficient)
coef_table <- coef_table[
order(
coef_table$Importance,
decreasing = TRUE
),
]
coef_table## Variable Coefficient Importance
## 6 marnot married 0.33353467 0.33353467
## 2 finyes -0.30560446 0.30560446
## 4 raceother -0.20255291 0.20255291
## 9 educ -0.19695614 0.19695614
## 5 wexpyes -0.08735927 0.08735927
## 8 prio 0.08227426 0.08227426
## 3 age -0.05763655 0.05763655
## 7 paroyes 0.00000000 0.000000004.2.2 Biểu đồ BMA
imageplot.bma(
modelBMA_testing
)
4.2.3 Biểu đồ Variable Importance – LASSO
barplot(
rev(coef_table$Importance),
names.arg = rev(coef_table$Variable),
horiz = TRUE,
las = 1,
main = "Variable Importance - LASSO",
xlab = "|LASSO coefficient|"
)
Các biến có hệ số LASSO khác 0 được giữ lại trong mô hình. Biến có trị tuyệt đối hệ số lớn hơn có ảnh hưởng mạnh hơn trong mô hình LASSO, trong khi các biến có hệ số bằng 0 bị loại bỏ.
Việc kết hợp kết quả BMA và LASSO giúp đánh giá mức độ nhất quán trong việc xác định các yếu tố dự báo quan trọng.
5 Bonus – Biểu đồ 3D dự báo xác suất bị bắt lại
Biểu đồ dưới đây thể hiện:
- Trục X: tuần theo dõi.
- Trục Y: hai đối tượng cần dự báo.
- Trục Z: xác suất bị bắt lại tích lũy.
- Điểm cuối biểu diễn kết quả tại tuần thứ 30.
library(plotly)
pred_all <- summary(
pred_cox,
times = 1:30
)
risk_all <- 1 - pred_all$surv
week_values <- 1:30
risk_person1 <-
risk_all[, 1] * 100
risk_person2 <-
risk_all[, 2] * 100plot_3D <- plot_ly()
plot_3D <- plot_3D %>%
add_trace(
x = week_values,
y = rep(1, 30),
z = risk_person1,
type = "scatter3d",
mode = "lines+markers",
name = "Người 1: age 30, prio 0",
text = paste0(
"Người 1",
"<br>Age = 30",
"<br>Prio = 0",
"<br>Tuần = ", week_values,
"<br>Xác suất = ",
round(risk_person1, 2),
"%"
),
hoverinfo = "text",
line = list(width = 5),
marker = list(size = 3)
)
plot_3D <- plot_3D %>%
add_trace(
x = week_values,
y = rep(2, 30),
z = risk_person2,
type = "scatter3d",
mode = "lines+markers",
name = "Người 2: age 35, prio 1",
text = paste0(
"Người 2",
"<br>Age = 35",
"<br>Prio = 1",
"<br>Tuần = ", week_values,
"<br>Xác suất = ",
round(risk_person2, 2),
"%"
),
hoverinfo = "text",
line = list(width = 5),
marker = list(size = 3)
)5.0.1 Đánh dấu kết quả tại tuần 30
plot_3D <- plot_3D %>%
add_markers(
x = c(30, 30),
y = c(1, 2),
z = pre_arrested * 100,
type = "scatter3d",
mode = "markers",
text = c(
paste0(
"Người 1 - Tuần 30",
"<br>Age = 30",
"<br>Prio = 0",
"<br>Xác suất bị bắt lại = ",
round(pre_arrested[1] * 100, 2),
"%"
),
paste0(
"Người 2 - Tuần 30",
"<br>Age = 35",
"<br>Prio = 1",
"<br>Xác suất bị bắt lại = ",
round(pre_arrested[2] * 100, 2),
"%"
)
),
hoverinfo = "text",
marker = list(size = 8),
name = "Kết quả tuần 30"
)5.0.2 Hiển thị biểu đồ 3D
plot_3D <- plot_3D %>%
layout(
title = list(
text = "Dự báo xác suất bị bắt lại đến tuần 30",
x = 0.5
),
scene = list(
xaxis = list(
title = "Tuần",
range = c(1, 30),
dtick = 5
),
yaxis = list(
title = "Đối tượng",
tickvals = c(1, 2),
ticktext = c(
"Người 1: 30 tuổi, prio 0",
"Người 2: 35 tuổi, prio 1"
)
),
zaxis = list(
title = "Xác suất bị bắt lại (%)",
ticksuffix = "%"
),
camera = list(
eye = list(
x = 1.5,
y = 1.5,
z = 1.1
)
)
)
)
plot_3D6 Kết luận
Phân tích dữ liệu sống còn cho phép đánh giá không chỉ việc một biến cố có xảy ra hay không mà còn xét đến thời gian xảy ra biến cố.
Kaplan-Meier và hồi quy Cox được sử dụng để mô tả và đánh giá nguy cơ theo thời gian. BMA và LASSO hỗ trợ lựa chọn và đánh giá vai trò của các biến dự báo. Bootstrap và calibration được sử dụng để đánh giá nội bộ khả năng dự báo của mô hình.
Ngoài ra, trực quan hóa 3D cho phép quan sát sự thay đổi của xác suất bị bắt lại theo thời gian đối với hai đối tượng cụ thể trong bài toán dự báo.