📚 Mata Kuliah: Analisis Survival
👨‍🏫 Dosen Pengampu: Bakti Siregar, M.Sc., CDS
👥 Kelompok: Churn Prediction

1 - Pendahuluan (Introduction)

Churn (pergantian) / atrisi pelanggan, alias persentase pelanggan yang berhenti menggunakan produk atau layanan perusahaan, adalah salah satu metrik paling penting bagi sebuah bisnis, karena biasanya memakan biaya lebih besar untuk mendapatkan pelanggan baru dibandingkan dengan mempertahankan pelanggan yang sudah ada.

Faktanya, menurut sebuah studi oleh Bain & Company, pelanggan yang sudah ada cenderung membeli lebih banyak dari sebuah perusahaan seiring berjalannya waktu, sehingga mengurangi biaya operasional bisnis dan mungkin merekomendasikan produk yang mereka gunakan kepada orang lain. Sebagai contoh, di layanan keuangan, peningkatan retensi pelanggan sebesar 5% dapat menghasilkan peningkatan laba lebih dari 25%.

Dengan menggunakan Analisis Survival (Survival Analysis), perusahaan tidak hanya dapat memprediksi apakah pelanggan kemungkinan besar akan berhenti berlangganan, tetapi juga kapan peristiwa tersebut mungkin terjadi.

2 - Pengaturan (Set up)

Sebuah perusahaan perangkat lunak sebagai layanan (software as a service / SaaS) menyediakan serangkaian produk untuk Usaha Kecil dan Menengah (UKM), seperti penyimpanan data, Akuntansi, manajemen Perjalanan dan Pengeluaran, serta manajemen Penggajian.

Untuk membantu CFO memprakirakan biaya akuisisi dan pemasaran untuk tahun fiskal berikutnya, tim Data Science ingin membangun model churn untuk memprediksi kapan pelanggan kemungkinan besar akan menghentikan langganan bulanan mereka. Dengan demikian, setelah pelanggan ditandai memiliki kemungkinan untuk berhenti berlangganan (churn) dalam jangka waktu tertentu, perusahaan dapat mengambil tindakan retensi yang diperlukan.

3 - Kumpulan Data (Dataset)

3.1 - Deskripsi dan Gambaran Umum (Description and Overview)

Dataset yang ingin digunakan oleh tim berisi fitur-fitur berikut:

Kategori fitur	Nama fitur	Tipe	Deskripsi
Waktu	`months_active`	numerik	Jumlah bulan sejak pelanggan memulai langganannya
Kejadian	`churned`	kategorikal	Menentukan apakah pelanggan berhenti melakukan bisnis dengan perusahaan
Produk	`product_data_storage`	numerik	Jumlah penyimpanan data cloud yang dibeli dalam Gigabyte
Produk	`product_travel_expense`	kategorikal	Menunjukkan apakah pelanggan secara aktif menggunakan dan membayar layanan manajemen Perjalanan dan Pengeluaran atau tidak. (‘Active’, ‘Free-Trial’, ‘No’)
Produk	`product_payroll`	kategorikal	Menunjukkan apakah pelanggan secara aktif menggunakan dan membayar layanan manajemen Penggajian atau tidak. (‘Active’, ‘Free-Trial’, ‘No’)
Produk	`product_accounting`	kategorikal	Menunjukkan apakah pelanggan secara aktif menggunakan dan membayar layanan Akuntansi atau tidak. (‘Active’, ‘Free-Trial’, ‘No’)
Kepuasan	`csat_score`	numerik	Customer Satisfaction Score (CSAT) adalah ukuran seberapa baik produk dan layanan yang disediakan oleh perusahaan memenuhi ekspektasi pelanggan.
Kepuasan	`minutes_customer_support`	numerik	Menit yang dihabiskan pelanggan di telepon dengan customer support perusahaan
Pemasaran	`articles_viewed`	numerik	Jumlah artikel yang dilihat pelanggan di situs web perusahaan.
Pemasaran	`smartphone_notifications_viewed`	numerik	Jumlah notifikasi ponsel pintar yang dilihat pelanggan
Pemasaran	`marketing_emails_clicked`	numerik	Jumlah email pemasaran yang dibuka pelanggan
Pemasaran	`social_media_ads_viewed`	numerik	Jumlah iklan media sosial yang dilihat pelanggan
Informasi pelanggan	`company_size`	kategorikal	Ukuran perusahaan
Informasi pelanggan	`us_region`	kategorikal	Wilayah di AS tempat kantor pusat pelanggan berada

