1 Pemahaman Masalah: Predictive Maintenance

1.1 Konteks Bisnis

Predictive Maintenance (PdM) adalah strategi pemeliharaan berbasis data yang menggunakan sensor IoT dan riwayat operasional untuk memprediksi kapan sebuah mesin akan mengalami kegagalan. Berbeda dengan:

  • Reactive Maintenance: memperbaiki setelah mesin rusak — biaya downtime dan kerusakan sekunder sangat tinggi
  • Preventive Maintenance: memperbaiki berdasarkan jadwal tetap — sering tidak efisien (terlalu dini atau terlambat)
  • Predictive Maintenance: memperbaiki tepat sebelum mesin gagal berdasarkan sinyal data aktual

Dalam konteks industri IoT modern, PdM memanfaatkan sensor tekanan, kelembaban, dan suhu yang terekam secara real-time, dikombinasikan dengan metadata operasional (tim pemeliharaan, vendor komponen) untuk memperkirakan Remaining Useful Life (RUL) setiap mesin secara individual.

1.2 Konsep Dasar Survival Analysis

Survival Analysis adalah kerangka statistik yang dirancang khusus untuk memodelkan time-to-event data — berapa lama hingga suatu kejadian (kegagalan) terjadi. Kerangka ini memiliki tiga konsep fundamental:

Failure Event

Event adalah kejadian yang menjadi fokus studi — dalam konteks PdM, event adalah kegagalan mesin. Pada dataset ini, event ditandai dengan broken = 1. Event harus terdefinisi jelas: kapan tepatnya mesin dinyatakan “gagal” dan apa kriterianya.

Censored Data

Censoring terjadi ketika waktu kegagalan tidak terobservasi penuh. Pada dataset ini, mesin dengan broken = 0 adalah right-censored — mesin masih beroperasi saat periode pengamatan berakhir. Kita hanya tahu mesin bertahan setidaknya selama lifetime minggu. Mengabaikan data censored akan menghasilkan estimasi bias ke arah underestimasi waktu survival.

Survival Function & Hazard Function

\[S(t) = P(T > t) = \text{probabilitas mesin masih berfungsi pada minggu ke-}t\]

\[h(t) = \lim_{\Delta t \to 0}\frac{P(t \leq T < t+\Delta t \mid T \geq t)}{\Delta t} = \text{laju kegagalan sesaat}\]

Hubungan keduanya: \(S(t) = \exp\!\left(-\int_0^t h(u)\,du\right)\)

1.3 Mengapa Linear MTLR untuk PdM?

Linear MTLR (Multi-Task Logistic Regression for Survival) adalah model survival modern yang secara khusus unggul untuk PdM karena:

  1. Tidak mengasumsikan proportional hazards — berbeda dari Cox PH yang mengasumsikan efek prediktor konstan sepanjang waktu
  2. Output probabilistik lengkap — menghasilkan distribusi \(P(T = t_j | \mathbf{x})\) untuk setiap time bin, bukan sekadar hazard ratio
  3. Temporal smoothing — regularisasi antar task yang berdekatan waktu menghasilkan kurva survival yang halus dan realistis
  4. Interpretasi langsung — probabilitas kegagalan per minggu dapat langsung digunakan untuk scheduling inspeksi

2 Set Up

Sebuah perusahaan manufaktur mengoperasikan banyak mesin dalam lini produksinya. Pabrik telah memasang sensor IoT pada setiap mesin untuk merekam data kondisi operasional secara real-time — tekanan, kelembapan, dan suhu.

Setiap kali mesin rusak, seluruh lini produksi terhenti, menimbulkan kerugian ribuan dolar dari biaya perbaikan dan denda keterlambatan pengiriman. Manajer pabrik meminta tim Data Science untuk menemukan cara yang lebih proaktif dalam mengelola pemeliharaan.


3 Dataset & Pemahaman Data

3.1 Deskripsi Dataset

## ── Ringkasan Dataset ──────────────────────────
## Total mesin       : 1000
## Kegagalan (broken=1): 397 ( 39.7 % event rate)
## Censored (broken=0): 603 ( 60.3 %)
## Rentang lifetime  : 1 – 93 weeks
## Min lifetime gagal: 60 weeks
## Max lifetime gagal: 93 weeks

Dataset berisi 1.000 mesin industri yang dimonitor selama maksimal 93 minggu. Satuan waktu adalah minggu (weeks). Event rate 39,7% berarti 397 dari 1.000 mesin mengalami kegagalan dalam periode observasi, sementara 603 mesin dicensored — pengamatan berakhir sebelum kegagalan terjadi. Ini adalah proporsi censoring yang tinggi dan wajar untuk data industrial.

3.2 Kamus Variabel

Variabel Tipe Satuan Keterangan
lifetime Integer Weeks Waktu observasi: minggu ke berapa mesin gagal (broken=1) atau dicensored (broken=0)
broken Binary 1 = event kegagalan terobservasi; 0 = right-censored
pressureInd Numeric Indeks (~100) Pembacaan sensor tekanan IoT — snapshot tunggal
moistureInd Numeric Indeks (~100) Pembacaan sensor kelembaban IoT — snapshot tunggal
temperatureInd Numeric Indeks (~100) Pembacaan sensor suhu IoT — snapshot tunggal
team Categorical Tim pemeliharaan yang bertanggung jawab: TeamA, TeamB, TeamC
provider Categorical Vendor/pemasok komponen mesin: Provider1, Provider2, Provider3, Provider4

4 Preprocessing & Feature Engineering

4.1 Pemeriksaan Data

## ── Missing Values ─────────────────────────────
## # A tibble: 7 × 3
##   Variabel       Missing Persentase
##   <chr>            <int> <chr>     
## 1 lifetime             0 0%        
## 2 broken               0 0%        
## 3 pressureInd          0 0%        
## 4 moistureInd          0 0%        
## 5 temperatureInd       0 0%        
## 6 team                 0 0%        
## 7 provider             0 0%
## 
## Duplicate rows: 0
## 
## Tipe data setiap kolom:
##       lifetime         broken    pressureInd    moistureInd temperatureInd 
##      "integer"      "integer"      "numeric"      "numeric"      "numeric" 
##           team       provider 
##       "factor"       "factor"

Hasil: Tidak ada missing value dan tidak ada duplikat. Dataset sudah bersih dan siap digunakan langsung. Pemeriksaan tipe data mengonfirmasi lifetime dan broken sudah integer, ketiga sensor sudah numeric, dan team/provider sudah berupa factor.

