# Load package
library(survival)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# Baca data
library(readxl)
data <- read_excel("C:/Users/62882/Documents/3SK3/SEMESTER 6/ATH/PROJEK/data_yokbisa.xlsx",sheet = "timur_fix")
## New names:
## • `` -> `...41`
## • `` -> `...42`
head(data)
## # A tibble: 6 × 94
## FWT R101 R102 R105 R301 R302 R303 R304 R305 R401 R403 R404
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 120. 74 4 2 6 0 6 6 2 4 3 1
## 2 281. 73 13 1 16 3 13 12 3 7 6 1
## 3 47.2 73 11 2 5 0 5 4 1 3 3 1
## 4 47.2 73 11 2 5 0 5 4 1 4 3 1
## 5 122. 73 22 2 3 0 3 3 0 2 6 1
## 6 83.6 73 18 2 4 0 4 3 1 3 3 1
## # ℹ 82 more variables: `R405(JK)` <dbl>, R406A <dbl>, R406B <dbl>, R406C <dbl>,
## # R407 <dbl>, R408 <dbl>, R409 <dbl>, R609 <dbl>, R612 <dbl>, R613 <dbl>,
## # R614 <dbl>, R615 <dbl>, Time <dbl>, Sensor <dbl>, R616 <dbl>, R617 <dbl>,
## # R618 <dbl>, R619 <dbl>, R620 <dbl>, R702_A <chr>, R702_B <chr>,
## # R702_C <chr>, R702_D <chr>, R702_X <chr>, `R703(Kat_kegiatan)` <dbl>,
## # kat_kegiatan <dbl>, R704 <dbl>, R705 <dbl>, ...41 <chr>, ...42 <lgl>,
## # R706 <dbl>, R707 <dbl>, R708 <dbl>, R801 <dbl>, R802 <dbl>, R807_A <chr>, …
# Objek survival
surv_obj <- Surv(data$Time, data$Sensor)
# Model awal (semua variabel X)
cox_model <- coxph(surv_obj ~ Klasifikasi + Kat_KIP + Internet +Kat_Kegiatan+JK, data = data)
# Ringkasan model
summary(cox_model)
## Call:
## coxph(formula = surv_obj ~ Klasifikasi + Kat_KIP + Internet +
## Kat_Kegiatan + JK, data = data)
##
## n= 1728, number of events= 80
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Klasifikasi -0.4311 0.6498 0.3121 -1.382 0.16711
## Kat_KIP -1.0529 0.3489 0.3977 -2.648 0.00811 **
## Internet -0.7144 0.4895 0.2405 -2.970 0.00298 **
## Kat_Kegiatan 1.9970 7.3670 0.2472 8.079 6.56e-16 ***
## JK -0.2023 0.8168 0.2335 -0.867 0.38614
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Klasifikasi 0.6498 1.5390 0.3525 1.1978
## Kat_KIP 0.3489 2.8660 0.1600 0.7607
## Internet 0.4895 2.0429 0.3055 0.7843
## Kat_Kegiatan 7.3670 0.1357 4.5381 11.9593
## JK 0.8168 1.2242 0.5169 1.2908
##
## Concordance= 0.753 (se = 0.033 )
## Likelihood ratio test= 80.16 on 5 df, p=8e-16
## Wald test = 100.9 on 5 df, p=<2e-16
## Score (logrank) test = 138.5 on 5 df, p=<2e-16
# Uji asumsi PH
ph_test <- cox.zph(cox_model)
print(ph_test)
## chisq df p
## Klasifikasi 2.2506 1 0.134
## Kat_KIP 0.0761 1 0.783
## Internet 0.8676 1 0.352
## Kat_Kegiatan 3.3648 1 0.067
## JK 2.2880 1 0.130
## GLOBAL 6.6430 5 0.249
plot(ph_test) # Visualisasi asumsi
## Warning in plot.cox.zph(ph_test): spline fit is singular, variable 1 skipped
## Warning in plot.cox.zph(ph_test): spline fit is singular, variable 2 skipped
## Warning in plot.cox.zph(ph_test): spline fit is singular, variable 3 skipped
## Warning in plot.cox.zph(ph_test): spline fit is singular, variable 4 skipped
## Warning in plot.cox.zph(ph_test): spline fit is singular, variable 5 skipped
# Uji simultan dan AIC
anova(cox_model)
## Analysis of Deviance Table
## Cox model: response is surv_obj
## Terms added sequentially (first to last)
##
## loglik Chisq Df Pr(>|Chi|)
## NULL -540.36
## Klasifikasi -536.41 7.9019 1 0.0049382 **
## Kat_KIP -531.78 9.2426 1 0.0023645 **
## Internet -526.27 11.0344 1 0.0008944 ***
## Kat_Kegiatan -500.66 51.2172 1 8.269e-13 ***
## JK -500.28 0.7604 1 0.3832169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(cox_model)
## [1] 1010.556
library(MASS)
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
backward_model <- stepAIC(cox_model, direction = "backward")
## Start: AIC=1010.56
## surv_obj ~ Klasifikasi + Kat_KIP + Internet + Kat_Kegiatan +
## JK
##
## Df AIC
## - JK 1 1009.3
## <none> 1010.6
## - Klasifikasi 1 1010.6
## - Internet 1 1016.9
## - Kat_KIP 1 1017.9
## - Kat_Kegiatan 1 1057.9
##
## Step: AIC=1009.32
## surv_obj ~ Klasifikasi + Kat_KIP + Internet + Kat_Kegiatan
##
## Df AIC
## <none> 1009.3
## - Klasifikasi 1 1009.4
## - Internet 1 1016.0
## - Kat_KIP 1 1017.4
## - Kat_Kegiatan 1 1058.5
summary(backward_model)
## Call:
## coxph(formula = surv_obj ~ Klasifikasi + Kat_KIP + Internet +
## Kat_Kegiatan, data = data)
##
## n= 1728, number of events= 80
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Klasifikasi -0.4344 0.6477 0.3117 -1.394 0.16340
## Kat_KIP -1.0846 0.3380 0.3962 -2.738 0.00619 **
## Internet -0.7263 0.4837 0.2400 -3.026 0.00248 **
## Kat_Kegiatan 2.0283 7.6008 0.2446 8.292 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Klasifikasi 0.6477 1.5440 0.3516 1.1930
## Kat_KIP 0.3380 2.9583 0.1555 0.7348
## Internet 0.4837 2.0674 0.3022 0.7742
## Kat_Kegiatan 7.6008 0.1316 4.7061 12.2762
##
## Concordance= 0.745 (se = 0.033 )
## Likelihood ratio test= 79.4 on 4 df, p=2e-16
## Wald test = 99.53 on 4 df, p=<2e-16
## Score (logrank) test = 137.5 on 4 df, p=<2e-16
stepwise_model <- stepAIC(cox_model, direction = "both")
## Start: AIC=1010.56
## surv_obj ~ Klasifikasi + Kat_KIP + Internet + Kat_Kegiatan +
## JK
##
## Df AIC
## - JK 1 1009.3
## <none> 1010.6
## - Klasifikasi 1 1010.6
## - Internet 1 1016.9
## - Kat_KIP 1 1017.9
## - Kat_Kegiatan 1 1057.9
##
## Step: AIC=1009.32
## surv_obj ~ Klasifikasi + Kat_KIP + Internet + Kat_Kegiatan
##
## Df AIC
## <none> 1009.3
## - Klasifikasi 1 1009.4
## + JK 1 1010.6
## - Internet 1 1016.0
## - Kat_KIP 1 1017.4
## - Kat_Kegiatan 1 1058.5
summary(stepwise_model)
## Call:
## coxph(formula = surv_obj ~ Klasifikasi + Kat_KIP + Internet +
## Kat_Kegiatan, data = data)
##
## n= 1728, number of events= 80
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Klasifikasi -0.4344 0.6477 0.3117 -1.394 0.16340
## Kat_KIP -1.0846 0.3380 0.3962 -2.738 0.00619 **
## Internet -0.7263 0.4837 0.2400 -3.026 0.00248 **
## Kat_Kegiatan 2.0283 7.6008 0.2446 8.292 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Klasifikasi 0.6477 1.5440 0.3516 1.1930
## Kat_KIP 0.3380 2.9583 0.1555 0.7348
## Internet 0.4837 2.0674 0.3022 0.7742
## Kat_Kegiatan 7.6008 0.1316 4.7061 12.2762
##
## Concordance= 0.745 (se = 0.033 )
## Likelihood ratio test= 79.4 on 4 df, p=2e-16
## Wald test = 99.53 on 4 df, p=<2e-16
## Score (logrank) test = 137.5 on 4 df, p=<2e-16