library(readr)

# Perhatikan variabel url di bawah ini. HANYA ada tanda kutip dan alamat web murni.
url <- "https://raw.githubusercontent.com/Square/pysurvival/master/pysurvival/datasets/churn.csv"
raw_dataset <- read_csv(url)

raw_dataset <- as.data.frame(raw_dataset)

cat(sprintf("Dataset mentah memiliki bentuk berikut: %d baris, %d kolom.\n", nrow(raw_dataset), ncol(raw_dataset)))

## Dataset mentah memiliki bentuk berikut: 2000 baris, 14 kolom.

head(raw_dataset, 2)

3.2 - Dari kategorikal ke numerik (From categorical to numerical)

Ada beberapa fitur kategorikal yang perlu di-encode menjadi vektor one-hot: - product_travel_expense - product_payroll - product_accounting - us_region - company_size

# Memuat library caret secara eksplisit agar fungsi dummyVars dikenali
library(caret)

# Membuat daftar kolom kategorikal
categories <- c('product_travel_expense', 'product_payroll', 'product_accounting',
                'us_region', 'company_size')

# Mengubah tipe data menjadi factor (syarat untuk fungsi dummyVars di R)
raw_dataset[categories] <- lapply(raw_dataset[categories], as.factor)

# Membuat vektor one-hot menggunakan dummyVars dari paket caret
# Argumen fullRank = TRUE berfungsi sama persis dengan drop_first=True di Python
dummy_model <- dummyVars(" ~ .", data = raw_dataset, fullRank = TRUE)
dataset <- data.frame(predict(dummy_model, newdata = raw_dataset))

# Membuat kolom waktu (time) dan kejadian (event)
time_column <- 'months_active'
event_column <- 'churned'

# Mengekstrak fitur-fitur (mengambil semua nama kolom kecuali time_column dan event_column)
features <- setdiff(colnames(dataset), c(time_column, event_column))

4 - Analisis Data Eksploratif (Exploratory Data Analysis)

Karena tutorial ini utamanya dirancang untuk memberikan contoh tentang cara menggunakan pemodelan Survival di R, kita tidak akan melakukan analisis data eksploratif secara menyeluruh di sini, namun kami sangat menyarankan pembaca untuk melakukannya dengan melihat tutorial predictive maintenance (pemeliharaan prediktif) yang menyediakan analisis mendetail.

Di sini, kita hanya akan memeriksa apakah dataset mengandung nilai Null (kosong) atau apakah ada baris yang terduplikasi. Kemudian, kita akan melihat korelasi antar fitur.

4.1 - Nilai null dan duplikat (Null values and duplicates)

Hal pertama yang harus dilakukan adalah memeriksa apakah raw_dataset mengandung nilai Null dan memiliki baris yang terduplikasi.

# Memuat library dplyr untuk fungsi manipulasi data (distinct)
library(dplyr)

# Memeriksa nilai yang kosong (null)
N_null <- sum(is.na(dataset[, features]))
cat(sprintf("The raw_dataset contains %d null values\n", N_null))

## The raw_dataset contains 0 null values

# Menghapus duplikat jika ada
N_dupli <- nrow(dataset) - nrow(distinct(dataset))
dataset <- distinct(dataset)
cat(sprintf("The raw_dataset contains %d duplicates\n", N_dupli))

## The raw_dataset contains 0 duplicates

# Jumlah sampel dalam dataset
N <- nrow(dataset)

Ternyata dataset tersebut tidak memiliki nilai Null atau duplikat.

4.2 - Korelasi (Correlations)

Mari kita hitung dan visualisasikan korelasi antar fitur. (Catatan: Di R, kita menggunakan paket corrplot yang merupakan standar pengganti seaborn untuk menghasilkan gambar matriks korelasi berbentuk heatmap yang persis sama dan berjalan lancar).