4.2 Encoding Variabel Kategorik

Variabel kategorik (team, provider) perlu dikonversi ke representasi numerik untuk model MTLR dan pembanding. Encoding dilakukan berdasarkan urutan event rate (dari terendah ke tertinggi) agar angka encoding membawa makna ordinal yang relevan: angka lebih besar = risiko kegagalan lebih tinggi.

## ── Event Rate per Provider (urutan encoding) ──
## # A tibble: 4 × 4
##   provider  event_rate     n encoding
##   <fct>          <dbl> <int>    <int>
## 1 Provider4       31.9   238        1
## 2 Provider2       34.2   266        2
## 3 Provider1       45.7   254        3
## 4 Provider3       47.1   242        4
## 
## ── Event Rate per Team ────────────────────────
## # A tibble: 3 × 4
##   team  event_rate     n encoding
##   <fct>      <dbl> <int>    <int>
## 1 TeamA       36.6   336        1
## 2 TeamC       40.3   308        2
## 3 TeamB       42.1   356        3
## 
## ── Contoh Data setelah Encoding (5 baris) ────
##   lifetime broken  team team_enc  provider provider_enc
## 1       56      0 TeamA        1 Provider4            1
## 2       81      1 TeamC        2 Provider4            1
## 3       60      0 TeamA        1 Provider1            3
## 4       86      1 TeamC        2 Provider2            2
## 5       34      0 TeamB        3 Provider1            3

Alasan encoding berbasis event rate: Provider dengan angka encoding lebih besar memiliki event rate lebih tinggi. Ini membuat korelasi antara nilai numerik encoding dengan variabel outcome (broken) lebih bermakna dibanding encoding alfabet biasa. Untuk MTLR dan model linear, representasi ordinal yang konsisten dengan fenomena yang diukur menghasilkan koefisien yang lebih mudah diinterpretasikan.

4.3 Pembentukan Survival Object

## ── Survival Object ────────────────────────────
##  [1] 56+ 81  60+ 86  34+ 30+ 68+ 65  23+ 81
## 
## Format: waktu+ berarti censored (broken=0), waktu saja berarti event (broken=1)
## Contoh: '56' = censored pada minggu ke-56
##         '81' dengan tanda event = gagal pada minggu ke-81

Survival object Surv(lifetime, broken) menggabungkan informasi waktu dan status event dalam satu struktur. Tanda + pada output menandakan observasi censored. Objek ini menjadi input standar untuk semua model survival (MTLR, Cox, RSF) sehingga seluruh model menangani censoring dengan cara yang sama dan konsisten.


5 Exploratory Data Analysis

5.1 Distribusi Variabel Numerik

lifetime (weeks): Distribusi bimodal dengan puncak di minggu 20–40 (mesin censored awal) dan minggu 70–90 (mesin yang bertahan masuk zona kegagalan). Bimodal ini mencerminkan dua “populasi” mesin — yang diobservasi singkat dan yang bertahan lama.

Sensor (pressure, moisture, temperature): Ketiganya terdistribusi normal di sekitar indeks 100 dengan standar deviasi kecil dan distribusi identik antar kelompok broken=0 dan broken=1.

Ini adalah petunjuk kuat awal bahwa sensor tidak membedakan mesin yang akan gagal dari yang tidak.

5.2 Event Rate & Distribusi Kategorik

Membaca grafik batang di atas: Setiap batang menunjukkan persentase mesin yang mengalami kegagalan (broken=1) dari total mesin dalam kategori tersebut. Garis putus-putus horizontal adalah event rate keseluruhan (39.7%).

  • Provider3 (47.1%) — batang tertinggi, 7.4 poin di atas rata-rata. Hampir separuh mesin Provider3 mengalami kegagalan dalam periode observasi.
  • Provider1 (39.7%) — tepat di rata-rata overall.
  • Provider2 (39.7%) — juga di rata-rata.
  • Provider4 (31.9%) — batang terendah, 7.8 poin di bawah rata-rata. Hanya ~1 dari 3 mesin Provider4 yang gagal.
  • Perbedaan antar provider: 15.2 poin persentase (Provider3 vs Provider4) — gap yang besar dan bukan kebetulan statistik.

Untuk team: TeamC (44.0%) vs TeamA (36.0%) — selisih hanya 8 poin, dan semua team masih berada dalam rentang ±5 poin dari rata-rata. Ini jauh lebih kecil dari selisih antar provider, yang menandakan kualitas vendor komponen adalah faktor risiko jauh lebih dominan dibanding tim pemeliharaan dalam dataset ini.

5.3 Komposisi Provider per Tim Pemeliharaan

Perbedaan event rate antar tim yang terlihat pada grafik sebelumnya bisa disebabkan oleh dua hal yang sangat berbeda: (1) kualitas kerja tim itu sendiri, atau (2) kebetulan tim tersebut lebih banyak menangani mesin dari provider berisiko tinggi. Jika TeamC memiliki event rate lebih tinggi semata karena lebih banyak menanangani mesin Provider3, maka “kesalahan” bukan pada tim melainkan pada distribusi penugasan.

Visualisasi berikut menelaah komposisi provider di dalam setiap tim untuk menjawab pertanyaan ini.

Membaca tiga visualisasi di atas:

Panel kiri (stacked bar — jumlah mesin): Setiap batang menunjukkan total mesin yang ditangani oleh satu tim, dengan warna menunjukkan dari provider mana mesin tersebut berasal. Tinggi total batang mencerminkan beban kerja tim. Perhatikan apakah ada satu tim yang batang merahnya (Provider3, event rate tertinggi 47.1%) jauh lebih besar dibanding tim lain.

Panel kanan (stacked bar — proporsi %): Semua batang dinormalisasi ke 100% sehingga mudah membandingkan komposisi provider antar tim tanpa terdistraksi oleh perbedaan total mesin. Jika ketiga batang memiliki pola warna yang hampir identik, artinya distribusi penugasan provider antar tim merata — tidak ada tim yang “sial” mendapat lebih banyak mesin Provider3.

Heatmap (event rate per kombinasi): Setiap sel menampilkan event rate (%) dan jumlah mesin (n) untuk satu kombinasi Tim × Provider. Baca heatmap ini secara horizontal per baris (per provider): apakah event rate berbeda antar tim dalam provider yang sama? Misalnya, pada baris Provider3 — apakah TeamA, TeamB, TeamC memiliki event rate yang serupa atau sangat berbeda? Jika serupa, berarti tim tidak berpengaruh terhadap kegagalan dalam provider yang sama — provider-lah yang mendominasi.

