Tabel Kehidupan

Tabel kehidupan (Life Table) merupakan salah satu pendekatan non-parametrik untuk penduga fungsi daya tahan, selain Kaplan Meier. Tabel kehidupan akan sangat berguna jika kita tidak memiliki data setiap individu, melainkan hanya berupa ringkasan atau tabulasi data. Untuk membentuk tabel kehidupan pengamatan dibagi menjadi beberapa interval waktu, selanjutnya fungsi daya tahan diduga untuk setiap interval waktu. Penduga fungsi daya tahan merupakan perkalian proporsi individu yang bertahan dan total individu beresiko pada tiap interval.

Dalam menyusun tabel kehidupan diasumsikan kejadian terjadi di akhir interval, dan pengamatan yang tersensor terjadi secara merata sepanjang interval, sehingga biasanya dilakukan penyesuaian (adjusment) banyaknya individu beresiko.

Penduga fungsi daya tahan dapat dihitung sebagai berikut:

\(\hat S(t+1)=S(t)\times p(t+1)\)

dengan,

\(p(t)=1-\frac{d_t}{n_t^*}\)

\(n^*_t=n_t-\frac{c_t}{2}\)

\(d_t\) merupakan banyaknya individu yang mengalammi kejadian pada interal waktu \(t\)

\(n_t\) merupakan banyaknya individu beresiko pada interval waktu \(t\)

\(c_t\) merupakan banyaknya individu yang tersensor pada interal waktu \(t\)

Ilustrasi Tabel Kehidupan

Pada ilustrasi ini akan kita gunakan data pasien kanker prostat. Data merupakan hasil simulasi dari data daya tahan dengan kasus competing risk oleh Lu-Yao et al. (2009), yaitu kematian akibat kanker prostat dan kematian yang lain. Pada ilustrasi ini, Kejadian yang kita amati adalah kematian akibat kanker prostat, sementara kematian yang disebabkan hal lain dianggap sebagai data tersensor. Data memuat peubah sebagai berikut:

  • grade: suatu peubah faktor dengan level mode (moderately differentiated) dan poor (poorly differentiated)

  • stage: jenis kanker dengan level T1ab (Stage T1, clinically diagnosed), T1c (Stage T1, diagnosed via a PSA test), T2 (Stage T2)

  • ageGroup: kelompok umur dengan level 66-69, 70-74, 75-79 80+

  • survTime: waktu dari diagnosis sampai meninggal atau terakhir kali diketahui hidup

  • status: indikator kejadian dengan 0 (tersensor), 1 (meninggal karna kanker prostat), 2 (meninggal karena sebab lain)

# required packages
# install.packages('asaur')
# install.packages('biostat3')
# install.packages('bshazard')
library(asaur)
library(biostat3)
library(bshazard)
library(tidyverse)
data("prostateSurvival")
str(prostateSurvival)
## 'data.frame':    14294 obs. of  5 variables:
##  $ grade   : Factor w/ 2 levels "mode","poor": 1 1 2 1 1 2 2 1 1 2 ...
##  $ stage   : Factor w/ 3 levels "T1ab","T1c","T2": 2 1 2 3 2 3 2 1 2 3 ...
##  $ ageGroup: Factor w/ 4 levels "66-69","70-74",..: 4 3 3 2 2 3 4 4 3 3 ...
##  $ survTime: int  18 23 37 27 42 38 18 78 47 13 ...
##  $ status  : int  0 0 0 0 0 2 0 0 0 0 ...

Sebelumnya kita harus mengubah peubah indikator kejadian menjadi bentuk dikotomus, yaitu 0 untuk tersensor dan 1 untuk kejadian.

prostat=prostateSurvival%>%
  mutate(status=recode(status,'0'=0, '1'=1, '2'=0))

head(prostat)

Untuk menentukan banyaknya interval waktu yang akan digunakan dalam tabel kehidupan, perlu kita ketahui terlebih dahulu nilai-nilai waktu bertahan yang ada dalam data

sort(unique(prostat$survTime))
##   [1]   0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17
##  [19]  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35
##  [37]  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53
##  [55]  54  55  56  57  58  59  60  61  62  63  64  65  66  67  68  69  70  71
##  [73]  72  73  74  75  76  77  78  79  80  81  82  83  84  85  86  87  88  89
##  [91]  90  91  92  93  94  95  96  97  98  99 100 101 102 103 104 105 106 107
## [109] 108 109 110 111 112 113 114 115 116 117 118 119

Kita akan bagi data ke dalam 12 rentang waktu.

breaks=seq(0,120,10)
lifetabel=lifetab2(Surv(survTime,status)~1,
                   data=prostat, breaks = breaks)