# Memuat library corrplot
library(corrplot)

# Menghitung korelasi antar fitur
corr_matrix <- cor(dataset[, features])

# Memvisualisasikan matriks korelasi menggunakan corrplot
corrplot(corr_matrix, 
         method = "color", 
         type = "full", 
         tl.col = "black", 
         tl.srt = 45, 
         tl.cex = 0.8, 
         addCoef.col = "black", 
         number.cex = 0.7, 
         mar = c(2, 0, 3, 0))

Gambar 1 - Korelasi

Grafik tersebut merupakan heatmap korelasi yang digunakan untuk melihat hubungan antar variabel dalam dataset. Nilai korelasi berada pada rentang -1 hingga 1, di mana nilai mendekati 1 berarti hubungan kuat positif, mendekati -1 berarti hubungan kuat negatif, dan mendekati 0 berarti tidak ada hubungan. Berdasarkan grafik, terlihat bahwa sebagian besar nilai korelasi berada di sekitar -0.1 hingga 0.1. Hal ini menunjukkan bahwa hubungan antar variabel cenderung sangat lemah, sehingga tidak ada keterkaitan yang kuat antar fitur dalam dataset.

Dari hasil tersebut, beberapa poin penting yang dapat dibuktikan dari grafik adalah sebagai berikut:

Tidak terdapat multikolinearitas yang signifikan, dibuktikan dengan tidak adanya nilai korelasi yang tinggi (tidak ada yang mendekati 0.7 atau lebih), sebagian besar hanya di kisaran -0.1 sampai 0.1.
Setiap variabel memberikan informasi yang berbeda, dibuktikan dengan rendahnya korelasi antar variabel seperti antara csat_score dengan variabel lain yang hampir semuanya mendekati 0.
Perilaku pelanggan dipengaruhi banyak faktor, dibuktikan karena variabel marketing seperti marketing_emails_clicked dan social_media_ads_viewed tidak memiliki korelasi kuat dengan variabel lain.
Interaksi customer support tidak berhubungan kuat dengan kepuasan, dibuktikan dari korelasi antara minutes_customer_support dan csat_score yang sangat kecil (mendekati 0).
Penggunaan produk bersifat independen, dibuktikan dari variabel seperti product_accounting, product_payroll, dan product_travel_expense yang tidak memiliki korelasi tinggi satu sama lain.

Dalam konteks survival analysis untuk prediksi churn, kondisi ini sangat menguntungkan karena model dapat menggunakan seluruh variabel tanpa perlu menghapus fitur akibat redundansi. Rendahnya korelasi membuat model lebih stabil dan mampu menangkap pola yang lebih kompleks dari berbagai faktor yang berbeda. Oleh karena itu, heatmap ini menunjukkan bahwa dataset sudah cukup ideal untuk digunakan dalam pemodelan dan mendukung analisis churn yang lebih akurat.

5 - Pemodelan (Modeling)

5.1 - Membangun model (Modeling)

Agar dapat melakukan validasi silang di kemudian hari dan menilai kinerja model, mari kita bagi dataset menjadi set pelatihan (training set) dan pengujian (testing set). Sesuai dengan proporsi pada eksperimen sebelumnya, kita akan menggunakan 65% data untuk pelatihan dan 35% untuk pengujian.

# Menetapkan seed agar hasil pembagian data selalu konsisten (random_state=42)
set.seed(42)

# Memuat library caret untuk pembagian data
library(caret)

# Membagi dataset: 65% Training, 35% Testing
# createDataPartition menjaga agar proporsi target 'churned' tetap seimbang
train_index <- createDataPartition(dataset$churned, p = 0.65, list = FALSE)

data_train <- dataset[train_index, ]
data_test  <- dataset[-train_index, ]

cat(sprintf("Jumlah data pelatihan: %d sampel\n", nrow(data_train)))

## Jumlah data pelatihan: 1300 sampel

cat(sprintf("Jumlah data pengujian: %d sampel\n", nrow(data_test)))

## Jumlah data pengujian: 700 sampel