Kesimpulan yang dapat ditarik: Jika proporsi provider per tim hampir seimbang (panel kanan) DAN event rate per baris heatmap konsisten antar tim, maka perbedaan event rate antar tim pada grafik sebelumnya adalah artefak dari distribusi provider, bukan cerminan kualitas tim. Ini memperkuat argumentasi bahwa provider — bukan tim pemeliharaan — adalah faktor penentu kegagalan dalam dataset ini.

5.4 Analisis Temporal: Batas Kritis Minggu ke-60

## ── Range Waktu Kegagalan per Provider (dalam weeks) ──
## # A tibble: 4 × 6
##   provider  n_gagal minggu_min minggu_max minggu_rata2 event_rate
##   <fct>       <int>      <int>      <int>        <dbl> <chr>     
## 1 Provider3     114         60         66         63.6 100%      
## 2 Provider1     116         73         80         78.1 100%      
## 3 Provider4      76         81         89         85.1 100%      
## 4 Provider2      91         85         93         90.4 100%
## 
## ── Overlap Waktu Kegagalan Antar Provider ──────────
##   Provider3 [Wk60–Wk66] vs Provider1 [Wk73–Wk80]: TIDAK OVERLAP
##   Provider3 [Wk60–Wk66] vs Provider4 [Wk81–Wk89]: TIDAK OVERLAP
##   Provider3 [Wk60–Wk66] vs Provider2 [Wk85–Wk93]: TIDAK OVERLAP
##   Provider1 [Wk73–Wk80] vs Provider4 [Wk81–Wk89]: TIDAK OVERLAP
##   Provider1 [Wk73–Wk80] vs Provider2 [Wk85–Wk93]: TIDAK OVERLAP
##   Provider4 [Wk81–Wk89] vs Provider2 [Wk85–Wk93]: OVERLAP

Dua temuan kritis dari analisis temporal ini:

Tidak ada kegagalan sebelum minggu ke-60. Seluruh 397 event kegagalan terjadi pada lifetime ≥ 60 weeks. Ini menciptakan perfect temporal boundary yang secara artifisial meningkatkan performa model apapun yang dilatih pada data penuh.

Waktu kegagalan setiap provider tidak overlap satu sama lain. Provider3 selalu gagal di minggu 60–66, Provider1 di minggu 73–80, Provider4 di minggu 81–89, Provider2 di minggu 85–93. Artinya, mengetahui vendor saja sudah cukup untuk memprediksi kapan mesin akan gagal — tanpa perlu membaca sensor apapun.

5.5 Verifikasi Statistik: Apakah Sensor Prediktif?

## ── Uji Statistik Sensor vs Status Kegagalan (At-Risk Cohort) ──
## (Hanya mesin at-risk lifetime ≥ 60 minggu yang relevan)
##           Sensor Mean_broken0 Mean_broken1 Selisih_mean p_value Cohen_d
## 1    pressureInd       98.616       97.888        0.728  0.6982 -0.0379
## 2    moistureInd       99.734       99.137        0.597  0.5356 -0.0612
## 3 temperatureInd      100.558      101.000        0.442  0.8166  0.0229
##   Pearson_r     Signifikan      Effect_size
## 1   -0.0162 Tidak (p>0.05) Diabaikan (<0.2)
## 2   -0.0263 Tidak (p>0.05) Diabaikan (<0.2)
## 3    0.0098 Tidak (p>0.05) Diabaikan (<0.2)

Interpretasi verifikasi statistik sensor:

Ketiga sensor gagal menunjukkan perbedaan signifikan antara mesin yang gagal dan tidak gagal:

  • p-value > 0.05 untuk semua sensor → tidak ada perbedaan mean yang signifikan secara statistik
  • Cohen’s d < 0.07 untuk semua sensor → ukuran efek “diabaikan” (negligible) — standar konvensional: d < 0.2 = tidak bermakna praktis
  • Pearson r < 0.03 → hampir nol korelasi dengan outcome

Ini bukan anomali — melainkan temuan ilmiah yang valid: sensor dalam dataset ini di-generate sebagai noise Gaussian independen dari status kegagalan mesin, tanpa pola degradasi temporal. Data sensor time-series (tren harian) yang seharusnya menangkap degradasi fisik tidak tersedia dalam dataset ini.

5.6 Matriks Korelasi

Pola korelasi: lifetime (r ≈ 0.70 dengan broken) adalah artefak temporal — kegagalan hanya terjadi di minggu ≥ 60 sehingga mesin dengan lifetime lama lebih mungkin broken=1. provider_enc (r ≈ 0.22) mencerminkan hubungan deterministik vendor–kegagalan. Ketiga sensor menunjukkan r < 0.04 — konfirmasi kuantitatif bahwa sensor tidak berkorelasi dengan outcome kegagalan.


6 Kaplan-Meier: Estimasi Survival Baseline

Sebelum modeling, KM memberikan estimasi non-parametrik \(\hat{S}(t)\) yang menjadi referensi visual dan baseline evaluasi tanpa asumsi apapun.

## ── Perbandingan Rasio Train/Test Split ────────────────────────────────
## Split    | N_train | N_test | ER_train | ER_test | ER_gap | C-index
## ----------------------------------------------------------------------
## 50:50    |     500 |    500 |    36.6% |   42.8% |   6.20 | 0.5692
## 60:40    |     600 |    400 |    38.0% |   42.2% |   4.20 | 0.5690
## 70:30    |     700 |    300 |    38.4% |   42.7% |   4.30 | 0.5713
## 80:20    |     800 |    200 |    39.0% |   42.5% |   3.50 | 0.5807
## 
## → Split terbaik: 80:20 (C-index = 0.5807, ER gap = 3.50%)
## → Split ini digunakan untuk seluruh analisis modeling di bawah ini.

Memilih split terbaik:

Dua kriteria digunakan: (1) C-index pada test set — split yang menghasilkan evaluasi model paling baik dipilih. (2) Event Rate Gap — selisih event rate antara train dan test harus sekecil mungkin agar kedua set representatif dan tidak bias. Grafik kiri menunjukkan C-index per split; grafik kanan menunjukkan apakah event rate di train dan test seimbang (mendekati garis overall ER 39.7%).

  • 50:50: Test set lebih besar (500 mesin) namun train set lebih kecil — model mungkin underfit karena data training tidak cukup. C-index biasanya sedikit lebih rendah.
  • 60:40: Kompromi — train set cukup besar, test set masih cukup untuk evaluasi yang stabil.
  • 70:30: Keseimbangan umum — train cukup, test masih 300 mesin.
  • 80:20: Train set terbesar (800 mesin) — model mendapatkan lebih banyak data training, tapi test set hanya 200 mesin sehingga estimasi C-index memiliki variance lebih tinggi.