print(lifetabel)
##         tstart tstop nsubs nlost   nrisk nevent      surv          pdf
## 0-10         0    10 14294  2481 13053.5     95 1.0000000 0.0007277742
## 10-20       10    20 11718  2466 10485.0    134 0.9927223 0.0012687151
## 20-30       20    30  9118  1767  8234.5    114 0.9800351 0.0013567794
## 30-40       30    40  7237  1206  6634.0    110 0.9664673 0.0016025234
## 40-50       40    50  5921  1072  5385.0     90 0.9504421 0.0015884826
## 50-60       50    60  4759  1009  4254.5     83 0.9345573 0.0018232049
## 60-70       60    70  3667   880  3227.0     68 0.9163252 0.0019308991
## 70-80       70    80  2719   736  2351.0     44 0.8970162 0.0016788053
## 80-90       80    90  1939   592  1643.0     24 0.8802282 0.0012857867
## 90-100      90   100  1323   517  1064.5     25 0.8673703 0.0020370369
## 100-110    100   110   781   449   556.5      9 0.8469999 0.0013698112
## 110-120    110   120   323   320   163.0      3 0.8333018 0.0015336843
## 120-Inf    120   Inf     0     0     0.0      0 0.8179650           NA
##               hazard      se.surv       se.pdf    se.hazard
## 0-10    0.0007304321 0.0000000000 7.439585e-05 7.494025e-05
## 10-20   0.0012862354 0.0007439585 1.089018e-04 1.111115e-04
## 20-30   0.0013940691 0.0013135023 1.262045e-04 1.305634e-04
## 30-40   0.0016719866 0.0018083900 1.515522e-04 1.594121e-04
## 40-50   0.0016853933 0.0023363714 1.660816e-04 1.776497e-04
## 50-60   0.0019700926 0.0028345156 1.982382e-04 2.162353e-04
## 60-70   0.0021296586 0.0034133314 2.317873e-04 2.582444e-04
## 70-80   0.0018892228 0.0040659992 2.508254e-04 2.847984e-04
## 80-90   0.0014714899 0.0047122039 2.606270e-04 3.003585e-04
## 90-100  0.0023764259 0.0053243590 4.027891e-04 4.752516e-04
## 100-110 0.0016304348 0.0065758001 4.530213e-04 5.434602e-04
## 110-120 0.0018575851 0.0078971731 8.774071e-04 1.072431e-03
## 120-Inf           NA 0.0117070066           NA           NA

Dari hasil di atas, pada rentang waktu pertama, 0-10, terdapat sebanyak 95 individu yang meninggal akibat kanker prostat dan 2481 individu tersensor (lost kontak atau meninggal karena hal lain). Total individu beresiko pada rentang waktu pertama adalah 14294, karena terdapat data tersensor sehingga dilakukan penyesuaian jumlah individu beresiko menjadi \(14294-\frac{2481}{2}=13053.5\). Peluang daya tahan untuk rentang waktu pertama selanjutnya dapat kita hitung sebagai berikut:

\(p(1)=1-\frac{d_1}{n_1^*}=1-\frac{95}{13053.5}=0.9927223\)

\(\hat S(1)=S(0)\times p(1)=1\times 0.9927221=0.9927223\)

Dengan cara yang sama, peluang daya tahan untuk rentang waktu kedua adalah:

\(n^*_2=11718-\frac{2466}{2}=10485\)

\(p(2)=1-\frac{d_2}{n_2^*}=1-\frac{134}{10485}=0.9872198\)

\(\hat S(2)=S(1)\times p(2)=0.9927223\times 0.9872198=0.9800351\)

Untuk rentang waktu selanjutnya dihitung dengan cara yang sama. Perhatikan hasil tabel kehidupan di atas. Pada baris pertama kolom surv bernilai 1. Nilai tersebut merupakan nilai \(S(0)\). Sedangkan nilai \(S(1)\) berada pada baris kedua, dimana nilainya 0.9927223 sesuai dengan perhitungan manual yang dilakukan di atas.

NB.: Untuk menyusun tabel kehidupan dengan melibatkan peubah kategorik tidak dapat dilakukan secara langsung dengan mengganti \(1\) di kanan notasi ~ dengan peubah kategorik. Tetapi harus dilakukan secara manual dengan memisahkan data berdasarkan kategori.

Fungsi Bahaya