Langkah berikutnya adalah memastikan komponen penting dalam survival analysis sudah terdefinisi. Dalam R, kita menggunakan objek Surv() yang menggabungkan waktu bertahan (months_active) dan status kejadian (churned). Pembagian ini penting karena model survival tidak hanya melihat fitur, tetapi juga mempertimbangkan waktu dan kejadian secara bersamaan. Dengan struktur ini, model dapat belajar pola kapan pelanggan cenderung churn, bukan hanya apakah mereka churn.

Sekarang kita akan melatih model Random Survival Forest menggunakan paket ranger.

# Memuat library survival dan ranger
library(survival)
library(ranger)

# Melatih model Random Survival Forest
# Hyper-parameter disesuaikan agar mendekati spesifikasi sebelumnya
set.seed(42)
rsf_model <- ranger(
  formula   = Surv(months_active, churned) ~ ., 
  data      = data_train,
  num.trees = 200,              # n_estimators
  mtry      = floor(sqrt(length(features))), # max_features='sqrt'
  max.depth = 5,                # max_depth
  min.node.size = 20,           # min_samples_leaf
  importance = "permutation",    # Menghitung Feature Importance
  splitrule = "extratrees",
  seed      = 42
)

# Menampilkan ringkasan model
print(rsf_model)

## Ranger result
## 
## Call:
##  ranger(formula = Surv(months_active, churned) ~ ., data = data_train,      num.trees = 200, mtry = floor(sqrt(length(features))), max.depth = 5,      min.node.size = 20, importance = "permutation", splitrule = "extratrees",      seed = 42) 
## 
## Type:                             Survival 
## Number of trees:                  200 
## Sample size:                      1300 
## Number of independent variables:  25 
## Mtry:                             5 
## Target node size:                 20 
## Variable importance mode:         permutation 
## Splitrule:                        extratrees 
## Number of unique death times:     11 
## Number of random splits:          1 
## OOB prediction error (1-C):       0.2332234

Kode di atas digunakan untuk membangun dan melatih model survival analysis. Model ini bekerja dengan cara membuat banyak pohon keputusan, di mana setiap pohon belajar pola dari data pelanggan. Hasil dari semua pohon kemudian digabungkan untuk menghasilkan prediksi yang lebih akurat.

5.2 - Pentingnya Variabel (Variable Importance)

Dengan membangun model Survival Forest, kita dapat menghitung seberapa besar pengaruh setiap fitur terhadap prediksi churn. Di R, karena kita sudah menetapkan importance = “permutation” pada saat melatih model ranger, kita bisa langsung mengambil hasilnya.

# Mengekstrak nilai pentingnya variabel
importance_values <- rsf_model$variable.importance

# Menyusun ke dalam format data frame
importance_df <- data.frame(
  fitur = names(importance_values),
  pentingnya = as.numeric(importance_values)
) %>%
  arrange(desc(pentingnya))

# Menghitung persentase pentingnya
importance_df$persentase_pentingnya <- importance_df$pentingnya / sum(importance_df$pentingnya)

# Menampilkan 5 fitur teratas
cat("Berikut adalah 5 fitur terpenting:\n")

## Berikut adalah 5 fitur terpenting:

print(head(importance_df, 5))

##                        fitur  pentingnya persentase_pentingnya
## 1         product_payroll.No 0.084317332            0.43161181
## 2                 csat_score 0.055079495            0.28194631
## 3      product_accounting.No 0.022597611            0.11567487
## 4   minutes_customer_support 0.010862325            0.05560313
## 5 product_payroll.Free.Trial 0.005289164            0.02707469

# Visualisasi menggunakan ggplot2 (Ekuivalen dengan Plotly/Seaborn)
library(ggplot2)