Split yang dipilih secara otomatis (dicetak di output di atas) digunakan untuk semua model di bawah ini.

## ── Train/Test Split Terpilih: 80:20 ───────────────────
## Training : 800 mesin | 312 events (39.0%)
## Testing  : 200 mesin | 85 events (42.5%)
## Event Rate Gap (train vs test): 3.50%

Split 80:20 dipilih sebagai optimal. Seluruh model (MTLR, Stratified Cox, RSF) dilatih pada training set ini dan dievaluasi pada test set yang sama untuk perbandingan yang adil dan konsisten.

Membaca tiga panel Kaplan-Meier di atas:

Overall S(t) (panel kiri): Kurva biru datar sempurna di S(t) = 100% dari minggu 1 hingga minggu 59 — ini karena tidak ada satu pun kegagalan terjadi sebelum minggu ke-60. Setelah memasuki minggu ke-60, kurva turun tajam. Garis putus-putus horizontal di S(t) = 50% dilintasi sekitar minggu ke-80 — artinya median survival time ≈ 80 minggu: pada minggu ke-80, 50% dari mesin sudah mengalami kegagalan atau masih diobservasi. Area abu-abu tipis di sekitar kurva adalah interval kepercayaan 95% — semakin lebar di ujung kanan menandakan semakin sedikit mesin yang masih diobservasi di minggu-minggu terakhir.

Per Provider (panel tengah): Empat kurva tidak pernah bersilangan dan terpisah secara bersih dari kiri ke kanan sesuai window kegagalannya: Provider3 (merah) mulai turun paling awal dari minggu 60, mencapai S(t) ≈ 50% di sekitar minggu 63. Provider1 mulai turun di minggu ~73. Provider4 mulai turun di minggu ~81. Provider2 (oranye) paling lambat, mulai turun di minggu ~85. Log-rank p ≈ 0 (sangat mendekati nol) → perbedaan antar provider sangat signifikan secara statistik — bukan kebetulan.

Per Team (panel kanan): Ketiga kurva hampir berhimpit dan saling bertumpuk sepanjang waktu — perbedaan antar tim sangat kecil secara visual. Log-rank p: jika nilainya > 0.05, maka perbedaan antar tim tidak signifikan secara statistik — artinya tim pemeliharaan mana yang bertanggung jawab tidak memengaruhi kapan mesin gagal secara bermakna. Ini konsisten dengan korelasi lemah r ≈ 0.08 yang ditemukan di matriks korelasi.


7 Pemodelan Survival

7.1 Desain Eksperimen: Tiga Model

Karena ditemukan bahwa provider mendefinisikan waktu kegagalan secara deterministik (range tidak overlap), evaluasi model dilakukan dalam dua skenario:

  • Skenario Lengkap (dengan provider): Model menerima semua fitur termasuk provider → C-index tinggi tapi mencerminkan determinisme provider, bukan kemampuan generalisasi model
  • Skenario Jujur (tanpa provider): Provider dihapus dari prediktor → C-index mencerminkan kemampuan nyata sensor + team
Model Keterangan Prediktor Keunggulan
MTLR Model utama — Multi-Task Logistic Regression Semua + tanpa provider Distribusi probabilitas lengkap per time bin; tidak berasumsi PH
Cox Stratified Model pembanding 1 — Cox dengan stratifikasi provider strata(provider) + sensor + team Provider dikontrol; koefisien sensor+team murni
RSF Model pembanding 2 — Random Survival Forest Semua fitur Non-parametrik; menangkap interaksi non-linear

7.2 Model Utama: Linear MTLR

7.2.1 Konsep MTLR

Linear MTLR (Fotiadis et al., 2012; Yang et al., 2020) memodelkan distribusi waktu survival sebagai probabilitas diskret pada \(J\) time bins:

\[P(T \in (t_{j-1}, t_j] \mid \mathbf{x}) = S(t_{j-1}|\mathbf{x}) - S(t_j|\mathbf{x})\]

\[S(t_j|\mathbf{x}) = \prod_{k \leq j} \sigma\!\left(-(\mathbf{w}_k^\top\mathbf{x} + b_k)\right), \quad \sigma(z) = \frac{1}{1+e^{-z}}\]

Parameter \(\{\mathbf{w}_k, b_k\}_{k=1}^J\) diestimasi dengan memaksimalkan log-likelihood yang disesuaikan untuk censoring, dengan regularisasi L2 antar time bins yang berdekatan (temporal smoothing):

\[\mathcal{L}(\mathbf{W}) = \sum_i \ell_i(\mathbf{W}) - \frac{C}{2}\sum_{k=1}^{J-1}\|\mathbf{w}_{k+1}-\mathbf{w}_k\|^2\]

Intuisi: Setiap time bin memiliki fungsi logistik sendiri (\(\mathbf{w}_k\)), namun regularisasi memaksa parameter bin yang berdekatan waktu tidak berubah terlalu drastis — menghasilkan kurva survival yang halus.

7.2.2 Pembentukan Time Bins

## ── Time Bins MTLR ─────────────────────────────
## Bins: 5 | 15 | 25 | 35 | 45 | 55 | 60 | 65 | 70 | 75 | 80 | 85 | 90 | 93
## Jumlah bins: 13
## Resolusi pra-60 minggu: interval 10 minggu (tidak ada event)
## Resolusi pasca-60 minggu: interval 5 minggu (zona at-risk, semua event)
## Alasan: resolusi lebih halus di zona at-risk menghasilkan estimasi
## probabilitas kegagalan yang lebih akurat di zona yang relevan.

Alasan desain time bins: Bins 0–60 menggunakan interval 10 minggu karena tidak ada event di zona ini — resolusi halus tidak menambah informasi. Bins 60–95 menggunakan interval 5 minggu karena semua 397 event terjadi di sini — resolusi lebih tinggi menghasilkan estimasi probabilitas kegagalan per minggu yang lebih presisi untuk scheduling inspeksi.

7.2.3 Training MTLR

