Data Loading
data("bmt")
bmt <- bmt
head(bmt)
## group t1 t2 d1 d2 d3 ta da tc dc tp dp z1 z2 z3 z4 z5 z6 z7 z8 z9
## 1 1 2081 2081 0 0 0 67 1 121 1 13 1 26 33 1 0 1 1 98 0 1
## 2 1 1602 1602 0 0 0 1602 0 139 1 18 1 21 37 1 1 0 0 1720 0 1
## 3 1 1496 1496 0 0 0 1496 0 307 1 12 1 26 35 1 1 1 0 127 0 1
## 4 1 1462 1462 0 0 0 70 1 95 1 13 1 17 21 0 1 0 0 168 0 1
## 5 1 1433 1433 0 0 0 1433 0 236 1 12 1 32 36 1 1 1 1 93 0 1
## 6 1 1377 1377 0 0 0 1377 0 123 1 12 1 22 31 1 1 1 1 2187 0 1
## z10
## 1 0
## 2 0
## 3 0
## 4 0
## 5 0
## 6 0
# variabel data BMT
str(bmt)
## 'data.frame': 137 obs. of 22 variables:
## $ group: int 1 1 1 1 1 1 1 1 1 1 ...
## $ t1 : int 2081 1602 1496 1462 1433 1377 1330 996 226 1199 ...
## $ t2 : int 2081 1602 1496 1462 1433 1377 1330 996 226 1199 ...
## $ d1 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ d2 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ d3 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ ta : int 67 1602 1496 70 1433 1377 1330 72 226 1199 ...
## $ da : int 1 0 0 1 0 0 0 1 0 0 ...
## $ tc : int 121 139 307 95 236 123 96 121 226 91 ...
## $ dc : int 1 1 1 1 1 1 1 1 0 1 ...
## $ tp : int 13 18 12 13 12 12 17 12 10 29 ...
## $ dp : int 1 1 1 1 1 1 1 1 1 1 ...
## $ z1 : int 26 21 26 17 32 22 20 22 18 24 ...
## $ z2 : int 33 37 35 21 36 31 17 24 21 40 ...
## $ z3 : int 1 1 1 0 1 1 1 1 0 1 ...
## $ z4 : int 0 1 1 1 1 1 0 0 1 1 ...
## $ z5 : int 1 0 1 0 1 1 1 0 0 0 ...
## $ z6 : int 1 0 0 0 1 1 1 0 0 1 ...
## $ z7 : int 98 1720 127 168 93 2187 1006 1319 208 174 ...
## $ z8 : int 0 0 0 0 0 0 0 0 0 0 ...
## $ z9 : int 1 1 1 1 1 1 1 1 1 3 ...
## $ z10 : int 0 0 0 0 0 0 0 0 0 1 ...
# mendefinisikan survival time dan status censored
time <- bmt$t2
status_d <- bmt$d3
y <- Surv(time, status_d)
# mendefinisikan variabel grup
bmt$group <- factor(
bmt$group,
levels = c(1, 2, 3),
labels = c("ALL", "AML Low Risk", "AML High Risk"))
table(bmt$group)
##
## ALL AML Low Risk AML High Risk
## 38 54 45
K-M Curve Overall
kmfit0 <- survfit(y~1)
kmfit0
## Call: survfit(formula = y ~ 1)
##
## n events median 0.95LCL 0.95UCL
## [1,] 137 83 481 381 1063
kmfit0 %>%
ggsurvfit() +
labs(
x = "Days",
y = "Overall Survival Probability",
title = "Kurva Kaplan-Meier Overall (Seluruh Pasien BMT)"
) +
add_confidence_interval()
K-M Curve Berdasarkan Kelompok Resiko
kmfit1 <- survfit(y~bmt$group)
kmfit1
## Call: survfit(formula = y ~ bmt$group)
##
## n events median 0.95LCL 0.95UCL
## bmt$group=ALL 38 24 418 194 NA
## bmt$group=AML Low Risk 54 25 2204 704 NA
## bmt$group=AML High Risk 45 34 183 115 456
kmfit1 %>%
ggsurvfit() +
labs(
x = "Days",
y = "Survival Probability",
title = "Kurva Kaplan-Meier Berdasarkan Kelompok Risiko"
) +
add_confidence_interval()
Uji Log Rank
lr1 <- survdiff(y ~ bmt$group)
lr1
## Call:
## survdiff(formula = y ~ bmt$group)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## bmt$group=ALL 38 24 21.9 0.211 0.289
## bmt$group=AML Low Risk 54 25 40.0 5.604 11.012
## bmt$group=AML High Risk 45 34 21.2 7.756 10.529
##
## Chisq= 13.8 on 2 degrees of freedom, p= 0.001
Pemodelan
Screening Univariabel
covariates <- c("z1", "z2", "z3", "z4", "z5", "z6", "z7", "z8", "z9", "z10")
univ_results <- data.frame()
for (var in covariates) {
formula <- as.formula(paste0("y ~", var))
cox_univ <- coxph(formula, data = bmt)
sum_univ <- summary(cox_univ)
hr <- round(sum_univ$coefficients[,2], 3)
ci_lower <- round(sum_univ$conf.int[,3], 3)
ci_upper <- round(sum_univ$conf.int[,4], 3)
pval <- round(sum_univ$coefficients[,5], 4)
univ_results <- rbind(univ_results,
data.frame(Variable = var,
HR_CI = paste0(hr, " (", ci_lower, " - ", ci_upper, ")"),
P_value = pval))}
print(univ_results)
## Variable HR_CI P_value
## 1 z1 1.011 (0.988 - 1.035) 0.3377
## 2 z2 1.014 (0.99 - 1.04) 0.2515
## 3 z3 0.795 (0.514 - 1.228) 0.3009
## 4 z4 0.991 (0.633 - 1.552) 0.9701
## 5 z5 1.167 (0.759 - 1.796) 0.4818
## 6 z6 1.047 (0.677 - 1.62) 0.8360
## 7 z7 1 (0.999 - 1.001) 0.7910
## 8 z8 1.892 (1.222 - 2.929) 0.0043
## 9 z9 0.835 (0.688 - 1.013) 0.0675
## 10 z10 1.489 (0.934 - 2.373) 0.0944
Model Penuh
cox_full <- coxph(y ~ z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z9 + z10, data = bmt)
summary(cox_full)
## Call:
## coxph(formula = y ~ z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z9 +
## z10, data = bmt)
##
## n= 137, number of events= 83
##
## coef exp(coef) se(coef) z Pr(>|z|)
## z1 6.990e-03 1.007e+00 2.025e-02 0.345 0.729997
## z2 1.526e-02 1.015e+00 1.912e-02 0.798 0.424735
## z3 -2.173e-01 8.047e-01 2.429e-01 -0.894 0.371150
## z4 9.792e-02 1.103e+00 2.416e-01 0.405 0.685249
## z5 -1.654e-01 8.475e-01 2.480e-01 -0.667 0.504841
## z6 1.287e-02 1.013e+00 2.399e-01 0.054 0.957215
## z7 -5.869e-05 9.999e-01 3.562e-04 -0.165 0.869125
## z8 7.883e-01 2.200e+00 2.384e-01 3.306 0.000946 ***
## z9 -4.570e-01 6.332e-01 1.444e-01 -3.165 0.001553 **
## z10 9.042e-01 2.470e+00 3.130e-01 2.889 0.003871 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## z1 1.0070 0.9930 0.9678 1.0478
## z2 1.0154 0.9849 0.9780 1.0542
## z3 0.8047 1.2427 0.4999 1.2955
## z4 1.1029 0.9067 0.6869 1.7708
## z5 0.8475 1.1799 0.5212 1.3781
## z6 1.0130 0.9872 0.6330 1.6210
## z7 0.9999 1.0001 0.9992 1.0006
## z8 2.1995 0.4546 1.3784 3.5098
## z9 0.6332 1.5793 0.4771 0.8403
## z10 2.4700 0.4049 1.3373 4.5621
##
## Concordance= 0.664 (se = 0.031 )
## Likelihood ratio test= 25.45 on 10 df, p=0.005
## Wald test = 23.05 on 10 df, p=0.01
## Score (logrank) test = 23.67 on 10 df, p=0.009
Seleksi Variabel
# Stepwise backward
cox_selection <- step(cox_full, direction = "backward", trace = 1)
## Start: AIC=741.15
## y ~ z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z9 + z10
##
## Df AIC
## - z6 1 739.15
## - z7 1 739.17
## - z1 1 739.26
## - z4 1 739.31
## - z5 1 739.59
## - z2 1 739.79
## - z3 1 739.94
## <none> 741.15
## - z10 1 747.77
## - z8 1 749.64
## - z9 1 752.05
##
## Step: AIC=739.15
## y ~ z1 + z2 + z3 + z4 + z5 + z7 + z8 + z9 + z10
##
## Df AIC
## - z7 1 737.18
## - z1 1 737.26
## - z4 1 737.31
## - z5 1 737.60
## - z2 1 737.84
## - z3 1 737.97
## <none> 739.15
## - z10 1 745.82
## - z8 1 747.88
## - z9 1 750.07
##
## Step: AIC=737.18
## y ~ z1 + z2 + z3 + z4 + z5 + z8 + z9 + z10
##
## Df AIC
## - z1 1 735.32
## - z4 1 735.37
## - z5 1 735.64
## - z2 1 735.85
## - z3 1 736.10
## <none> 737.18
## - z10 1 743.84
## - z8 1 746.06
## - z9 1 748.08
##
## Step: AIC=735.32
## y ~ z2 + z3 + z4 + z5 + z8 + z9 + z10
##
## Df AIC
## - z4 1 733.45
## - z5 1 733.69
## - z3 1 734.27
## <none> 735.32
## - z2 1 735.94
## - z10 1 742.11
## - z8 1 744.42
## - z9 1 746.18
##
## Step: AIC=733.45
## y ~ z2 + z3 + z5 + z8 + z9 + z10
##
## Df AIC
## - z5 1 731.80
## - z3 1 732.33
## <none> 733.45
## - z2 1 734.05
## - z10 1 740.16
## - z8 1 742.46
## - z9 1 744.18
##
## Step: AIC=731.8
## y ~ z2 + z3 + z8 + z9 + z10
##
## Df AIC
## - z3 1 730.87
## <none> 731.80
## - z2 1 732.21
## - z10 1 738.16
## - z8 1 740.49
## - z9 1 742.18
##
## Step: AIC=730.87
## y ~ z2 + z8 + z9 + z10
##
## Df AIC
## <none> 730.87
## - z2 1 731.31
## - z10 1 737.79
## - z8 1 739.74
## - z9 1 740.34
summary(cox_selection)
## Call:
## coxph(formula = y ~ z2 + z8 + z9 + z10, data = bmt)
##
## n= 137, number of events= 83
##
## coef exp(coef) se(coef) z Pr(>|z|)
## z2 0.01946 1.01965 0.01247 1.561 0.118620
## z8 0.76300 2.14469 0.22473 3.395 0.000686 ***
## z9 -0.41293 0.66171 0.13965 -2.957 0.003107 **
## z10 0.88806 2.43041 0.30081 2.952 0.003155 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## z2 1.0197 0.9807 0.9950 1.045
## z8 2.1447 0.4663 1.3806 3.332
## z9 0.6617 1.5112 0.5033 0.870
## z10 2.4304 0.4115 1.3478 4.383
##
## Concordance= 0.665 (se = 0.03 )
## Likelihood ratio test= 23.72 on 4 df, p=9e-05
## Wald test = 21.69 on 4 df, p=2e-04
## Score (logrank) test = 22.33 on 4 df, p=2e-04
Model Hasil Seleksi
cox_pars <- coxph(y ~ z2 + z8 + z9 + z10, data = bmt)
summary(cox_pars)
## Call:
## coxph(formula = y ~ z2 + z8 + z9 + z10, data = bmt)
##
## n= 137, number of events= 83
##
## coef exp(coef) se(coef) z Pr(>|z|)
## z2 0.01946 1.01965 0.01247 1.561 0.118620
## z8 0.76300 2.14469 0.22473 3.395 0.000686 ***
## z9 -0.41293 0.66171 0.13965 -2.957 0.003107 **
## z10 0.88806 2.43041 0.30081 2.952 0.003155 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## z2 1.0197 0.9807 0.9950 1.045
## z8 2.1447 0.4663 1.3806 3.332
## z9 0.6617 1.5112 0.5033 0.870
## z10 2.4304 0.4115 1.3478 4.383
##
## Concordance= 0.665 (se = 0.03 )
## Likelihood ratio test= 23.72 on 4 df, p=9e-05
## Wald test = 21.69 on 4 df, p=2e-04
## Score (logrank) test = 22.33 on 4 df, p=2e-04
anova(cox_full, cox_pars, test = "LRT")
## Analysis of Deviance Table
## Cox model: response is y
## Model 1: ~ z1 + z2 + z3 + z4 + z5 + z6 + z7 + z8 + z9 + z10
## Model 2: ~ z2 + z8 + z9 + z10
## loglik Chisq Df Pr(>|Chi|)
## 1 -360.57
## 2 -361.44 1.7273 6 0.943
Model dengan Interaksi
cox_interaksi <- coxph(y ~ z2 + z8 + z9 + z10 + z8:z10, data=bmt)
summary(cox_interaksi)
## Call:
## coxph(formula = y ~ z2 + z8 + z9 + z10 + z8:z10, data = bmt)
##
## n= 137, number of events= 83
##
## coef exp(coef) se(coef) z Pr(>|z|)
## z2 0.02329 1.02357 0.01233 1.889 0.05890 .
## z8 1.16617 3.20968 0.27115 4.301 1.70e-05 ***
## z9 -0.44783 0.63901 0.14174 -3.159 0.00158 **
## z10 1.42202 4.14550 0.35850 3.967 7.29e-05 ***
## z8:z10 -1.44256 0.23632 0.56952 -2.533 0.01131 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## z2 1.0236 0.9770 0.9991 1.0486
## z8 3.2097 0.3116 1.8865 5.4609
## z9 0.6390 1.5649 0.4840 0.8437
## z10 4.1455 0.2412 2.0531 8.3702
## z8:z10 0.2363 4.2315 0.0774 0.7216
##
## Concordance= 0.685 (se = 0.028 )
## Likelihood ratio test= 31.1 on 5 df, p=9e-06
## Wald test = 27.05 on 5 df, p=6e-05
## Score (logrank) test = 28.11 on 5 df, p=3e-05
anova(cox_pars, cox_interaksi, test = "LRT")
## Analysis of Deviance Table
## Cox model: response is y
## Model 1: ~ z2 + z8 + z9 + z10
## Model 2: ~ z2 + z8 + z9 + z10 + z8:z10
## loglik Chisq Df Pr(>|Chi|)
## 1 -361.44
## 2 -357.74 7.3848 1 0.006578 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Stratified Cox
cox_strat <- coxph(y ~ z2 + z8 + z9 + z8:z10 + strata(z10), data = bmt)
summary(cox_strat)
## Call:
## coxph(formula = y ~ z2 + z8 + z9 + z8:z10 + strata(z10), data = bmt)
##
## n= 137, number of events= 83
##
## coef exp(coef) se(coef) z Pr(>|z|)
## z2 0.02024 1.02044 0.01211 1.671 0.09474 .
## z8 1.22625 3.40842 0.27281 4.495 6.96e-06 ***
## z9 -0.45672 0.63336 0.14009 -3.260 0.00111 **
## z8:z10 -1.48226 0.22712 0.57100 -2.596 0.00943 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## z2 1.0204 0.9800 0.99650 1.0450
## z8 3.4084 0.2934 1.99681 5.8179
## z9 0.6334 1.5789 0.48128 0.8335
## z8:z10 0.2271 4.4029 0.07417 0.6955
##
## Concordance= 0.674 (se = 0.031 )
## Likelihood ratio test= 30 on 4 df, p=5e-06
## Wald test = 27.95 on 4 df, p=1e-05
## Score (logrank) test = 29.5 on 4 df, p=6e-06
Uji Asumsi PH Model Stratified
ph_test_strat <- cox.zph(cox_strat)
print(ph_test_strat)
## chisq df p
## z2 1.3396 1 0.25
## z8 0.1115 1 0.74
## z9 0.0659 1 0.80
## z8:z10 0.6503 1 0.42
## GLOBAL 2.8187 4 0.59
