Pada postingan kali ini saya ingin memberikan sedikit ilustrasi tentang analisis survival yakni:
1. Pengetahuan tentang dasar-dasar analisis survival,
2. Membiasakan diri dengan data frame, grafik, dan model linear di R, dan
3. Mencoba untuk menerapkan analisis survival menggunakan R.
Saya akan menggunakan R versi terbaru untuk saat ini yaitu (R.3.1.3), sehingga bagi anda yang ingin mencoba atau ingin mensimulasikan program yang terlampir pada tulisan ini dapat menggunakan R versi yang sama atau versi terbaru. Dalam hal ini saya juga akan menggunakan paket (package) yang sudah ada di dalam R itu sendiri seperti halnya:
survival, OIsurv, dan paket KMsurv
Paket "survival" akan digunakan dalam setiap contoh dalam dokumen ini dan paket data yang digunakan dapat ditemukan pada "KMsurv", yang meliputi paket data dari buku ke 5 oleh Klein dan Moeschberger. Fungsi tambahan yang digunakan dapat ditemukan pada "OIsurv". Paket-paket ini dapat diinstal dengan menggunakan install.packages dalam R sebagai berikut:
> install.packages("OIsurv")
Dengan Instalasi ini akan secara otomatis memuat tiga paket termasuk "survival" dan "KMsurv". Untuk memuat atau mempersiakan paket sehingga dapat digunakan, lakukan pemanggilan paket
> library(OIsurv) # pemanggilan paket
Loading required package: survival
Loading required package: splines
Loading required package: KMsurv
Untuk melihat paket data yang tersedia dalam "KMsurv" anda dapat mengunjungi http://cran.r-project.org/web/packages/KMsurv/KMsurv.pdf . Andaikan anda ingin menggunakan salah satu data tersebut, anda dapat mengunakan fungsi "data(data yang anda inginkan)". Sebagai contoh,
> data(aids)
> aids # menampilkan seluruh data aids
> head(aids) # menampilkan data aids paling depan saja
infect induct adult
1 0.00 5.00 1
2 0.25 6.75 1
3 0.75 5.00 1
4 0.75 5.00 1
5 0.75 7.25 1
6 1.00 4.25 1
Selanjutnya, membuat kolom data yang tersedia untuk digunakan sebagai variabel adalah,
> attach(aids) # membuat kolom data sebagai variabel
> infect # menampilkan variabel infect
[1] 0.00 0.25 0.75 0.75 0.75 1.00 1.00 1.00 1.00 1.25 1.25 1.25 1.25 1.50 1.50 .... dan seterusnya.
Sebagai programer yang baik anda juga dapat membuat kolom data yang tersedia untuk digunakan sebagai variabel dengan,
> aids$infect
[1] 0.00 0.25 0.75 0.75 0.75 1.00 1.00 1.00 1.00 1.25 1.25 1.25 1.25 1.50 1.50 .... dan seterusnya.
> detach(aids) # mengembalikan data(aids) ke data awal sebelum observasi dilakukan. Sebagai programer yang baik ada baiknya anda lakukan hal ini di akhir observasi.
> detach(aids) # mengembalikan data(aids) ke data awal sebelum observasi dilakukan. Sebagai programer yang baik ada baiknya anda lakukan hal ini di akhir observasi.
Target Survival:
Surv(time,event), Surv(time,time2,event, type)
Banyak fungsi dalam paket R untuk mengaplikasikan analisis survival seperti fungsi "Surv ()". Kita akan bahas data survival tersensor dikanan (right censor), data yang belum adanya suatu kejadian baik\buruk pada masa obsevasi usai. Lihat referensi http://www.weibull.com/LifeDataWeb/data_classification.htm untuk deskripsi tipe data survival. Untuk data tersensor dikanan, ada dua argumen yang perlu diperhatikan dalam penggunaan fungsi Surv (): dengan menunjukkan vektor waktu yang diamati dan disensor.
> data(tongue)
> attach(tongue)
> head(tongue)
type time delta
1 1 1 1
2 1 3 1
3 1 3 1
4 1 4 1
5 1 10 1
6 1 13 1
Menentukan himpunan bagian dari kelompok pertama dengan menggunakan [tipe == 1]
> my.surv.object<-Surv(time[type==1],delta[type==1])
> my.surv.object
[1] 1 3 3 4 10 13 13 16 16 24 26 27 28 30 30
[16] 32 41 51 65 67 70 72 73 77 91 93 96 100 104 157
[31]167 61+ 74+ 79+ 80+ 81+ 87+ 87+ 88+ 89+ 93+ 97+ 101+ 104+ 108+
[46]109+ 120+ 131+ 150+ 231+ 240+ 400+
> detach(tongue)
Untuk mengidentifikasi pengamatan yang data tersensor dikanan. Argumen pertama di "Surv()" harus dimasukan sebagai vektor waktu pengamatan dan tersensor dikanan. Vektor indikator yang digunakan dalam argumen kedua untuk menandakan apakah kejadian itu diamati (1) atau tidak (0). Seperti yang anda perhatikan diatas bahwa data yang bertanda (+) adalah data yang tersensor di kanan (hasil belum deketahui saat periode penelitian atau observasi di hentikan).
Anda juga dapat mempertimbangkan data yang terpotong dan tersensor dikanan tetapi dilakukan observasi dibagian kiri perhatikan contoh berikut:
> data(psych)
> attach(psych)
> head(psych)
sex age time death
1 2 51 1 1
2 2 58 1 1
3 2 55 2 1
4 2 28 22 1
5 1 21 30 0
6 1 19 28 1
> my.surv.object<-Surv(age,age+time,death)
> my.surv.object
[1](51,52] (58,59] (55,57] (28,50] (21,51+] (19,47] (25,57] (48,59] (47,61] (25,61+] (31,62+] (24,57+] (25,58+] (30,67+] (33,68+] (36,61] (30,61+] (41,63] (43,69] (45,69] (35,70+] (29,63+] (35,65+] (32,67] (36,76] (32,71+]
> detach(psych)
Catatan: Ada banyak jenis analisis survival yang bisa anda kreasikan, tetapi mereka tidak tercakup dalam tulisan ini. Selain itu, beberapa fungsi survival di R hanya menerima beberapa jenis data survival.
Estimasi Kaplan-Meier dan batas terbaik (pointwise interval):
survfit(formula, conf.int = 0.95, conf.type = "log")
Estimasi Kaplan-Meier adalah memprediksi maksimum kemungkinan (maximum likelihood estimate) nonparametrik dari fungsi survival, $S(t).$ Estiamasi ini adalah langkah untuk memprediksi waktu yang diamati, $t_i.$ Dalam formula matematika di bawah ini, diasumsikan $t_i$ sebagai: $0 <t_1 <t_2 < \dots < t_D$. Jika urutan individu dinyatakan dengan waktu $t_i$ pengamatan adalah $di$, dan nilai $Y_i $ merupakan jumlah individu yang berisiko saat $t_i$ (di mana risiko berarti orang yang meninggal pada waktu $t_i$ atau lambat), sehingga perkiraan Kaplan-Meier adalah fungsi survival dan estimasi varians dinyatakan dengan
Batas kekauratan (Confidence Interval) untuk linear sederhana dan "log-log" yang juga disediakan dalam R adalah,
Estimasi Kaplan-Meier adalah "fit" di R dengan menggunakan fungsi "survfit ()". Memprediksi atau estimasi batas kelayakan adalah dengan mengimput fungsi survival terhadap intercept:
> data(tongue)
> attach(tongue)
> my.surv<-Surv(time[type==1],delta[type==1])
> survfit(my.surv~1)
Call: survfit(formula = my.surv ~ 1)
records n.max n.start events median 0.95LCL 0.95UCL
52 52 52 31 93 67 NA
Survfit () juga memiliki sejumlah argumen opsional. Sebagai contoh, tingkat keakuratan memungkinkan untuk diubah (misalnya conf.int = 0.90 untuk 90% batas keakuratan) dengan menggunakan argumen, "conf.int ". Argumen conf.int menjelaskan jenis titik keakuratan. Lebih spesifiknya lagi, hal itu menggambarkan transformasi untuk membangun titik keakuratan. Standarnya adalah "log" $g (t) = log (t),$ yang setara untuk fungsi transformasi. Selain itu anda juga dapat menggunkan fungsi "log-log" $g (t) = log (-log (t))$ sebagai pilihan anda mentrasformasikan kembali fungsi tersebut.
Seperti banyak fungsi dalam R, survfit () adalah salah satu cara untuk mengembalikan fungsi informasi tersembunyi yang dapat diakses. Berikut diperlihatkan beberapa elemen informasi tersembunyi, yang disimpan dalam daftar fungsi survfit (). Untuk ringkasan lengkap dari objek, mengunakan fungsi "str" untuk my.fit dan "summary()".
> my.fit <- survfit(my.surv)
> summary(my.fit)$surv # estimasi Kaplan-Meier pada setiap t_i
> summary(my.fit)$time # {t_i}
> summary(my.fit)$n.risk # {Y_i}
> summary(my.fit)$n.event # {d_i}
> summary(my.fit)$std.err # standar kesalahan K-M estimasi {t_i}
> summary(my.fit)$lower # batas bawah estimasi
> str(my.fit) # ringakasan lengkap {my.fit}
> str(summary(my.fit)) # ringakasan lengkap {my.fit}
Fungsi "str" berguna untuk melihat lebih rincian tentang apa yang terkandung dalam daftar fungsi, dan seperti yang dinyatakan di atas, kita dapat mengakses setiap item dalam daftar menggunakan operator $. Estimasi Kaplan-Meier dapat dinyatakann dalam grafik menggunakan fungsi "plot(data yang sudah di estimasi)". Argumen standar dalam fungsi plot dapat digunakan untuk meningkatkan estetika (karakter) grafik:
> plot(my.fit,main="Kaplan-Meier estimate with 95% confidence bounds",xlab="time",ylab="survival function")
Kadang-kadang perbedaan kelompok yang terkandung dalam sebuah objek tunggal survival. Misalnya, jenis variabel dalam kumpulan data "tongue" menggambarkan propil DNA pasien. Kita dapat memperoleh estimasi Kaplan-Meier untuk masing-masing kelompok dengan regresi objek survival pada jenis variabel:
> my.fit1<-survfit(Surv(time,delta)~type) #fit dengan pemisahan type
> detach(tongue)
Dengan menggunakan code ini kita dapat mempartisi data set dengan memperhatikan type data tersebut. (Hal ini menjadi alasan untuk menggunakan beberapa variabel di sisi kanan persamaan.) Ringkasan my.fit1 akan berisi daftar item tambahan-strata, dapat diakses melalui summary statistik (my.fit1)$strata yang menunjukkan komponen output sesuai dengan kelompok lain.
Interval Keakuratan Simultan (Simultaneous Confidence Interval)
confBands(x, confType="plain", confLevel=0.9, type="ep")
Interval kekakuratan, seperti yang diperkenalkan di halaman sebelumnya, berlaku untuk satu titik dalam skala waktu. Sekarang kita mengalihkan perhatian kita keakuratan simultan (untuk jangka pendek atau kurun waktu tertentu), yang berlaku untuk seluruh rentang nilai waktu secara bersamaan. Suatu keakuratan 95% misalnya, akan menangkap seluruh kurva survival yang benar sekitar 19 dari 20 kali.
Sementara paket yang digunakan fungsi survival tidak memenuhi batasan, jadi dilakukan perhitungan dengan menggunakan "confBands" dari perpustakaan "OIsurv". Argumen pertama adalah objek survival untuk $x,$ dan argumen lainnya memungkinkan kustomisasi. kode untuk "confType" adalah, "log-log", atau "asin-sqrt"; dengan confLevel mungkin 0,90, 0,95, atau 0,99; dan jenisnya mungkin "ep" atau "hall" (Hall-Wellner). Ada juga dua argumen opsional, tL (batas bawah) dan tU (batas atas), yang membatasi dukungan dari batasan keakuratan ini. Batasan kekuaratan dapat ditambahkan ke plot menggunakan fungsi garis.
> data(bmt); attach(bmt)
> my.surv <- Surv(t2[group==1], d3[group==1])
> my.cb<-confBands(my.surv,confLevel=0.95,type="hall")
>plot(survfit(my.surv~1),xlim=c(100,600),xlab="waktu",ylab="Estimasi fungsi Survival",main="Batas keakuratan")
Perhatikan bahwa hasil dari batas keakuratan (atas/bawah) adalah transformasi dari log dengan melanjutkan kode R dibawah ini, yang akan menghasilkan interval batas keakuratan asimetris.
> lines(my.cb$time,my.cb$lower,lty=3,type="s")
> legend(100,0.3,legend=c("Etimasi survival K-M","batas atas","batas bawah"),lty=1:3)
> detach(bmt)
Fungsi Kumulatif Resiko (Hazard Function)
Fungsi Kumulatif Resiko dan fungsi kelangsungan hidup (survival) terkait dengan cara berikut untuk Data kontinu:
Estimasi maksimal kemungkinan (MLE) fungsi resiko dapat diperoleh dengan transformasi kebalikan dari Kaplan-Meier Estimasi:
Anda juga dapat menggunakan cara lain untuk memperkirakan H (t) seperti estimator Nelson-Aalen:
Sementara tidak ada fungsi dalam menghitung paket survival secara otomatis, pengembalian objek dengan ringkasan (survfit ()) dapat digunakan untuk menghitung estimasi:
> data(tongue); attach(tongue)
> my.surv<-Surv(time[type==1],delta[type==1])
> my.fit<-summary(survfit(my.surv~1))
> my.fit
Call: survfit(formula = my.surv ~ 1)
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 52 1 0.981 0.0190 0.944 1.000
3 51 2 0.942 0.0323 0.881 1.000
4 49 1 0.923 0.0370 0.853 0.998
10 48 1 0.904 0.0409 0.827 0.988
13 47 2 0.865 0.0473 0.777 0.963
.......... dan seterusnya
> H.hat <- -log(my.fit$surv)
> H.hat <- c(H.hat, tail(H.hat, 1))
Suatu ringkasan plot atau tabel dapat dibuat dengan menggunakan H.hat dengan my.fit$waktu. Estimasi Nelson-Aalen juga dapat dibangun dengan:
> h.sort.of<-my.fit$n.event/my.fit$n.risk
> H.tilde<-cumsum(h.sort.of)
> H.tilde<-c(H.tilde,tail(H.tilde,1))
> plot(c(my.fit$time,250),H.hat,xlab="waktu",ylab="Resiko kumulatif resiko",main="Perbandingan resiko kumulatif",ylim=range(c(H.hat, H.tilde)),type="s")
> points(c(my.fit$time,250),H.tilde,lty=2,type="s",col=2)
>
> detach(tongue)
Perkiraan rata-rata dan variansi
Median waktu survival adalah waktu $t_{0.5}$ sehingga $S (t_{0. 5}) = 0. 5.$ Hal ini divisualisasikan oleh grafik perkiraan fungsi survival dan menggambar garis horizontal di $0.5.$ Diperkirakan median sama waktu $\hat{ t}_{0.5}$ di mana fungsi dan garis berpotongan. Interval con dence untuk $t_{0. 5}$ diberikan oleh titik di mana garis horizontal ini menyilang interval batas keakuratan $\hat{ S} (t).$
Rata-rata waktu survival dapat di estimasi dengan
$$\mu=\int_0^\infty S(t)dt, \quad \hat{\mu}=\int_0^\infty\hat{S}(t)dt$$
Jika $S(t)$ atau $\hat{S}(t)$ tidak konvergen di nol, maka integral divergen. Properti ini menciptakan tantangan untuk banyak data, dan satu resolusi adalah dengan menggunakan nilai terbatas $\tau$ sebagai terikat untuk integral, di mana $\tau$ dapat mewakili waktu survival maksimum dianggap memungkinkan. Pilihan lain yang masuk akal untuk $\tau$ adalah maksimum yang dapat diamati atau disensor waktu. Menggunakan nite $\tau$ hasil dalam statistik baru dan sesuai estimasi: $\mu_\tau= \int_0^\tau S (t)dt.$ Andaikan $t_i, Y_i, d_i,$ dan $D$ adalah sebagaimana dijelaskan dalam Kaplan-Meier, estimasi varians $\hat{\mu}_\tau$ adalah
$$\hat{V}(\hat{\mu}_\tau)=\sum_{i=1}^D\biggr[\int_{t_i}^\tau\hat{S}(t)dt\biggl]^2{d_i\over Y_i(Y_i-d_i)}$$
Median dan 95% interval keakuratan dapat diestimasi dengan fungsi "survfit()":
> data(drug6mp)
> attach(drug6mp)
> head(drug6mp)
pair remstat t1 t2 relapse
1 1 1 1 10 1
2 2 2 22 7 1
3 3 2 3 32 0
4 4 2 12 23 1
5 5 2 8 22 1
6 6 1 17 6 1
> my.surv<-Surv(t1,rep(1,21)) # semua pasien placebo dioservasi
> survfit(my.surv~1)
Call: survfit(formula=my.surv~1)
records n.max n.start events median 0.95LCL 0.95UCL
21 21 21 21 8 4 12
Menggunakan "print(survfit())" akan menampilkan waktu survival dan standar eror.
> data(drug6mp)
> attach(drug6mp)
> head(drug6mp)
pair remstat t1 t2 relapse
1 1 1 1 10 1
2 2 2 22 7 1
3 3 2 3 32 0
4 4 2 12 23 1
5 5 2 8 22 1
6 6 1 17 6 1
> my.surv<-Surv(t1,rep(1,21)) # semua pasien placebo dioservasi
> survfit(my.surv~1)
Call: survfit(formula=my.surv~1)
records n.max n.start events median 0.95LCL 0.95UCL
21 21 21 21 8 4 12
Menggunakan "print(survfit())" akan menampilkan waktu survival dan standar eror.
> print(survfit(my.surv~1),print.rmean=TRUE)
Call: survfit(formula=my.surv~1)
records n.max n.start events *rmean *se(rmean) median
21.00 21.00 21.00 21.00 8.67 1.38 8.00
0.95LCL 0.95UCL
4.00 12.00
* restricted mean with upper limit = 23
Call: survfit(formula=my.surv~1)
records n.max n.start events *rmean *se(rmean) median
21.00 21.00 21.00 21.00 8.67 1.38 8.00
0.95LCL 0.95UCL
4.00 12.00
* restricted mean with upper limit = 23
> detach(drug6mp)
Kode "print.rmean =TRUE" adalah argumen yang digunakan untuk mendapatkan mean dan standard error-nya, dan $\tau$ secara otomatis ditetapkan sebagai terbesar diamati atau disensor waktu. Atau, $\tau$ ditentukan menggunakan argumen "rmean".
Uji dengan dua atau lebih sampel
Diberikan dua atau lebih sampel, apakah ada diferinsial statistik antara waktu survival? kita dapat melakukan analisis dengan formula (hipotesis fungsi resiko) sebagai,
Diberikan dua atau lebih sampel, apakah ada diferinsial statistik antara waktu survival? kita dapat melakukan analisis dengan formula (hipotesis fungsi resiko) sebagai,
- $H_0 : h_1(t)=h_2(t)=\cdots=h_n(t)$ untuk semua $ t,$ dan
- $H_A : h_i(t_0)\ne h_j(t_0)$ untuk setidaknya satu pasangan $i; j$ dan waktu $t_0$
- $t_i$ adalah waktu di mana peristiwa yang diamati (menganggap $D$ sebagai waktu),
- $d_{ik}$ nilai kejadian obbservasi dari grup $k$ pada waktu $t_i$,
- $Y_{ik}$ nilai subjek dalam grup $k$ terdapat resiko pada $t_i$ (lihat K-M section),
- $d_i=\sum_{j=1}^n d_{ij},$
- $Y_i=\sum_{j=1}^n Y_{ij},$ dan
- $W(t_i)$ menjadi titik berat pengamatan pada waktu $t_i.$
Kemudian Uji hipotesis diatas, suatu kalkulasi vector $Z$, dimana elemen $k^{th}$ adalah $$Z_k=\sum_{i=1}^D W(t_i)\biggr[d_{ik}-Y_{ik} {d_i\over Y_i}\biggl]$$
Kovarians matriks $\hat{\sum}$ juga perhitungan data. Di bawah nol hipotesis uji statistik $X^2=Z'\hat{\sum}^{-1}Z$ melalui distribusi $\chi^2$ dengan $n$ derajat kebebasan. Dalam hal ini, jika $X^2>\chi_{1-\alpha,df-n}^2,$ data ini memberikan bukti kuat terhadap nol hipotesis dan kita menolak $H_0.$
Fungsi "survdiff" digunakan dalam nol hipotesis ini. Argumen pertama objek survival suatu kategori variabel kovarians adalah biasanya menunjuk variabel kelompok yang sesuai dengan waktu survival. Objek dikembalikan oleh "survdiff ()" memberikan informasi yang berguna.
> data(btrial); attach(btrial)
> head(btrial)
time death im
1 19 1 1
2 25 1 1
3 30 1 1
4 34 1 1
5 37 1 1
6 46 1 1
Fungsi "survdiff" digunakan dalam nol hipotesis ini. Argumen pertama objek survival suatu kategori variabel kovarians adalah biasanya menunjuk variabel kelompok yang sesuai dengan waktu survival. Objek dikembalikan oleh "survdiff ()" memberikan informasi yang berguna.
> data(btrial); attach(btrial)
> head(btrial)
time death im
1 19 1 1
2 25 1 1
3 30 1 1
4 34 1 1
5 37 1 1
6 46 1 1
> survdiff(Surv(time,death)~im) # hasil
Call:
survdiff(formula = Surv(time, death) ~ im)
N Observed Expected (O-E)^2/E (O-E)^2/V
im=1 36 16 20.19 0.869 5.49
im=2 9 8 3.81 4.599 5.49
Chisq= 5.5 on 1 degrees of freedom, p= 0.0191
Call:
survdiff(formula = Surv(time, death) ~ im)
N Observed Expected (O-E)^2/E (O-E)^2/V
im=1 36 16 20.19 0.869 5.49
im=2 9 8 3.81 4.599 5.49
Chisq= 5.5 on 1 degrees of freedom, p= 0.0191
Kode kedua "survdiff ()", rho, menunjuk bobot menurut $\hat{S} (t)^\rho$ dan mungkinkan untuk nilai numerik. Standarnya adalah rho = 0 dan sesuai dengan uji log-rank. Peto & Peto modi fikasi dari perhitungan uji Gehan-Wilcoxon dengan menggunakan rho=1:
> survdiff(Surv(time,death)~im,rho=1)
Call:
survdiff(formula=Surv(time,death)~im,rho=1)
N Observed Expected (O-E)^2/E (O-E)^2/V
im=1 36 12.1 14.96 0.558 4.35
im=2 9 5.8 2.91 2.867 4.35
Chisq= 4.4 on 1 degrees of freedom, p= 0.037
> detach(btrial)
Untuk memberikan bobot yang lebih besar untuk pertama bagian dari kurva survival, gunakan rho lebih besar dari 0. Untuk memberikan bobot ke bagian akhir dari kurva survival, gunakan rho lebih kecil dari 0. Output dari "survdiff()" adalah umumnya cukup jelas. Suatu $\chi^2$ statistik dihitung bersama dengan nilai-p.
Model resiko proposional Cox , konstanta kovariat
coxph(formula, method)
Model resiko proporsional Cox (Cox PH) kelayakan survival (fit) data dengan kovariat $z $ untuk fungsi resiko dibentuk
$$h(t|s)=h_0(t)exp\{\beta'z\}$$
diman $\beta$ suatu vektor yang tidak diketahui dan $h_0(t)$ adalah resiko baseline nonparametrik. Hal medasar untuk memperkirakan parameter $\beta$ menggunakan kemungkinan parsial:
> survdiff(Surv(time,death)~im,rho=1)
Call:
survdiff(formula=Surv(time,death)~im,rho=1)
N Observed Expected (O-E)^2/E (O-E)^2/V
im=1 36 12.1 14.96 0.558 4.35
im=2 9 5.8 2.91 2.867 4.35
Chisq= 4.4 on 1 degrees of freedom, p= 0.037
> detach(btrial)
Untuk memberikan bobot yang lebih besar untuk pertama bagian dari kurva survival, gunakan rho lebih besar dari 0. Untuk memberikan bobot ke bagian akhir dari kurva survival, gunakan rho lebih kecil dari 0. Output dari "survdiff()" adalah umumnya cukup jelas. Suatu $\chi^2$ statistik dihitung bersama dengan nilai-p.
Model resiko proposional Cox , konstanta kovariat
coxph(formula, method)
Model resiko proporsional Cox (Cox PH) kelayakan survival (fit) data dengan kovariat $z $ untuk fungsi resiko dibentuk
$$h(t|s)=h_0(t)exp\{\beta'z\}$$
diman $\beta$ suatu vektor yang tidak diketahui dan $h_0(t)$ adalah resiko baseline nonparametrik. Hal medasar untuk memperkirakan parameter $\beta$ menggunakan kemungkinan parsial:
dimana $R(t_i)$ adalah himpunan resiko pada $t_i$
Kemugkinan resiko terbesar (MLE) $\hat{\beta}$ sebagai suatu vektor adalah asimtotik $N(\beta,I^{-1})$, dimana $I$ adalah menyatakan informasi dari Fihser. Aproksimasi normal ini memungkinkan melakukan uji lokal. Uji lokal adalah memeriksa himpunan bagian unsur-unsur yang terkandung dalam $\beta$, pastikan uji dendan $C\beta=d$ dimana $C$ adalah matriks $q\times p$ dengan semua urutan dan $d$ adalah suatu vektor dengan panjang $q.$ Sebagai contoh anda ingin mencoba uji sederhana untuk $\beta_1=0$ dengan memilih $C_{1 \times p} = (1,0,0,\dots,0)$ dan $d_{1\times 1}=0.$ Secara umum $C$ dan $d$ uji statistik adalah,
$$X_W^2=\biggl(C\hat{\beta}-d\biggr)'\biggl[C\hat{I^{-1}}C'\biggr]^{-1}\biggl(C\hat{\beta}-d\biggr), $$
dimana nol hipotesis mengguanakan $\chi_q^2.$ Hal ini dikenal sebagai uji Wald. Selain memperoleh uji nilai-p, mungkin ada kepentingan dalam fungsi survival untuk kovariat tertentu. Jika perkiraan fungsi survival dasar $\hat{ S}_0 (t)$ diberikan, maka estimasi fungsi survival bagi suatu individu (obsevasi) dengan kovariat $z_k$ dapat diperoleh
$$\hat{S}(t|z_k)=\biggl[\hat{ S}_0 (t)\biggr]^{exp(\hat{\beta'}z_k)}$$
Fungsi "coxph ()" menentukan model terbaik Cox PH ke data yang diberikan. Argumen pertama adalah formula, di mana respon adalah objek survival.
> data(burn); attach(burn)
> my.surv<-Surv(T1,D1)
> coxph.fit<-coxph(my.surv~Z1+as.factor(Z11),method="breslow")
> coxph.fit
Call:
coxph(formula=my.surv~Z1+as.factor(Z11),method="breslow")
coef exp(coef) se(coef) z p
Z1 0.497 1.644 0.208 2.38 0.017
as.factor(Z11)2 -0.877 0.416 0.498 -1.76 0.078
as.factor(Z11)3 -1.650 0.192 0.802 -2.06 0.040
as.factor(Z11)4 -0.407 0.666 0.395 -1.03 0.300
Likelihood ratio test=14.6 on 4 df, p=0.00569 n= 154, number of events= 99
Dua kovariat telah digunakan dalam contoh ini. Argumen kedua yang terdaftar, metode spesifik bagaimana ikatan ditangani dengan defaultnya adalah "efron" dan pilihan lain "Breslow" dan "exact". banyak informasi yang berguna diperoleh dalam ringkasan "coxph ()", termasuk
- Memperkiraan $\beta_k$, termasuk standar kesalahan (standar diviasi) dan nilai-p untuk setiap uji $H_0: \beta_k = 0,$
- Meperkirakan rasio risiko dan selang confidence nya, dan
- Nilai-p untuk rasio kemungkinan, Wald, dan nilai uji untuk nol hipotesis secara umum, $H_0: \beta_i = 0$ untuk semua $i.$
> co<-coxph.fit$coefficients # menampilkan koefisien model
> va<-coxph.fit$var # estimasi kovarians matriks
> va
[,1] [,2] [,3] [,4]
[1,] 4.341949e-02 0.007746416 0.0009624982 6.528035e-05
[2,] 7.746416e-03 0.248005461 0.1437885487 1.442857e-01
[3,] 9.624982e-04 0.143788549 0.6439844004 1.435196e-01
[4,] 6.528035e-05 0.144285704 0.1435196176 1.563852e-01
> ll<-coxph.fit$loglik # estimasi log-likelihood (MLEs)
> ll
[1] -423.3526 -416.0686
Untuk mendapatkan fungsi dasar survival dari model PH Cox, gunakan survfit () untuk coxph ():
> my.survfit.objek<- survfit(coxph.fit)
> my.survfit.objek
Call: survfit(formula = coxph.fit)
records n.max n.start events median 0.95LCL 0.95UCL
154 154 154 99 14 12 18
> my.survfit.objek
Call: survfit(formula = coxph.fit)
records n.max n.start events median 0.95LCL 0.95UCL
154 154 154 99 14 12 18
Pengembalian obejek dengan "survfit ()" memiliki karakteristik khas dan sifat seperti sebelumnya. Misalnya, fungsi survival dasar dapat diplot menggunakan plot () fungsi. Di sini kita akan menjalankan satu tes lokal. Di atas, $Z11$ adalah variabel faktor direpresentasikan sebagai koefisien kedua, ketiga, dan keempat dari koefisien coxph.fit. Saat kita membangun $C$ dan $d,$ maka kita menghitung uji statistik $X_W^2$. Dalam kode di bawah, diperlihatkan obyek $t_1$ dan $t_2$ dihitung sebagai cara untuk mendapatkan uji statistik.
> C <- matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),nrow=3,byrow=TRUE)
> d <- rep(0, 3)
> t1 <- C %*% co - d
> t2 <- C %*% va %*% t(C)
> XW2 <- c(t(t1) %*% solve(t2) %*% t1)
> pchisq(XW2, 3, lower.tail=FALSE)
[1] 0.1025018
> detach(burn)
Rangkaian fungsi c () digunakan dalam perhitungan XW2 membuat XW2 vektor dari pada matriks $1\times 1$. Nilai-p dari uji ini adalah 0,103, yang berarti Z11 tidak uji statistik yang signifikan untuk tingkat signifikan dari $\alpha= 0.05.$
Model resiko proportional Cox kovariat tergantung waktu
Menggunakan kovariat tergantung waktu dalam R adalah latihan dalam mengorganisasi. Sebelumnya dianggap kovariat adalah atribut yang tidak berubah, seperti perlakuan terhadap kelompok obsevasi atau pengontrolan kelompok atau ras pasien. Saat Ini kita mempertimbangkan tergantung waktu kovariat, seperti intervensi atau faktor lingkungan yang mungkin mengakibatkan perubahan pertengahan studi.
Untuk menggunakan kovariat tergantung waktu di R, kita akan menerapkan pemotongan periode disebelah kiri (melakukan obsevasi sebelum periode selesai) secara bebas. Sebagai contoh, jika ada intervensi untuk pasien $i,$ maka kita membagi pasien $i$ menjadi dua pengamatan terpisah: pra dan pasca-intervensi. Lebih eksplisit, misalkan intervensi pasien terjadi pada waktu $t_i $= 45 dan acara untuk pasien saya diamati pada $t_{kejadian}$= 58. Kita akan membagi catatan pasien ini menjadi dua potongan : 0-45 dan 45-58. kovariat intervensi dapat diberikan nilai berbeda untuk setiap interval. Menerapkan metode memecah waktu di R membutuhkan sedikit ketelitian. Memulai dan mengakhiri untuk masing-masing interval yang dibangun,
dan data tersensor dikanan harus dilacak untuk setiap pasangan dari awal hingga akhir.
Asumsikan data berikut adalah contoh simulasi. Catatan pasien dari klinik penyalahgunaan alkohol diperoleh selama 150 hari setelah pasien sudah dinyatakan sembuh. Hal yang menjadi perhatian adalah kadar alkohol pasien dan empat variabel yang tersedia sebagai berikut:
> data(relapse) #memanggil data base
> relapse #menampilkan data
event delta gender inter
1 150 FALSE 0 84
2 53 TRUE 1 50
3 12 TRUE 1 NA
4 150 FALSE 0 89
5 150 FALSE 1 77
6 135 TRUE 1 7
Pada output, (Intercept) dan Log (skala) sesuai dengan perkiraan $\mu$ dan $log\sigma$. Estimasi yang lainnya akan sesuai dengan koefisien kovariat. Kita juga mungkin mempertimbangkan model eksponensial:
> srFitExp <- survreg(Surv(time, delta) ~ as.factor(stage) + age, dist="exponential")
> summary(srFitExp)
Call:
survreg(formula = Surv(time, delta) ~ as.factor(stage) + age,
dist = "exponential")
Value Std. Error z p
(Intercept) 3.7550 0.9902 3.792 1.49e-04
as.factor(stage)2 -0.1456 0.4602 -0.316 7.52e-01
as.factor(stage)3 -0.6483 0.3552 -1.825 6.80e-02
as.factor(stage)4 -1.6350 0.3985 -4.103 4.08e-05
age -0.0197 0.0142 -1.388 1.65e-01
Scale fixed at 1
Exponential distribution
Loglik(model)= -141.9 Loglik(intercept only)= -151.1
Chisq= 18.44 on 4 degrees of freedom, p= 0.001
Number of Newton-Raphson Iterations: 4
n= 90
> C <- matrix(c(0,1,0,0,0,0,1,0,0,0,0,1),nrow=3,byrow=TRUE)
> d <- rep(0, 3)
> t1 <- C %*% co - d
> t2 <- C %*% va %*% t(C)
> XW2 <- c(t(t1) %*% solve(t2) %*% t1)
> pchisq(XW2, 3, lower.tail=FALSE)
[1] 0.1025018
> detach(burn)
Rangkaian fungsi c () digunakan dalam perhitungan XW2 membuat XW2 vektor dari pada matriks $1\times 1$. Nilai-p dari uji ini adalah 0,103, yang berarti Z11 tidak uji statistik yang signifikan untuk tingkat signifikan dari $\alpha= 0.05.$
Model resiko proportional Cox kovariat tergantung waktu
Menggunakan kovariat tergantung waktu dalam R adalah latihan dalam mengorganisasi. Sebelumnya dianggap kovariat adalah atribut yang tidak berubah, seperti perlakuan terhadap kelompok obsevasi atau pengontrolan kelompok atau ras pasien. Saat Ini kita mempertimbangkan tergantung waktu kovariat, seperti intervensi atau faktor lingkungan yang mungkin mengakibatkan perubahan pertengahan studi.
Untuk menggunakan kovariat tergantung waktu di R, kita akan menerapkan pemotongan periode disebelah kiri (melakukan obsevasi sebelum periode selesai) secara bebas. Sebagai contoh, jika ada intervensi untuk pasien $i,$ maka kita membagi pasien $i$ menjadi dua pengamatan terpisah: pra dan pasca-intervensi. Lebih eksplisit, misalkan intervensi pasien terjadi pada waktu $t_i $= 45 dan acara untuk pasien saya diamati pada $t_{kejadian}$= 58. Kita akan membagi catatan pasien ini menjadi dua potongan : 0-45 dan 45-58. kovariat intervensi dapat diberikan nilai berbeda untuk setiap interval. Menerapkan metode memecah waktu di R membutuhkan sedikit ketelitian. Memulai dan mengakhiri untuk masing-masing interval yang dibangun,
dan data tersensor dikanan harus dilacak untuk setiap pasangan dari awal hingga akhir.
Asumsikan data berikut adalah contoh simulasi. Catatan pasien dari klinik penyalahgunaan alkohol diperoleh selama 150 hari setelah pasien sudah dinyatakan sembuh. Hal yang menjadi perhatian adalah kadar alkohol pasien dan empat variabel yang tersedia sebagai berikut:
- kejadian (event) = menjelaskan diamati atau disensor kekambuhan,
- delta = menjelaskan apakah kejadian diamati (TRUE) atau disensor (FALSE),
- gender = kovariat waktu-independen, dan
- int = menyatakan kovariat tergantung waktu pasien memiliki intervensi, di mana beberapa pasien kambuh sebelum intervensi mereka dijadwalkan dan dinyatakan dengan NA.
> data(relapse) #memanggil data base
> relapse #menampilkan data
event delta gender inter
1 150 FALSE 0 84
2 53 TRUE 1 50
3 12 TRUE 1 NA
4 150 FALSE 0 89
5 150 FALSE 1 77
6 135 TRUE 1 7
....... dan seterusnya
Data ini dapat dimodelkan dengan menggunakan dua langkah:
- Buatlah catatan hidup yang mungkin termasuk kiri dan kanan pengamatan disensor. kelangsungan hidup catatan setiap pasien dengan intervensi ini dibagi menjadi dua survival: satu sebelum intervensi dan satu setelah. Langkah ini sebagian besar sudah ada dalam pembukuan dan pemrograman R.
- Catatan hidup baru dapat dijalankan melalui coxph () fungsi.
Pertama-tama menginisialisasi vektor baru untuk mewakili variabel catatan hidup. Variabel $t_1$ dan $t_2$ merupakan awal dan akhir, masing-masing $d$ merupakan apakah kambuh yang diamati (TRUE) atau benar-disensor (FALSE), $g$ merupakan jenis kelamin, dan saya mewakili apakah pasien menjalani perawatan intervensi (kita harus berhati-hati untuk tidak menimpa variabel ini dalam satu lingkaran).
> N<-dim(relapse)[1] # menanyakan R berapa banyak data
> t1<-rep(0,N+sum(!is.na(relapse$int))) #inisial waktu dimulai 0
> t2<-rep(-1,length(t1)) # membangun vektor untuk akhir waktu
> d<-rep(-1,length(t1)) # apakah data disensor
> g<-rep(-1,length(t1)) # jenis kelamin kovariate
> i<-rep(FALSE,length(t1)) # inisialisasi intervensi pada FALSE
Selanjutnya, setiap pasien diperiksa apakah mereka memiliki intervensi. Jika mereka tidak, maka rekor mereka hanya disalin ke variabel baru yang akan digunakan. Jika ada intervensi, observasi dibagi menjadi dua bagian: periode sebelum intervensi dan periode setelah intervensi. Periode berikut adalah setelah intervensi yang tersisa-dipotong untuk waktu intervensi.
> j<-1
> for(ii in 1:dim(relapse)[1]){if(is.na(relapse$int[ii]))
{ # bukan intervensi, copy survival rekord
t2[j]<-relapse$event[ii]
d[j]<-relapse$delta[ii]
g[j]<-relapse$gender[ii]
j<- j+1
} else { # intervensi, catatan perpecahan
g[j+0:1]<-relapse$gender[ii] # jenis kelamin sama setiap waktu
d[j] <- 0 # bukan relapse pra-intervensi
d[j+1] <- relapse$delta[ii] # relapse terjadi pasca-intervensi?
i[j+1] <- TRUE # intervensi kovariat, pasca-intervensi
t2[j] <- relapse$int[ii]-1 # akhir pra-intervensi
t1[j+1] <- relapse$int[ii]-1 # awal pasca-intervensi
t2[j+1] <- relapse$event[ii] # akhir of pasca-intervensi
j <- j+2 # penambahan dua rekord
}
}
Arti dari variabel adalah sebagai berikut:
- $t_1$ merupakan awal setiap interval,
- $t_2$ akhir setiap interval,
- $e$ menunjukkan jika acara diamati,
- $g$ adalah kovariat gender, dan
- $i$ adalah indikator variabel untuk apakah interval ini dikaitkan dengan pasca-intervensi pasien.
Sementara banyak pasien memiliki waktu bervariasi kovariat (yaitu, intervensi), masing-masing interval baru memiliki kovariat yang tidak berubah, yang memungkinkan untuk objek survival dan Model PH Cox yang akan dibangun:
> mySurv<-Surv(t1,t2,d)
> myCPH<-coxph(mySurv~g+i)
Contoh ini adalah kasus sederhana; tergantung pada kovariat waktu tunggal, dan hal itu paling banyak berubah setiap kasus. Dalam beberapa kasus, mungkin ada banyak kovariat waktu bervariasi, bahkan ada juga yang berubah setiap satuan waktu. Dalam sebagian besar kasus ini, metode yang sama dapat digunakan. Setiap kali perubahan kovariat dari satu satuan waktu ke depan, membagi interval menjadi dua dan menggunakan pemotongan waktu observasi data tersensor dikanan seseuai yang diperlukan. Untuk lebih jelas lagi silahkan anda membaca kembali pada bagian "Target Survival" diatas.
Model kesalahan percepatan waktu
survreg(formula, dist='weibull')
Suatu model kesalahan percepatan waktu atau Accelerated failure-time (AFT) models adalah model parametrik dengan kovariat dan kesalahan waktu. Berikut fungsi kelangsungan hidup bentuk $S (x|Z) = S_0 (x*exp [\theta' Z]),$ dimana $S_0$ adalah fungsi untuk tingkat survival awal. Istilah $exp[\theta'Z]$ disebut sebagai faktor akselerasi. Model AFT menggunakan kovariat untuk menempatkan orang pada skala waktu berbeda. Perhatikan skala dengan kovariat $S(t|Z)$ melalui $exp [\theta'Z]$. Model AFT dapat ditulis dalam bentuk log-linear, dimana log waktu kegagalan (disebut $logX$) berhubungan linier dengan rata-rata $\mu$, faktor akselerasi, dan standar kesalahan dengan istilah $\sigma W$:
$$log X= \mu-\theta'Z+\sigma W$$
Fungsi "survreg ( )" dalam paket survival digunakan untuk pemodelan AFT. Argumen pertama adalah rumus, di mana objek kemunduran survival pada prediksi. Argumen "dist" memiliki beberapa pilihan untuk menggambarkan model parametrik digunakan ("weibull", "eksponensial", "gaussian", "logistik", "lognormal", atau "loglogistic").
> data(larynx)
> attach(larynx)
> srFit<-survreg(Surv(time,delta)~as.factor(stage)+age, dist="weibull")
> summary(srFit)
Call:
survreg(formula = Surv(time, delta) ~ as.factor(stage) + age,
dist = "weibull")
Value Std. Error z p
(Intercept) 3.5288 0.9041 3.903 9.50e-05
as.factor(stage)2 -0.1477 0.4076 -0.362 7.17e-01
as.factor(stage)3 -0.5866 0.3199 -1.833 6.68e-02
as.factor(stage)4 -1.5441 0.3633 -4.251 2.13e-05
age -0.0175 0.0128 -1.367 1.72e-01
Log(scale) -0.1223 0.1225 -0.999 3.18e-01
Scale= 0.885
Weibull distribution
Loglik(model)= -141.4 Loglik(intercept only)= -151.1
Chisq= 19.37 on 4 degrees of freedom, p= 0.00066
Number of Newton-Raphson Iterations: 5
n= 90
Pada output, (Intercept) dan Log (skala) sesuai dengan perkiraan $\mu$ dan $log\sigma$. Estimasi yang lainnya akan sesuai dengan koefisien kovariat. Kita juga mungkin mempertimbangkan model eksponensial:
> srFitExp <- survreg(Surv(time, delta) ~ as.factor(stage) + age, dist="exponential")
> summary(srFitExp)
Call:
survreg(formula = Surv(time, delta) ~ as.factor(stage) + age,
dist = "exponential")
Value Std. Error z p
(Intercept) 3.7550 0.9902 3.792 1.49e-04
as.factor(stage)2 -0.1456 0.4602 -0.316 7.52e-01
as.factor(stage)3 -0.6483 0.3552 -1.825 6.80e-02
as.factor(stage)4 -1.6350 0.3985 -4.103 4.08e-05
age -0.0197 0.0142 -1.388 1.65e-01
Scale fixed at 1
Exponential distribution
Loglik(model)= -141.9 Loglik(intercept only)= -151.1
Chisq= 18.44 on 4 degrees of freedom, p= 0.001
Number of Newton-Raphson Iterations: 4
n= 90
Ketika $\sigma$ = 1, model Weibull setara dengan model eksponensial. Anda dapat mempertimbangkan dua strategi ini untuk memilih model terbaik:
- Uji rasio kemungkinan, yang mengevaluasi nol hipotesis $\sigma$= 1 melawan dua sisi alternatif, dan
- Pemeriksaan signifikansi dari koefisien Log (skala) (lihat output ringkasan (srFit)).
Dalam contoh ini, kedua pendekatan menghasilkan kesimpulan yang sama: tidak ada alasan yang memadai untuk menolak hipotesis $\sigma$= 1 $(H_0).$ Untuk alasan ini, akan saya perlihatkan model eksponesial sederhan. Jika tertarik anda dapat menjelajahi banyak komponen yang disimpan dalam "survreg ( )" berikut:
> srFitExp$coeff # koefisien kovariat
> srFitExp$icoef # koefisien intercept and skala
> srFitExp$var # matriks varians-kovarians
> srFitExp$loglik # log-likelihood
Akhirnya, sebagai programer yang baik jangan pernah lupa untuk melepaskan kumpulan data.
> detach(larynx)
Referensi
- Terry Therneau and original Splus\R port by Thomas Lumley (2011). survival: Survival analysis, including penalised likelihood. R package version 2.36-10. http://CRAN.R-project.org/package=survival.
- R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
- David M Diez (2012). OIsurv: Supplement to survival analysis tutorial. R package version 1.0. http://CRAN.R-project.org/package=OIsurv
- Original by Klein, Moeschberger and modi cations by Jun Yan (2010). KMsurv: Data sets from Klein and Moeschberger (1997), Survival Analysis. R package version 0.1-4. http://CRAN.R-project.org/package=KMsurv.
- Klein, John P., and Melvin L. Moeschberger. Survival Analysis: Techniques for Censored and Truncated Data. New York: Springer, 2003.
- Fox, John 2002. Cox Proportional-Hazards Regression for Survival Data. Appendix to An R and S-PLUS Companion to Applied Regression.