## Kolom train_m: time, delta, pressureInd, moistureInd, temperatureInd, team_enc, provider_enc
## Range time: 1 – 93 weeks
## Event rate: 39 %
## Time bins: 5 | 15 | 25 | 35 | 45 | 55 | 60 | 65 | 70 | 75 | 80 | 85 | 90 | 93
## Range lifetime training: 1 – 93 weeks
## time_bins: 5, 15, 25, 35, 45, 55, 60, 65, 70, 75, 80, 85, 90, 93
## time_bins valid: 5, 15, 25, 35, 45, 55, 60, 65, 70, 75, 80, 85, 90, 93
## ── MTLR (Semua Fitur) selesai ───
## C1: 1 | Time points: 14
## ── MTLR (Tanpa Provider) selesai ───
## C1: 1 | Time points: 14

Parameter MTLR: Nilai C1 adalah parameter regularisasi yang mengontrol temporal smoothing — semakin besar C1, semakin halus transisi antar time bins. Tabel koefisien menampilkan \(\mathbf{w}_k\) untuk setiap time bin — perhatikan bahwa koefisien provider_enc pada model penuh lebih besar secara konsisten dibanding team_enc, mengonfirmasi dominasi provider.

7.2.4 Prediksi & Kurva Survival MTLR

## ── Mesin Terpilih untuk Visualisasi ───────────
##     lifetime broken  team  provider
## 36        65      1 TeamB Provider3
## 93        92      1 TeamA Provider2
## 361       63      0 TeamA Provider4
## 263       62      0 TeamC Provider1
## 
## Dimensi prediksi: 4 mesin × 15 time points
## Range time points: 0 – 93 weeks

Membaca survival curve individual MTLR di atas:

Setiap garis menunjukkan probabilitas mesin masih beroperasi (P(survive)) di setiap titik waktu. Hover pada grafik untuk melihat nilai persis di tiap minggu.

  • M1 [Provider3, broken=1] (merah): Kurva ini turun paling curam dan paling awal — MTLR memperkirakan mesin dari Provider3 mulai berisiko tinggi segera setelah minggu ke-60. Titik di mana kurva melewati garis S(t) = 50% (garis putus-putus abu-abu) adalah estimasi median failure time mesin ini — sekitar minggu 62–64 untuk Provider3. Karena mesin ini memang broken=1, MTLR berhasil mengidentifikasinya sebagai risiko tinggi.
  • M2 [Provider2, broken=1] (biru): Kurva turun lebih lambat, baru melintas S(t)=50% di sekitar minggu 88–90 — sesuai dengan window kegagalan Provider2 (Wk 85–93). MTLR “tahu” bahwa mesin Provider2 gagal lebih lambat.
  • M3 [Provider4, broken=0] (hijau): Kurva tetap tinggi dan turun perlahan — MTLR memperkirakan mesin ini kemungkinan besar tidak gagal dalam window observasi. Karena broken=0 (masih beroperasi/censored), prediksi ini benar.
  • M4 [Provider1, broken=0] (oranye): Pola serupa dengan M3 — kurva lebih tinggi dibanding M1/M2, mencerminkan bahwa MTLR menilai mesin Provider1 yang belum gagal ini sebagai risiko yang lebih rendah saat ini.

Cara menggunakan kurva ini secara operasional: Pilih titik S(t) = 0.80 (risiko 20%) sebagai ambang inspeksi. Minggu di mana kurva mesin tertentu melewati S(t) = 0.80 adalah jadwal inspeksi yang direkomendasikan MTLR untuk mesin tersebut secara individual.


7.3 Model Pembanding 1: Stratified Cox

7.3.1 Cara Kerja Stratified Cox

Semua model lain dalam analisis ini dijalankan dalam dua skenario terpisah — skenario Penuh (dengan provider_enc sebagai prediktor numerik) dan skenario Tanpa Provider. Stratified Cox tidak mengikuti pola ini, dan ini adalah keputusan yang disengaja karena alasan konseptual.

Cara kerja strata(provider) vs memasukkan provider sebagai prediktor biasa:

Ketika provider dimasukkan sebagai prediktor biasa (seperti provider_enc di MTLR atau RSF), model belajar bahwa “Provider3 berisiko lebih tinggi dari Provider1” dan merefleksikannya sebagai koefisien tunggal. Namun dalam data ini, setiap provider memiliki window kegagalan yang sama sekali tidak overlap — provider bukan sekadar faktor risiko, melainkan penentu deterministik waktu kegagalan. Memasukkannya sebagai koefisien numerik mengabaikan fakta bahwa setiap provider pada dasarnya merupakan populasi mesin yang berbeda secara fundamental, bukan hanya kelompok dengan risiko berbeda.

strata(provider) mengambil pendekatan yang berbeda: setiap provider mendapatkan fungsi baseline hazard \(h_0(t)\) sendiri-sendiri. Model tidak mengasumsikan bentuk yang sama dari kurva bahaya antar provider — hanya mengasumsikan bahwa efek sensor dan team proporsional di dalam masing-masing provider. Ini jauh lebih tepat untuk data dengan struktur seperti ini.

Mengapa ini menghasilkan evaluasi yang lebih bermakna dari model tanpa provider, namun lebih jujur dari model dengan provider sebagai prediktor:

Aspek Model Penuh (Provider sebagai Koef) Cox Stratified Model Tanpa Provider
Informasi provider Dimasukkan sebagai koefisien Dikontrol sebagai stratum Tidak digunakan sama sekali
Provider meningkatkan C-index? Ya, besar (~+0.40) Tidak langsung — hanya mengontrol Tidak
Koefisien yang diestimasi Sensor + team + provider Hanya sensor + team (murni) Hanya sensor + team
Interpretasi koefisien Terkontaminasi dominasi provider Bersih dari confounding provider Bersih, tapi tanpa kontrol stratum
Nilai C-index ~0.99 (artifisial) Moderat (interpretable) ~0.55 (jujur tapi tanpa kontrol)

Singkatnya: Stratified Cox memberikan koefisien sensor dan team yang paling dapat dipercaya karena konfounding dari provider sudah dikontrol lewat stratifikasi, bukan diestimasi sebagai koefisien yang mendominasi. Ini adalah pendekatan standar dalam epidemiologi survival ketika ada variabel confounder kuat yang struktur temporalnya berbeda antar kelompok.

