library(survival)
data(ovarian)
## Warning in data(ovarian): data set 'ovarian' not found
head(ovarian) #untuk menampilkan 6 data pertama
##   futime fustat     age resid.ds rx ecog.ps
## 1     59      1 72.3315        2  1       1
## 2    115      1 74.4932        2  1       1
## 3    156      1 66.4658        2  1       2
## 4    421      0 53.3644        2  2       1
## 5    431      1 50.3397        2  1       1
## 6    448      0 56.4301        1  1       2

Keterangan : futime = waktu survival (terjadi event ataupun tidak) fustat = indikator yang menyatakan apakah futime teramati (fustat=1) atau tersensor (fustat=0) age = usia pasien (tahun) resid.ds = apakah masih terdapat residual disease (1 jika tdk, 2 jika ya) rx = treatment grup ecog.ps = status peforma ECOG, bernilai 1 jika kondisi lebih baik

Interpretasi: -Pasien1 berusia 72,3315 thn, terdapat residual disease (resid.ds=2), diberi perlakuan jenis 1 (rx=1), dengan kondisi ecog yg lebih baik (ecog.ps=1). Pasien mengalami event (fustat=1) pada waktu T=59(futime=59)

-Pasien6 berusia 53,4301 thn, tdk terdapat residual disease (resid.ds=1), diberi perlakuan jenis 1 (rx=1), dengan kondisi ecog tdk lebih baik (ecog.ps=2). Setelah selang waktu T=448 pasien blm mengalami event (fustat=0)

# Statistika Deskriptif
boxplot(futime~rx, data = ovarian, main = "Boxplot", xlab = "Treatment", 
        ylab = "Waktu Survival", 
        col=c("light blue","yellow"), border = "black")

#Estimasi regresi cox (serentak dan partial)
cox = coxph(Surv(futime, fustat)~rx+ecog.ps+resid.ds, method="breslow", data=ovarian)
summary(cox)
## Call:
## coxph(formula = Surv(futime, fustat) ~ rx + ecog.ps + resid.ds, 
##     data = ovarian, method = "breslow")
## 
##   n= 26, number of events= 12 
## 
##             coef exp(coef) se(coef)      z Pr(>|z|)  
## rx       -0.7493    0.4727   0.5949 -1.260   0.2078  
## ecog.ps   0.4528    1.5727   0.5903  0.767   0.4431  
## resid.ds  1.3466    3.8444   0.6803  1.980   0.0478 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## rx          0.4727     2.1156    0.1473     1.517
## ecog.ps     1.5727     0.6359    0.4945     5.002
## resid.ds    3.8444     0.2601    1.0134    14.584
## 
## Concordance= 0.741  (se = 0.074 )
## Likelihood ratio test= 6.03  on 3 df,   p=0.1
## Wald test            = 5.44  on 3 df,   p=0.1
## Score (logrank) test = 5.83  on 3 df,   p=0.1

Interpretasi: Terlihat bahwa koefisien untuk ketiga variabel, secara berurutan adalah -0.7493(beta1), 0.4528(beta2), 1.3466(beta3). Berdasarkan uji likelihood ratio, di peroleh p-value = 0.1, df sebesar 3 dan nilai likelihood 6,03. dimana nilai pvalue kurang dari alfa 10%, sehingga ddapat diartikan minimal ada satu beta yang signifikan terhadap model.

pada uji partial diketahui bahwa hanya variabel residual dieses yang signifikan dengan alfa 10%, maka perlu dilakukan estimasi parameter atau pemilihan model terbaik.

penulisan model model= 0,4727RX+1,5727Ecog.ps+3,8444Resid.ds

#pemilihan model terbaik

cox1 = coxph(Surv(futime, fustat)~resid.ds, method="breslow", data=ovarian)
cox1
## Call:
## coxph(formula = Surv(futime, fustat) ~ resid.ds, data = ovarian, 
##     method = "breslow")
## 
##            coef exp(coef) se(coef)     z      p
## resid.ds 1.2092    3.3507   0.6724 1.798 0.0721
## 
## Likelihood ratio test=3.76  on 1 df, p=0.05251
## n= 26, number of events= 12

Interpretasi: Terlihat bahwa koefisien untuk variabel residual deases 1.2092. Berdasarkan uji likelihood, di peroleh p-value = 0.0721, nilai tersebut kurang dari alfa 10%, sehingga diartikan bahwa variabel residual deases merupakan variabel yang mempengaruhi lama waktu pasien bertahan pada kanker ovarium

model = 1,2092 Resid.ds HR 3,3507

#menghasilkan residual Schoenfeld dan memeriksa asumsi PH:
model=coxph(Surv(futime, fustat)~rx+age+resid.ds, method="breslow", data=ovarian)
Asumsi=cox.zph(model)
Asumsi
##          chisq df    p
## rx       0.631  1 0.43
## age      0.210  1 0.65
## resid.ds 1.120  1 0.29
## GLOBAL   1.958  3 0.58

Interpretasi: Terlihat bahwa nilai p dari uji global (0.58) lebih besar dari nilai taraf sig. (0.05), sehingga dapat disimpulkan bahwa asumsi hazard terpenuhi.namun hasil estimasi ini masih dengan variabel yang tida signifikan.. sehingga dilakukan pengujian kembali

#menghasilkan residual Schoenfeld dan memeriksa asumsi PH:
model1=coxph(Surv(futime, fustat)~resid.ds, method="breslow", data=ovarian)
Asumsi=cox.zph(model1)
Asumsi
##          chisq df    p
## resid.ds   2.2  1 0.14
## GLOBAL     2.2  1 0.14

Interpretasi: Terlihat bahwa nilai p dari uji global (0.14) lebih besar dari nilai taraf sig. (0.05), sehingga dapat disimpulkan bahwa asumsi hazard terpenuhi.