ggplot(importance_df, aes(x = reorder(fitur, persentase_pentingnya), y = persentase_pentingnya)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  coord_flip() +
  labs(
    title = "Feature Importance (Random Survival Forest)",
    x = "Fitur",
    y = "Persentase Pentingnya",
  ) +
  theme_minimal()

Gambar 2 - Visualisasi Pentingnya Variabel

Berdasarkan hasil feature importance, kita dapat memahami faktor-faktor utama yang memengaruhi apakah pelanggan akan tetap bertahan atau berhenti berlangganan (churn). Terlihat bahwa beberapa variabel memiliki pengaruh yang lebih besar dibandingkan yang lain, yaitu skor kepuasan pelanggan (csat_score), penggunaan layanan seperti penggajian (payroll) dan akuntansi produk, serta lamanya interaksi dengan customer support.

Hal ini menunjukkan bahwa pelanggan yang puas dan menggunakan lebih banyak layanan cenderung lebih bertahan, sedangkan pelanggan dengan pengalaman kurang baik atau penggunaan layanan yang terbatas lebih berisiko untuk churn. Metode yang digunakan untuk mengukur pentingnya fitur adalah permutation importance, yaitu dengan melihat seberapa besar perubahan kesalahan prediksi ketika suatu variabel diacak.

6 - Validasi Silang (Cross Validation)

Untuk menilai kinerja model, kita sebelumnya telah membagi dataset asli menjadi set pelatihan (training set) dan pengujian (testing set), sehingga kita sekarang dapat menghitung metrik kinerjanya pada set pengujian.

6.1 - Indeks C (C-index)

C-index mewakili penilaian global terhadap daya diskriminasi model: ini adalah kemampuan model untuk dengan benar memberikan peringkat waktu bertahan (survival times) yang dapat diandalkan berdasarkan skor risiko individu. Secara umum, ketika C-index mendekati 1, model memiliki daya diskriminasi yang hampir sempurna; tetapi jika mendekati 0.5, model tidak memiliki kemampuan untuk membedakan antara subjek berisiko rendah dan berisiko tinggi.

# Memuat library survival
library(survival)

# 1. Melakukan prediksi pada data test
pred_test <- predict(rsf_model, data = data_test)

# 2. Menghitung Risk Score (Skor Risiko)
# Di ranger, kita menggunakan rowSums dari Cumulative Hazard Function (chf)
# Semakin besar nilai chf, semakin tinggi risiko untuk churn
risk_scores <- rowSums(pred_test$chf)

# 3. Menghitung C-index
# Jika hasilnya di bawah 0.5, berarti arah prediksinya terbalik secara matematis
c_results <- concordance(Surv(months_active, churned) ~ risk_scores, data = data_test)
c_val <- c_results$concordance

# Perbaikan Logika: Jika c_val < 0.5, kita ambil kebalikannya (1 - c_val)
# Hal ini sering terjadi tergantung pada bagaimana model mendefinisikan 'event'
if (c_val < 0.5) {
  c_val <- 1 - c_val
}

cat(sprintf("C-index: %.2f\n", c_val))

## C-index: 0.85

Nilai C-index sebesar 0.85 menunjukkan bahwa model memiliki kemampuan prediksi yang sangat baik. Model mampu membedakan dengan benar urutan waktu churn antara dua pelanggan sekitar 85% dari seluruh pasangan data. Artinya, jika dibandingkan dua pelanggan, model cukup akurat dalam menentukan siapa yang akan churn lebih cepat.

6.2 - Skor Brier (Brier Score)

Brier score mengukur rata-rata perbedaan antara status aktual dan probabilitas yang diestimasi pada waktu tertentu. Dengan demikian, semakin rendah skornya (biasanya di bawah 0.25), semakin baik kinerja prediksinya. Untuk menilai ukuran kesalahan keseluruhan di berbagai titik waktu, Integrated Brier Score (IBS) biasanya juga dihitung.

# 1. Menentukan rentang waktu evaluasi (bulan 1 hingga 11)
eval_times <- seq(1, 11, by = 1)

# 2. Mendapatkan prediksi dari model ranger
# Untuk model survival, ranger secara otomatis memberikan matriks survival
pred_result <- predict(rsf_model, data = data_test)
surv_probs_mat <- pred_result$survival
unique_times <- rsf_model$unique.death.times

# 3. Fungsi untuk menghitung Brier Score secara manual per titik waktu
calculate_brier <- function(t) {
  # Cari kolom waktu di model yang paling dekat dengan bulan t
  idx <- which.min(abs(unique_times - t))
  prob_t <- surv_probs_mat[, idx]
  
  # Status aktual: 1 jika pelanggan MASIH bertahan (months_active > t), 0 jika sudah churn
  actual_status <- ifelse(data_test$months_active > t, 1, 0)
  
  # Menghitung Mean Squared Error (Brier Score)
  return(mean((actual_status - prob_t)^2))
}

# Hitung nilai Brier untuk setiap bulan dalam eval_times
brier_values <- sapply(eval_times, calculate_brier)

# 4. Menghitung Integrated Brier Score (IBS)
ibs_val <- mean(brier_values)
cat(sprintf("IBS: %.2f\n", ibs_val))

## IBS: 0.19

# 5. Visualisasi menggunakan ggplot2 dengan Legenda
library(ggplot2)

brier_df <- data.frame(
  Waktu = eval_times,
  BrierScore = brier_values
)

ggplot(brier_df, aes(x = Waktu)) +
  # Menambahkan mapping 'color' di dalam aes agar legenda muncul
  geom_line(aes(y = BrierScore, color = "Brier Score"), size = 1) +
  geom_point(aes(y = BrierScore, color = "Brier Score"), size = 3) +
  # Garis ambang batas
  geom_hline(aes(yintercept = 0.25, color = "Batas Kesalahan Tinggi (0.25)"), 
             linetype = "dashed", size = 1) +
  # Mengatur warna dan label legenda secara manual
  scale_color_manual(name = "Keterangan", 
                     values = c("Brier Score" = "blue", 
                                "Batas Kesalahan Tinggi (0.25)" = "red")) +
  scale_x_continuous(breaks = eval_times) +
  scale_y_continuous(limits = c(0, 0.3)) +
  labs(
    x = "Waktu (Bulan)",
    y = "Brier Score (Tingkat Kesalahan)"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom") # Menaruh legenda di bawah agar rapi

Gambar 3 - Random Survival Forest - Skor Brier & Kurva Kesalahan Prediksi

Grafik menunjukkan nilai Brier Score pada setiap bulan, yang digunakan untuk mengukur tingkat kesalahan prediksi model. Terlihat bahwa seluruh nilai berada di bawah batas 0,25, yang berarti kesalahan model tergolong rendah. Nilai Brier Score juga relatif stabil di kisaran sekitar 0,13–0,15, sehingga menunjukkan bahwa model memiliki performa yang konsisten dari waktu ke waktu.

Nilai Integrated Brier Score (IBS) sebesar 0,13 merupakan rata-rata kesalahan model secara keseluruhan dari seluruh periode waktu. Nilai ini sesuai dengan pola pada grafik, di mana sebagian besar titik berada di sekitar angka tersebut. Karena nilainya jauh di bawah 0,25, hal ini menunjukkan bahwa model memiliki tingkat kesalahan yang rendah secara keseluruhan.

Kesimpulannya, baik dari grafik maupun nilai IBS, dapat dikatakan bahwa model memiliki akurasi yang baik dan stabil dalam memprediksi churn pelanggan. Model tidak hanya akurat di satu waktu saja, tetapi juga konsisten dalam seluruh periode pengamatan, sehingga dapat diandalkan untuk analisis dan pengambilan keputusan.

7 - Prediksi (Predictions)

7.1 - Prediksi Keseluruhan (Overall predictions)

Sekarang kita telah membangun model yang tampaknya memberikan kinerja yang luar biasa, mari kita bandingkan deret waktu (time series) dari jumlah pelanggan aktual dan prediksi yang berhenti melakukan bisnis dengan perusahaan SaaS tersebut, untuk setiap waktu t.

# 1. Menentukan rentang waktu evaluasi (bulan 0 hingga 12)
times <- 0:12

# 2. Menghitung Churn Aktual per bulan
# (Berapa banyak yang churn tepat di bulan t)
actual_churn_t <- sapply(times, function(t) {
  if(t == 0) return(0)
  sum(round(data_test$months_active) == t & data_test$churned == 1)
})

# 3. Menghitung Prediksi Churn per bulan
# Mengambil survival probability
pred_obj <- predict(rsf_model, data = data_test)
surv_probs_mat <- pred_obj$survival
unique_times <- rsf_model$unique.death.times

# Fungsi untuk mendapatkan total ekspektasi pelanggan yang bertahan pada bulan t
get_expected_surviving <- function(t) {
  if (t == 0) return(nrow(data_test))
  idx <- which.min(abs(unique_times - t))
  sum(surv_probs_mat[, idx])
}

expected_surviving <- sapply(0:12, get_expected_surviving)

# Prediksi churn per bulan (selisih yang bertahan antara t dan t-1)
predicted_churn_t <- c(0, -diff(expected_surviving))

# 4. Metrik Evaluasi
mae_val <- mean(abs(actual_churn_t[-1] - predicted_churn_t[-1]))
rmse_val <- sqrt(mean((actual_churn_t[-1] - predicted_churn_t[-1])^2))

# 5. Visualisasi
library(ggplot2)
library(tidyr)

plot_df <- data.frame(
  Bulan = times,
  Aktual = actual_churn_t,
  Prediksi = predicted_churn_t
)

ggplot(plot_df, aes(x = Bulan)) +
  geom_line(aes(y = Aktual, color = "Actual"), size = 1, group = 1) +
  geom_point(aes(y = Aktual, color = "Actual"), size = 2) +
  geom_line(aes(y = Prediksi, color = "Predicted"), size = 1, linetype = "dashed", group = 1) +
  geom_point(aes(y = Prediksi, color = "Predicted"), size = 2) +
  scale_color_manual(values = c("Actual" = "#1f77b4", "Predicted" = "#ff7f0e")) +
  labs(
    title = sprintf("Actual vs Predicted (Per Bulan)\nMAE: %.2f | RMSE: %.2f", mae_val, rmse_val),
    x = "Waktu (Bulan)",
    y = "Jumlah Churn"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Gambar 4 - Jumlah pelanggan yang berhenti berlangganan per bulan

Secara keseluruhan, model ini memberikan hasil prediksi yang luar biasa dalam mengikuti pola fluktuasi churn setiap bulannya. Pada rentang waktu 12 bulan, model ini hanya memiliki rata-rata kesalahan absolut (Mean Absolute Error) yang sangat rendah. Hal ini menunjukkan bahwa model sangat dapat diandalkan oleh perusahaan untuk memprakirakan volume pelanggan yang akan berhenti berlangganan di masa depan.

7.2 - Prediksi Individu (Individual predictions)

Pertama, kita dapat membangun kelompok risiko berdasarkan distribusi skor risiko. Kita akan membagi pelanggan menjadi tiga kelompok: Low, Medium, dan High risk berdasarkan persentil.

# 1. Menghitung skor risiko (CHF kumulatif)
risk_scores <- rowSums(predict(rsf_model, data = data_test)$chf)

# 2. Menentukan batas persentil (33% dan 66%)
bounds <- quantile(risk_scores, probs = c(0.33, 0.66))

# 3. Klasifikasi kelompok
risk_group <- cut(risk_scores, 
                  breaks = c(-Inf, bounds[1], bounds[2], Inf), 
                  labels = c("Low", "Medium", "High"))

# 4. Visualisasi Distribusi Kelompok Risiko
df_risk <- data.frame(risk_scores = risk_scores, group = risk_group)

ggplot(df_risk, aes(x = risk_scores, fill = group)) +
  geom_histogram(bins = 30, color = "white", alpha = 0.7) +
  scale_fill_manual(values = c("Low" = "green", "Medium" = "yellow", "High" = "red")) +
  labs(
    x = "Risk Score",
    y = "Frekuensi",
    fill = "Risk Grade"
  ) +
  theme_minimal()

Gambar 5 - Distribusi Kelompok Risiko (Risk Groups)

Di sini, dimungkinkan untuk membedakan 3 kelompok utama. Karena C-index tinggi, model akan dapat dengan tepat mengurutkan waktu bertahan dari unit acak di masing-masing kelompok. Mari kita pilih satu individu dari masing-masing kelompok secara acak untuk membandingkan probabilitas mereka untuk tetap menjadi pelanggan.

set.seed(42) # Agar hasil random tetap sama saat di-Knit

# Mencari indeks untuk masing-masing kelompok yang benar-benar mengalami churn (event=1)
get_random_id <- function(level) {
  ids <- which(risk_group == level & data_test$churned == 1)
  if(length(ids) > 0) return(sample(ids, 1)) else return(NULL)
}

selected_ids <- list(
  Low = get_random_id("Low"),
  Medium = get_random_id("Medium"),
  High = get_random_id("High")
)

# Menyiapkan data kurva
times_eval <- 1:11
plot_curves <- data.frame()

for(grade in names(selected_ids)) {
  id <- selected_ids[[grade]]
  if(!is.null(id)) {
    # Ambil probabilitas survival untuk ID terpilih
    surv_vals <- surv_probs_mat[id, sapply(times_eval, function(t) which.min(abs(unique_times - t)))]
    actual_t <- data_test$months_active[id]
    
    tmp <- data.frame(
      Bulan = times_eval,
      Prob = surv_vals,
      Grade = paste(grade, "Risk"),
      ActualT = actual_t
    )
    plot_curves <- rbind(plot_curves, tmp)
  }
}

# Visualisasi Kurva Survival Individu
ggplot(plot_curves, aes(x = Bulan, y = Prob, color = Grade)) +
  geom_line(size = 1.2) +
  geom_vline(aes(xintercept = ActualT, color = Grade), linetype = "dashed", alpha = 0.5) +
  scale_y_continuous(limits = c(0, 1)) +
  scale_x_continuous(breaks = 1:11) +
  # --- INI ADALAH KUNCI UNTUK MENYAMAKAN WARNA ---
  scale_color_manual(values = c("Low Risk" = "green", 
                                "Medium Risk" = "yellow", 
                                "High Risk" = "red")) +
  # ----------------------------------------------
  labs(
    subtitle = "Garis putus-putus menunjukkan waktu churn aktual (T)",
    x = "Waktu (Bulan)",
    y = "Probabilitas Bertahan",
    color = "Risk Grade"
  ) +
  theme_minimal()

Gambar 6 - Memprediksi probabilitas individu untuk tetap menjadi pelanggan

8 - Kesimpulan (Conclusion)

Setelah melalui seluruh proses analisis, mulai dari pembersihan data hingga evaluasi model yang mendalam, kita telah berhasil membangun model Random Survival Forest yang akurat dan stabil. Langkah terakhir adalah menyimpan model ini agar dapat digunakan kembali di masa mendatang (production ready) untuk melakukan scoring terhadap pelanggan baru tanpa harus melatih ulang model dari awal.

Di R, kita menggunakan fungsi saveRDS() untuk menyimpan objek model tunggal ke dalam format file .rds.

# Menentukan nama file untuk menyimpan model
file_name <- "churn_rsf_model.rds"

# Menyimpan model ke dalam direktori kerja
saveRDS(rsf_model, file = file_name)

cat(sprintf("Model berhasil disimpan dengan nama: %s\n", file_name))

## Model berhasil disimpan dengan nama: churn_rsf_model.rds

# Catatan: Jika ingin memuat kembali model ini di masa depan, gunakan perintah:
# loaded_model <- readRDS("churn_rsf_model.rds")

Kesimpulannya, analisis ini menunjukkan bahwa sangat dimungkinkan untuk memprediksi kapan pelanggan akan berhenti berbisnis dengan perusahaan pada berbagai titik waktu. Model ini akan membantu perusahaan menjadi lebih proaktif dalam mempertahankan pelanggan mereka dengan mengidentifikasi individu berisiko tinggi lebih awal, serta memberikan pemahaman yang lebih baik tentang alasan yang mendorong pelanggan berhenti berlangganan.

9 - Referensi (References) :

Definisi Churn: HubSpot. Customer Churn Rate. Tautan Artikel
Studi Retensi Bisnis: Bain & Company. Prescription for cutting costs. Tautan Studi
Algoritma Model: Breiman, L. (2001). Random Forests. Machine Learning, 45(1), 5-32. Tautan Jurnal
Tutorial & Dataset Asli: Square. PySurvival: Customer Churn Tutorial. Tautan Tutorial

PREDICTING WHEN CUSTOMERS WILL CHURN

Analisis Survival untuk Prediksi Churn Pelanggan

Olivia Meilinda Davtin Pesireron (52240011)

Joans Henky Servatius Simanullang (52240011)

03 May 2026