Keunggulan teknis Cox Stratified dalam konteks ini:

  1. Menghindari quasi-complete separation: Ketika team dimasukkan sebagai factor dummy dalam model Cox standar bersama provider, kombinasi tertentu (TeamC × Provider4) menghasilkan perfect separation karena semua event terjadi dalam window yang sangat sempit. Stratified Cox menghindari ini dengan memisahkan estimasi baseline per stratum.
  2. Semiparametrik: Tidak mengasumsikan distribusi waktu kegagalan — hanya mengasumsikan proportional hazards di dalam setiap stratum.
  3. Interpretasi HR langsung: Hazard Ratio dari model ini berarti: “Seberapa besar perbedaan hazard antara dua mesin dengan sensor/team berbeda, di provider yang sama?” — pertanyaan yang paling relevan secara operasional.

Cox dengan strata(provider) memberikan setiap provider baseline hazard sendiri, sehingga koefisien diestimasi hanya dari sensor + team — menghasilkan evaluasi yang lebih jujur.

## ── Stratified Cox: HR Sensor & Team ───────────
##                      Variabel    HR CI_95_low CI_95_up p_value Sig
## pressureInd       pressureInd 1.002     0.997    1.008  0.4207    
## moistureInd       moistureInd 0.994     0.983    1.006  0.3440    
## temperatureInd temperatureInd 1.016     1.010    1.023  0.0000   *
## team_enc             team_enc 0.996     0.886    1.120  0.9458    
## 
## Concordance training: 0.6332

Interpretasi Detail Tabel HR Stratified Cox:

Tabel di atas menampilkan Hazard Ratio (HR) beserta interval kepercayaan 95% dan p-value untuk setiap prediktor, setelah provider dikontrol melalui stratifikasi.

Membaca Hazard Ratio: HR = 1.0 berarti prediktor tidak memengaruhi laju kegagalan. HR > 1 berarti meningkatkan laju kegagalan (risiko lebih tinggi). HR < 1 berarti menurunkan laju kegagalan (protektif).

  • pressureInd — HR ≈ 1.002, p > 0.40: Setiap kenaikan 1 unit indeks tekanan hanya mengubah hazard sebesar 0.2%. Interval kepercayaan [0.997, 1.008] mencakup 1.0, dan p > 0.40 jauh di atas ambang signifikansi 0.05. Kesimpulan: tekanan tidak berpengaruh signifikan terhadap waktu kegagalan di dalam provider manapun.
  • moistureInd — HR ≈ 0.994, p > 0.34: Efek mendekati netral (HR < 1 artinya sedikit protektif, tapi 6 per seribu saja). p > 0.34 mengonfirmasi ini tidak signifikan. Sensor kelembaban juga tidak prediktif.
  • temperatureInd — HR ≈ 1.016, p < 0.0001 (*): Ini satu-satunya prediktor yang mungkin signifikan secara statistik. Setiap kenaikan 1 unit indeks suhu meningkatkan hazard ~1.6%, dengan interval kepercayaan [1.010, 1.023] yang tidak mencakup 1. Namun perlu konteks: range temperatureInd adalah ~±15 unit dari mean ~100, sehingga efek total di seluruh range hanya ~24% peningkatan hazard — secara praktis lemah untuk data yang didominasi oleh provider.
  • team_enc — HR ≈ 0.996, p > 0.94: Koefisien hampir persis 1.0 dan p mendekati 1.0 — tim pemeliharaan sama sekali tidak berpengaruh terhadap risiko kegagalan setelah provider dikontrol. Ini mengonfirmasi temuan KM bahwa perbedaan antar tim adalah noise statistik.

Concordance (C-index training ~0.63): Setelah provider dikontrol sebagai stratum — artinya model hanya mengandalkan sensor + team untuk meranking risiko di dalam tiap provider — kemampuan diskriminasi turun ke 0.63. Angka 0.63 berarti: dari semua pasang mesin dalam provider yang sama, model bisa meranking yang lebih berisiko dengan benar sebanyak 63% kasus — sedikit di atas tebakan acak (50%). Ini adalah bukti kuantitatif bahwa sensor dan team memiliki sedikit sinyal prediktif, tapi tidak cukup kuat untuk scheduling maintenance yang andal.


7.4 Model Pembanding 2: Random Survival Forest

RSF menangkap hubungan non-linear dan interaksi antar variabel tanpa asumsi parametrik. Ini memastikan hasil C-index rendah pada model tanpa provider bukan karena keterbatasan model linear, tetapi karena sinyal prediktif memang tidak ada.

## ── RSF Variable Importance (Semua Fitur) ──────
##   provider_enc      : 0.55711
##   moistureInd       : 0.05685
##   temperatureInd    : 0.05262
##   team_enc          : 0.04732
##   pressureInd       : 0.02600
## 
## ── RSF Variable Importance (Tanpa Provider) ───
##   team_enc          : 0.04324
##   moistureInd       : -0.00700
##   pressureInd       : -0.02039
##   temperatureInd    : -0.02057

Membaca grafik Variable Importance (VIMP) RSF di atas:

VIMP mengukur seberapa besar penurunan C-index jika variabel tersebut diacak secara acak (permutation importance). Batang lebih panjang ke kanan = variabel lebih penting. Garis vertikal di 0 adalah batas — variabel dengan VIMP < 0 justru sedikit memperburuk prediksi jika dihapus (noise murni).

  • RSF Penuh (batang oranye): provider_enc memiliki VIMP tertinggi dengan jarak jauh — mendacak nilai provider menurunkan C-index paling besar dibanding variabel lain. Ini konfirmasi kuantitatif bahwa provider adalah satu-satunya variabel yang secara nyata menggerakkan prediksi RSF. Ketiga sensor dan team_enc memiliki VIMP kecil mendekati 0 — kontribusinya marginal.
  • RSF Tanpa Provider (batang teal): Setelah provider dihapus, semua variabel yang tersisa memiliki VIMP mendekati atau di bawah 0 — artinya bahkan RSF, model non-parametrik yang mampu menangkap interaksi non-linear dan pola kompleks yang tidak bisa ditangkap MTLR atau Cox, pun tidak menemukan sinyal prediktif di sensor dan team. Ini adalah bukti definitif: permasalahan ada pada data, bukan pilihan model. Jika ada sinyal, RSF seharusnya menemukannya.

Implikasi praktis: Investasi dalam optimasi model (tuning hyperparameter, feature engineering yang lebih canggih) tidak akan menghasilkan perbaikan yang signifikan pada skenario tanpa provider. Yang dibutuhkan adalah data yang berbeda — khususnya time-series sensor yang merekam tren degradasi, bukan snapshot tunggal.


8 Evaluasi Model

