Regresi Cox-PH
Alfisyahrina Hapsery
2023-10-04
install.packages(“KMsurv”)
# Memuat paket survival dan KMsurv
library(survival)
library(KMsurv)# Mengakses dataset larynx
data(larynx)
head(larynx)## stage time age diagyr delta
## 1 1 0.6 77 76 1
## 2 1 1.3 53 71 1
## 3 1 2.4 45 71 1
## 4 1 2.5 57 78 0
## 5 1 3.2 58 74 1
## 6 1 3.2 51 77 0??larynx #untuk mengetahui informasi data larynx## starting httpd help server ... done# Statistika Deskriptif
boxplot(time~stage, data = larynx, main = "Boxplot", xlab = "Stadium Penyakit",
ylab = "Waktu hingga meninggal (bulan))",
col=c("light blue","pink","orange","yellow"), border = "black")
# Fungsi Survival
fit.larynx <- survfit(Surv(time,delta)~stage, data=larynx, conf.type="log-log")
summary(fit.larynx)## Call: survfit(formula = Surv(time, delta) ~ stage, data = larynx, conf.type = "log-log")
##
## stage=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.6 33 1 0.970 0.0298 0.804 0.996
## 1.3 32 1 0.939 0.0415 0.779 0.984
## 2.4 31 1 0.909 0.0500 0.744 0.970
## 3.2 29 1 0.878 0.0573 0.706 0.952
## 3.3 27 1 0.845 0.0637 0.667 0.933
## 3.5 25 2 0.778 0.0744 0.588 0.888
## 4.0 23 2 0.710 0.0819 0.515 0.838
## 4.3 21 1 0.676 0.0847 0.481 0.811
## 5.3 18 1 0.639 0.0879 0.441 0.782
## 6.0 14 1 0.593 0.0927 0.391 0.748
## 6.4 11 1 0.539 0.0987 0.331 0.708
## 6.5 10 1 0.485 0.1025 0.277 0.665
## 7.4 6 1 0.404 0.1129 0.191 0.610
##
## stage=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.2 17 1 0.941 0.0571 0.650 0.991
## 1.8 16 1 0.882 0.0781 0.606 0.969
## 2.0 15 1 0.824 0.0925 0.547 0.939
## 3.6 11 1 0.749 0.1103 0.456 0.899
## 4.0 9 1 0.665 0.1255 0.364 0.849
## 6.2 5 1 0.532 0.1557 0.209 0.776
## 7.0 4 1 0.399 0.1641 0.110 0.683
##
## stage=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.3 27 2 0.926 0.0504 0.735 0.981
## 0.5 25 1 0.889 0.0605 0.694 0.963
## 0.7 24 1 0.852 0.0684 0.652 0.942
## 0.8 23 1 0.815 0.0748 0.611 0.918
## 1.0 22 1 0.778 0.0800 0.571 0.893
## 1.3 21 1 0.741 0.0843 0.532 0.867
## 1.6 20 1 0.704 0.0879 0.494 0.839
## 1.8 19 1 0.667 0.0907 0.457 0.811
## 1.9 18 2 0.593 0.0946 0.386 0.750
## 3.2 16 1 0.556 0.0956 0.352 0.718
## 3.5 15 1 0.519 0.0962 0.319 0.685
## 5.0 10 1 0.467 0.0995 0.267 0.644
## 6.3 7 1 0.400 0.1053 0.200 0.593
## 6.4 6 1 0.333 0.1068 0.143 0.538
## 7.8 4 1 0.250 0.1078 0.078 0.471
##
## stage=4
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.1 13 1 0.923 0.0739 0.56636 0.989
## 0.3 12 1 0.846 0.1001 0.51220 0.959
## 0.4 11 1 0.769 0.1169 0.44214 0.919
## 0.8 10 2 0.615 0.1349 0.30834 0.818
## 1.0 8 1 0.538 0.1383 0.24766 0.760
## 1.5 7 1 0.462 0.1383 0.19161 0.696
## 2.0 6 1 0.385 0.1349 0.14054 0.628
## 2.3 5 1 0.308 0.1280 0.09498 0.554
## 3.6 3 1 0.205 0.1196 0.03845 0.463
## 3.8 2 1 0.103 0.0940 0.00666 0.355plot(fit.larynx, lty = 1:4, col = 1:4,
main = expression(paste("Estimasi Kaplan-Meier", hat(S)(t))),
xlab = "Waktu survival,T(bulan)", ylab = "Fungsi Survival, S(t)")
legend(8,1,c("Stage 1", "Stage 2", "Stage 3", "Stage 4"), lty = 1:4, col = 1:4,
bty="n")
survdiff(Surv(time,delta)~stage, data=larynx)## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=1 33 15 22.57 2.537 4.741
## stage=2 17 7 10.01 0.906 1.152
## stage=3 27 17 14.08 0.603 0.856
## stage=4 13 11 3.34 17.590 19.827
##
## Chisq= 22.8 on 3 degrees of freedom, p= 5e-05# Fungsi Survival
fit.larynx <- survfit(Surv(time,delta)~stage, data=larynx, conf.type="log-log")
summary(fit.larynx)## Call: survfit(formula = Surv(time, delta) ~ stage, data = larynx, conf.type = "log-log")
##
## stage=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.6 33 1 0.970 0.0298 0.804 0.996
## 1.3 32 1 0.939 0.0415 0.779 0.984
## 2.4 31 1 0.909 0.0500 0.744 0.970
## 3.2 29 1 0.878 0.0573 0.706 0.952
## 3.3 27 1 0.845 0.0637 0.667 0.933
## 3.5 25 2 0.778 0.0744 0.588 0.888
## 4.0 23 2 0.710 0.0819 0.515 0.838
## 4.3 21 1 0.676 0.0847 0.481 0.811
## 5.3 18 1 0.639 0.0879 0.441 0.782
## 6.0 14 1 0.593 0.0927 0.391 0.748
## 6.4 11 1 0.539 0.0987 0.331 0.708
## 6.5 10 1 0.485 0.1025 0.277 0.665
## 7.4 6 1 0.404 0.1129 0.191 0.610
##
## stage=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.2 17 1 0.941 0.0571 0.650 0.991
## 1.8 16 1 0.882 0.0781 0.606 0.969
## 2.0 15 1 0.824 0.0925 0.547 0.939
## 3.6 11 1 0.749 0.1103 0.456 0.899
## 4.0 9 1 0.665 0.1255 0.364 0.849
## 6.2 5 1 0.532 0.1557 0.209 0.776
## 7.0 4 1 0.399 0.1641 0.110 0.683
##
## stage=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.3 27 2 0.926 0.0504 0.735 0.981
## 0.5 25 1 0.889 0.0605 0.694 0.963
## 0.7 24 1 0.852 0.0684 0.652 0.942
## 0.8 23 1 0.815 0.0748 0.611 0.918
## 1.0 22 1 0.778 0.0800 0.571 0.893
## 1.3 21 1 0.741 0.0843 0.532 0.867
## 1.6 20 1 0.704 0.0879 0.494 0.839
## 1.8 19 1 0.667 0.0907 0.457 0.811
## 1.9 18 2 0.593 0.0946 0.386 0.750
## 3.2 16 1 0.556 0.0956 0.352 0.718
## 3.5 15 1 0.519 0.0962 0.319 0.685
## 5.0 10 1 0.467 0.0995 0.267 0.644
## 6.3 7 1 0.400 0.1053 0.200 0.593
## 6.4 6 1 0.333 0.1068 0.143 0.538
## 7.8 4 1 0.250 0.1078 0.078 0.471
##
## stage=4
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.1 13 1 0.923 0.0739 0.56636 0.989
## 0.3 12 1 0.846 0.1001 0.51220 0.959
## 0.4 11 1 0.769 0.1169 0.44214 0.919
## 0.8 10 2 0.615 0.1349 0.30834 0.818
## 1.0 8 1 0.538 0.1383 0.24766 0.760
## 1.5 7 1 0.462 0.1383 0.19161 0.696
## 2.0 6 1 0.385 0.1349 0.14054 0.628
## 2.3 5 1 0.308 0.1280 0.09498 0.554
## 3.6 3 1 0.205 0.1196 0.03845 0.463
## 3.8 2 1 0.103 0.0940 0.00666 0.355plot(fit.larynx, lty = 1:4, col = 1:4,
main = expression(paste("Estimasi Kaplan-Meier", hat(S)(t))),
xlab = "Waktu survival,T(bulan)", ylab = "Fungsi Survival, S(t)")
legend(8,1,c("Stage 1", "Stage 2", "Stage 3", "Stage 4"), lty = 1:4, col = 1:4,
bty="n")
survdiff(Surv(time,delta)~stage, data=larynx)## Call:
## survdiff(formula = Surv(time, delta) ~ stage, data = larynx)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## stage=1 33 15 22.57 2.537 4.741
## stage=2 17 7 10.01 0.906 1.152
## stage=3 27 17 14.08 0.603 0.856
## stage=4 13 11 3.34 17.590 19.827
##
## Chisq= 22.8 on 3 degrees of freedom, p= 5e-05# Membuat Model Regresi Cox:
reg.larynx <- coxph(Surv(time,delta)~factor(stage), data=larynx)
summary(reg.larynx)## Call:
## coxph(formula = Surv(time, delta) ~ factor(stage), data = larynx)
##
## n= 90, number of events= 50
##
## coef exp(coef) se(coef) z Pr(>|z|)
## factor(stage)2 0.06481 1.06696 0.45843 0.141 0.8876
## factor(stage)3 0.61481 1.84930 0.35519 1.731 0.0835 .
## factor(stage)4 1.73490 5.66838 0.41939 4.137 3.52e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## factor(stage)2 1.067 0.9372 0.4344 2.62
## factor(stage)3 1.849 0.5407 0.9219 3.71
## factor(stage)4 5.668 0.1764 2.4916 12.90
##
## Concordance= 0.668 (se = 0.037 )
## Likelihood ratio test= 16.49 on 3 df, p=9e-04
## Wald test = 19.24 on 3 df, p=2e-04
## Score (logrank) test = 22.88 on 3 df, p=4e-05