Fungsi bahaya memberikan informasi yang berlawanan dengan fungsi daya tahan. Jika dari fungsi daya tahan (Kaplan Meier atau Tabel Kehidupan) kita memperoleh informasi terkait peluang individu untuk bertahan, dari fungsi bahaya kita memperoleh informasi laju kegagalan untuk bertahan. Fungsi bahaya dapat juga kita definisikan sebagai potensi individu untuk mengalami kejadian pada waktu tertentu jika diketahui individu tersebut telah bertahan sampai waktu t, sehingga fungsi bahaya menerapkan konsep peluang bersyarat.

\(h(t)=\lim_{\Delta t\to 0}\frac{P(t\le T\le t+\Delta t|T\ge t)}{\Delta t}\)

atau dapat dinyatakan sebagai:

\(h(t)=-\frac{d}{dt}ln\;S(t)\)

Ilustrasi Fungsi Bahaya

Untuk ilustrasi fungsi bahaya kita gunakan data yang sama dengan ilustrasi tabel kehidupan di atas.

hazard=bshazard(Surv(survTime,status)~1, data=prostat,
                verbose=F)
plot(hazard, xlim=c(0,120),ylim=c(0.0005,0.0035))

rbind(head(summary(hazard)$HazardEstimates,10),
      tail(summary(hazard)$HazardEstimates,10))
##        time       hazard        lower        upper
##           1 0.0007657078 0.0006206511 0.0009446667
##           2 0.0007908849 0.0006501300 0.0009621137
##           3 0.0008168899 0.0006803379 0.0009808496
##           4 0.0008437500 0.0007110871 0.0010011630
##           5 0.0008714053 0.0007420848 0.0010232620
##           6 0.0008986961 0.0007727644 0.0010451498
##           7 0.0009268415 0.0008037032 0.0010688463
##           8 0.0009558684 0.0008346371 0.0010947087
##           9 0.0009854458 0.0008651789 0.0011224307
##          10 0.0010135387 0.0008948861 0.0011479234
## [110,]  110 0.0017792733 0.0011334166 0.0027931598
## [111,]  111 0.0017749808 0.0011089857 0.0028409358
## [112,]  112 0.0017709497 0.0010845916 0.0028916533
## [113,]  113 0.0017669662 0.0010599784 0.0029455032
## [114,]  114 0.0017629917 0.0010352374 0.0030023449
## [115,]  115 0.0017590262 0.0010104816 0.0030620777
## [116,]  116 0.0017551721 0.0009857883 0.0031250414
## [117,]  117 0.0017513338 0.0009610708 0.0031914090
## [118,]  118 0.0017475039 0.0009364251 0.0032610936
## [119,]  119 0.0017436824 0.0009119364 0.0033340354

Hasil di atas menampilkan nilai fungsi bahaya untuk 10 waktu awal dan 10 waktu akhir. Pada saat \(t=1\), potensi seorang individu untuk mengalami kejadian (dalam kasus ini meninggal akibat kanker prostat) jika individu tersebut diketahui bertahan hingga waktu $t=1 $adalah 0.000766, begitupun untuk \(t=110\), nilai 0.001779 dapat diinterpretasikan sebagai potensi seorang individu untuk meninggal akibat kanker prostat jika individu tersebut diketahui bertahan hingga waktu $t=110 $adalah 0.001779. Jika kita perhatikan grafik laju bahaya di atas, potensi individu untuk meninggal akibat kanker prostat tertinggi adalah ketika \(t=64\) dengan laju bahaya 0.001988.

summary(hazard)$HazardEstimates[
  summary(hazard)$HazardEstimates[,2]==max(summary(
    hazard)$HazardEstimates[,2]),]
##         time       hazard        lower        upper 
## 64.000000000  0.001988382  0.001732070  0.002282622

Selanjutnya kita akan bandingkan laju bahaya untuk tiga kategori stage kanker prostat.

hazard2=bshazard(Surv(survTime,status)~stage, data=prostat,
                verbose=F)
## NOTE: entry.status has been set to 0 for all.
## NOTE: Dropping  229  rows with duration of follow up < tol
plot(hazard2,
     conf.int = F, overall = F,
     lty = 1:3, lwd=2, ylim = c(0.0005,0.0035))
legend(0, 0.0035,
       c("T1ab", "T1c", "T2")[1:3],
       lty = 1:3, title = "Stage")

Tampak bahwa individu dengan kategori kanker prostat T2 memiliki laju bahaya paling tinggi dibanding dua kategori yang lain. Sementara ketgori T1c memiliki laju bahaya yang paling rendah dibanding dua kategori lain. Akan tetapi jika kita perhatikan nilai laju bahaya ketiganya sangat kecil, dan mungkin saja perbedaan laju bahaya untuk ketiga kategori tidak signifikan.