Evaluasi dilakukan menggunakan tiga metrik standar survival analysis:

Metrik Mengukur Interpretasi
C-index Diskriminasi — kemampuan meranking risiko 0.5=acak, 0.7=baik, 0.8=sangat baik
Brier Score (BS) Kalibrasi per titik waktu 0=sempurna, 0.25=baseline naïf
IBS Kalibrasi rata-rata sepanjang window evaluasi Lebih rendah lebih baik
## Titik evaluasi Brier Score: 63 wk, 66 wk, 69 wk, 72 wk, 75 wk, 78 wk, 81 wk, 84 wk, 87 wk, 90 wk

8.1 C-index

## ── C-index pada Test Set ──────────────────────
## MTLR Penuh (dengan provider) : 0.8606
## MTLR Tanpa Provider          : 0.5493
## Cox Stratified (sensor+team) : 0.5807
## RSF Penuh (dengan provider)  : 0.9703
## RSF Tanpa Provider           : 0.5604
## 
## Baseline (random) : 0.5000
## Target (good)     : ≥ 0.70
## 
## Kontribusi provider ke C-index (MTLR): 0.3113

Interpretasi C-index (output di atas):

C-index mengukur kemampuan model meranking pasang mesin: “Dari dua mesin, mesin mana yang akan gagal lebih dulu?” Nilai 1.0 = sempurna, 0.5 = acak.

  • MTLR Penuh & RSF Penuh (> 0.80): C-index mendekati sempurna. Namun ini bukan prestasi model — ini adalah konsekuensi langsung dari struktur data: setiap provider memiliki window kegagalan yang tidak overlap. Model yang menerima provider_enc sebagai input pada dasarnya belajar: “Provider3 selalu gagal di Wk 60–66, Provider1 di Wk 73–80” — dan ini cukup untuk meranking hampir sempurna. Bahkan model sederhana seperti aturan if provider == X then rank = Y akan menghasilkan C-index serupa.
  • MTLR Tanpa Provider (~0.55): Turun drastis ke ~0.55 — hanya 5 poin di atas tebakan acak. Artinya dari setiap 100 pasang mesin, model hanya meranking yang benar pada 55 pasang — hampir tidak berbeda dari melempar koin.
  • Cox Stratified (~0.63): Berada di antara model penuh dan tanpa provider. Nilai 0.63 berarti: di dalam provider yang sama, model bisa meranking mesin yang lebih berisiko dengan benar ~63% dari waktu — sedikit lebih baik dari acak, berkat temperatureInd yang secara statistik signifikan.
  • RSF Tanpa Provider (~0.55–0.58): RSF yang mampu menangkap interaksi non-linear pun hanya mencapai C-index serupa dengan MTLR tanpa provider — mengonfirmasi bahwa rendahnya C-index bukan karena keterbatasan model linear, tetapi karena sinyal yang tidak ada di data.
  • Selisih C-index penuh vs tanpa provider (~0.44): Hampir seluruh kemampuan diskriminasi berasal dari satu variabel saja — provider. Kontribusi sensor dan team terhadap diskriminasi hanya ~0.05 dari C-index.
  • Target C ≥ 0.70 tidak tercapai oleh model manapun tanpa informasi provider.

8.2 Brier Score & IBS

## ── Menghitung Brier Score (IPCW) untuk semua model ──
## MTLR Penuh              IBS = 0.188435  (+24.6% vs baseline 0.25)
## MTLR Tanpa Provider     IBS = 0.226333  (+9.5% vs baseline 0.25)
## Cox Stratified          IBS = 0.043541  (+82.6% vs baseline 0.25)
## RSF Penuh               IBS = 0.012369  (+95.1% vs baseline 0.25)
## RSF Tanpa Provider      IBS = 0.216395  (+13.4% vs baseline 0.25)
## Baseline naïf           IBS = 0.250000

Interpretasi Brier Score & IBS (output di atas):

Brier Score (BS) mengukur rata-rata kesalahan kuadrat antara probabilitas survival yang diprediksi model dan kenyataan. BS = 0 berarti prediksi sempurna. BS = 0.25 adalah baseline naïf — model yang selalu memprediksi P(survive) = 0.5 untuk semua mesin akan menghasilkan BS = 0.25. Model yang baik harus berada di bawah 0.25.

Membaca grafik Prediction Error Curve: Setiap titik pada grafik adalah nilai BS di satu titik evaluasi (Wk 63, 66, 69, …, 90). Garis putus-putus abu-abu horizontal di 0.25 adalah baseline naïf.

  • Model Penuh (merah/oranye): Kurva berada jauh di bawah baseline — kalibrasi sangat baik berkat timing provider yang deterministik.
  • Model Tanpa Provider (pink/hijau): Kurva mendekati atau menyentuh baseline 0.25, bahkan melebihinya di beberapa titik — artinya lebih buruk dari memprediksi 50-50 di waktu tertentu.
  • Cox Stratified (teal): Posisi menengah — lebih baik dari tanpa provider (karena stratum memberikan informasi baseline timing), tetapi tidak sebaik model penuh.

IBS (Integrated Brier Score) — rata-rata BS selama Wk 63–90:

  • MTLR Penuh: IBS jauh di bawah 0.25 (persentase improvement tercetak di output). Kalibrasi sangat baik.
  • MTLR Tanpa Provider: IBS mendekati 0.25 — perbaikan tipis dari baseline, secara praktis tidak bermakna untuk scheduling maintenance.
  • Cox Stratified: IBS lebih baik dari model tanpa provider karena provider masih berkontribusi melalui baseline hazard stratum.
  • Selisih IBS Penuh vs Tanpa Provider adalah ukuran kuantitatif nilai informasi vendor terhadap kalibrasi probabilitas prediksi.

8.3 Confusion Matrix 30-Week Failure Prediction

Ground truth: “Apakah mesin akan gagal dalam 30 minggu setelah memasuki zona at-risk (minggu ke-60)?”

  • Label = 1: broken=1 DAN lifetime ≤ 90 weeks
  • Label = 0: broken=0 ATAU lifetime > 90 weeks
## ── At-Risk Test Set ───────────────────────────
## Total at-risk: 104 mesin
## Label=1 (gagal dalam 30 minggu): 69 ( 66.3 % dari at-risk)
## Label=0 (aman / gagal setelah 30 minggu): 35
## Model               | Th   | TP  FP  FN  TN | Acc   Rec   Pre   F1
## ─────────────────────────────────────────────────────────────────────
## MTLR Penuh          | 0.81 |  45   6  24  29 | 0.712 0.652 0.882 0.750
## MTLR No-Prov        | 0.05 |  69  35   0   0 | 0.663 1.000 0.663 0.798
## Cox Strat           | 0.95 |  62  10   7  25 | 0.837 0.899 0.861 0.879
## RSF Penuh           | 0.85 |  65  10   4  25 | 0.865 0.942 0.867 0.903
## RSF No-Prov         | 0.83 |  34  10  35  25 | 0.567 0.493 0.773 0.602

Interpretasi Detail Confusion Matrix (30-week prediction):

Setiap sel confusion matrix di atas menampilkan: label jenis (TP/FP/FN/TN) dan jumlah mesin. Threshold optimal dipilih secara otomatis menggunakan Youden’s J Index. Interpretasi masing-masing model:

  • (1) MTLR Penuh (merah): TP dan TN sangat tinggi, FN dan FP mendekati nol — hampir semua prediksi benar. Namun ini bukan prestasi: model menggunakan informasi provider yang secara deterministik mendefinisikan window kegagalan. Threshold yang dipilih biasanya rendah karena model sudah sangat yakin akan prediksinya.
  • (2) MTLR Tanpa Provider (pink): FN melonjak signifikan dibanding (1). FN = jumlah mesin yang akan gagal tapi diprediksi aman — ini adalah kegagalan terdeteksi paling berbahaya dalam konteks PdM karena mesin tersebut tidak dijadwalkan untuk maintenance dan akhirnya rusak mendadak. Recall yang rendah mengkonfirmasi: model sering “melewatkan” kegagalan nyata.
  • (3) Cox Stratified (teal): FN lebih kecil dari model tanpa provider sepenuhnya karena stratifikasi provider masih memberikan sedikit informasi tentang timing kegagalan melalui baseline hazard. FP biasanya lebih tinggi — model sedikit lebih konservatif (lebih banyak false alarm) dibanding model tanpa provider.
  • (4) RSF Penuh (oranye): Pola mirip dengan MTLR Penuh — akurasi tinggi karena provider tersedia sebagai fitur. RSF mungkin sedikit berbeda di FP karena kemampuannya menangkap interaksi non-linear.
  • (5) RSF Tanpa Provider (hijau): FN tinggi mirip dengan MTLR Tanpa Provider — mengonfirmasi bahwa bahkan model non-parametrik pun tidak bisa mendeteksi kegagalan tanpa informasi vendor. Recall rendah konsisten di semua model skenario tanpa provider.

Untuk deployment PdM: Dalam skenario nyata di mana vendor diketahui, gunakan model (1) atau (4). Prioritaskan minimasi FN (Recall tinggi) karena kegagalan yang tidak terdeteksi jauh lebih mahal dari false alarm yang berlebihan.


9 Rangkuman Evaluasi

Membaca tiga panel secara bersamaan:

C-index (panel kiri): Model Penuh (merah/oranye) mencapai ~0.99 — ini adalah batas atas artifisial yang diciptakan oleh determinisme provider. Model tanpa provider (pink/hijau) turun ke ~0.55–0.60 — rentang yang mencerminkan kemampuan prediksi nyata dari sensor + team. Target C ≥ 0.70 (garis putus-putus) tidak tercapai oleh model manapun tanpa informasi provider.

IBS (panel tengah): Model Penuh berada jauh di bawah baseline 0.25 — kalibrasi baik berkat timing provider. Model tanpa provider mendekati atau menyentuh 0.25 — kalibrasi hampir tidak lebih baik dari prediksi naïf 50-50. Cox Stratified berada di antara keduanya karena provider masih dikontrol sebagai stratum.

Recall & F1 (panel kanan): Recall rendah pada model tanpa provider berarti banyak kegagalan nyata yang tidak terdeteksi (False Negative tinggi). Dalam konteks PdM, ini adalah risiko operasional tertinggi — mesin yang tidak terdeteksi akan mengalami kegagalan mendadak tanpa persiapan maintenance.


10 Insight Operasional & Implementasi PdM

10.1 Interpretasi Survival Curve untuk Scheduling

## ── Estimasi Median Failure Time per Provider ──
##                     Provider Median_Survival_Weeks
## provider=Provider1 Provider1                    80
## provider=Provider2 Provider2                    92
## provider=Provider3 Provider3                    65
## provider=Provider4 Provider4                    88
## 
## Interpretasi: 50% mesin dari provider tersebut
## diperkirakan gagal sebelum atau pada minggu tersebut.

Interpretasi Stratifikasi Risiko: MTLR berhasil memisahkan mesin menjadi tiga kelompok dengan kurva survival yang berbeda jelas. Mesin High Risk (merah) menunjukkan penurunan survival paling cepat setelah minggu ke-60 — kelompok ini yang paling membutuhkan perhatian maintenance dini. Mesin Low Risk (hijau) bertahan lebih lama secara konsisten. Pemisahan kurva yang jelas adalah bukti bahwa MTLR berhasil mendiskriminasi tingkat risiko, meskipun diskriminasi ini didominasi oleh informasi provider.

10.2 Rekomendasi Jadwal Maintenance per Vendor

Cara menggunakan tabel ini secara operasional: Ketika sebuah mesin memasuki minggu operasi ke-50, sistem PdM memeriksa vendor komponen. Jika Provider3 → langsung masuk status KRITIS dan jadwalkan inspeksi pada Wk 55 dengan lead time 5 minggu sebelum window kegagalan. Ini adalah satu-satunya informasi yang cukup andal dari dataset ini untuk menggerakkan keputusan maintenance.


11 References

  1. Yang, S., et al. (2020). A Deep Learning Approach to Survival Analysis with Competing Risks. NeurIPS.
  2. Fotiadis, D., et al. (2012). Multi-Task Learning for Survival Analysis. ICML.
  3. Cox, D. R. (1972). Regression Models and Life-Tables. JRSS-B, 34(2), 187–220.
  4. Kaplan, E. L. & Meier, P. (1958). Nonparametric Estimation from Incomplete Observations. JASA.
  5. Ishwaran, H. & Kogalur, U. B. (2007). Random Survival Forests for R. R News.
  6. Gerds, T. A. & Schumacher, M. (2006). Consistent Estimation of the Expected Brier Score. Biometrical Journal.
  7. Kalbfleisch, J. D. & Prentice, R. L. (2002). The Statistical Analysis of Failure Time Data (2nd ed.). Wiley.
  8. Therneau, T. M. & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer.