Import Data

Data yang digunakan merupakan data profil pelanggan dari sebuah perusahaan telekomunikasi yang diperoleh dari Kaggle. Dataset tersebut berisikan data untuk 7043 pelanggan yang meliputi demografis pelanggan, informasi pembayaran akun, serta produk layanan yang didaftarkan oleh tiap pelanggan. Dari informasi tersebut, kita ingin memprediksi apakah seorang pelanggan akan Churn atau tidak.

customer <- read.csv("data_input/Telco-Customer-Churn.csv", stringsAsFactors = T)
head(customer)

Berikut ini merupakan deskripsi untuk setiap variabel:

• CustomerID: Customer ID
• Gender: Gender pelanggan yaitu Female dan Male
• SeniorCitizen: Apakah pelanggan merupakan senio citizen (0: No, 1: Yes)
• Partner: Apakah pelanggan memiliki partner atau tidak (Yes, No)
• Dependents: Apakah pelanggan memiliki tanggungan atau tidak (Yes, No)
• Tenure: Jumlah bulan dalam menggunakan produk perusahaan
• MultipleLines: Apakah pelanggan memiliki banyak saluran atau tidak (Yes, No, No phone service)
• OnlineSecurity: Apakah pelanggan memiliki keamanan online atau tidak
• OnlineBackup: Apakah pelanggan memiliki cadangan online atau tidak
• DeviceProtection: Apakah pelanggan memiliki perlindungan perangkat atau tidak
• TechSupport: Apakah pelanggan memiliki dukungan teknis atau tidak
• StreamingTV: Apakah pelanggan berlangganan TV streaming atau tidak
• StreamingMovies: Apakah pelanggan berlangganan movies streaming atau tidak
• Contract: Ketentuan kontrak berlangganan (Month-to-month, One year, Two year)
• PaperlessBilling: Apakah pelanggan memiliki tagihan tanpa kertas atau tidak (Yes, No)
• PaymentMethod: Metode pembayaran (Electronic check, Mailed check, Bank transfer (automatic), Credit card (automatic))
• MonthlyCharges: Jumlah pembayaran yang dilakukan setiap bulan
• TotalCharges: Jumlah total yang dibebankan oleh pelanggan
• Churn: Apakah pelanggan Churn atau tidak (Yes or No)

Data Cleansing

Sebelum masuk ke tahap modeling, mari kita membersihkan datanya terlebih dahulu.

Pertama, cek kelengkapan data, dari tahap ini kita akan memperoleh informasi apakah data kita sudah lengkap.

colSums(is.na(customer))

#>       customerID           gender    SeniorCitizen          Partner 
#>                0                0                0                0 
#>       Dependents           tenure     PhoneService    MultipleLines 
#>                0                0                0                0 
#>  InternetService   OnlineSecurity     OnlineBackup DeviceProtection 
#>                0                0                0                0 
#>      TechSupport      StreamingTV  StreamingMovies         Contract 
#>                0                0                0                0 
#> PaperlessBilling    PaymentMethod   MonthlyCharges     TotalCharges 
#>                0                0                0               11 
#>            Churn 
#>                0

Dari 7043 observasi ternyata terdapat sebanyak 11 observasi pada kolom TotalCharges yang merupakan missing values (NA). Dikarenakan jumlah NA yang cukup sedikit, kita dapat membuang observasi tersebut.

Kedua, kita perlu buang variabel yang tidak relevan dengan pemodelan, yaitu CustomerID.

Ketiga, kita menyesuaikan tipe data kolom SeniorCitizen yang sebelumnya numerik menjadi kategorik.

customer <- customer %>% 
            select(-customerID) %>% 
            na.omit() %>% 
            mutate(SeniorCitizen = as.factor(SeniorCitizen)) 

Exploratory Data Analysis

Selanjutnya mari lakukan eksplorasi data baik untuk kolom kategorik maupun numerik.

Untuk mengetahui proporsi kelas pada setiap variabel kategori, kita dapat menggunakan fungsi inspect_cat dari package inspectdf seperti berikut:

customer %>% inspect_cat() %>% show_plot()

Dari visualisasi di atas dapat diketahui proporsi kelas untuk variabel target Churn lebih banyak di kategori No dibandingkan Yes. Lalu, untuk proporsi variabel lainnya mayoritas seimbang.

Berikutnya kita dapat eksplorasi persebaran untuk variabel data numerik dengan fungsi inspect_num dari package inspectdf seperti berikut:

customer %>% inspect_num() %>% show_plot()

Dari visualisasi di atas dapat disimpulkan bahwa persebaran data numerik cukup beragam untuk setiap variabelnya.

Train-Test Splitting

Setelah kita melakukan data cleansing dan eksplorasi data, tahap berikutnya adalah train-test splitting yaitu membagi data menjadi data train dan test dengan proporsi 80:20. Data train digunakan untuk membuat model sedangkan data test digunakan untuk mengevaluasi performa model.

set.seed(100)
idx <- initial_split(data = customer,
                     prop = 0.8,
                     strata = "Churn")
data_train <- training(idx)
data_test <- testing(idx)

Modeling

Selanjutnya kita akan melakukan modeling menggunakan algoritma Random Forest (package caret) dengan menentukan banyaknya cross validation, repetisi, serta mencantumkan nama target variabel dan juga prediktor yang digunakan dari data train.

set.seed(100)
ctrl <- trainControl(method = "repeatedcv",
                     number = 5,
                     repeats = 3)
model_forest <- train(Churn ~ .,
                      data = data_train,
                      method = "rf",
                      trControl = ctrl)
# saveRDS(model_forest, "assets/model_forest.rds")

Chunk di atas membutuhkan waktu yang cukup lama untuk dieksekusi. Untuk mempersingkat waktu, mari load model yang sebelumnya sudah disimpan ke dalam bentuk file RDS.

model_forest <- readRDS("assets/model_forest.rds")
model_forest

#> Random Forest 
#> 
#> 5627 samples
#>   19 predictor
#>    2 classes: 'No', 'Yes' 
#> 
#> No pre-processing
#> Resampling: Cross-Validated (5 fold, repeated 3 times) 
#> Summary of sample sizes: 4501, 4502, 4501, 4502, 4502, 4501, ... 
#> Resampling results across tuning parameters:
#> 
#>   mtry  Accuracy   Kappa    
#>    2    0.7837817  0.3252122
#>   16    0.7750746  0.3779712
#>   30    0.7731203  0.3727503
#> 
#> Accuracy was used to select the optimal model using the largest value.
#> The final value used for the model was mtry = 2.

Untuk saat ini, kita memperoleh model Random Forest dengan tingkat akurasi di data train sebesar 78,38% dengan nilai mtry optimum sebanyak 2.

Selanjutnya kita akan melakukan tuning model dengan melakukan upsampling, yang artinya kita akan menyetarakan proporsi target variabel menjadi sama besar.

data_train_up <- upSample(x = data_train[, -20],
                          y = data_train$Churn,
                          yname = "Churn")

# cek proporsi
prop.table(table(data_train_up$Churn))

#> 
#>  No Yes 
#> 0.5 0.5

Dari data yang sudah dilakukan upsampling, kita akan membuat ulang model Random Forestnya.

set.seed(100)
ctrl <- trainControl(method = "repeatedcv",
                     number = 5,
                     repeats = 3)
model_forest_up <- train(Churn ~ .,
                         data = data_train_up,
                         method = "rf",
                         trControl = ctrl)
# saveRDS(model_forest_up, "assets/model_forest_up.rds")

Untuk mempersingkat waktu, mari load model yang sebelumnya sudah disimpan ke dalam bentuk file RDS.

model_forest_up <- readRDS("assets/model_forest_up.rds")
model_forest_up

#> Random Forest 
#> 
#> 8262 samples
#>   19 predictor
#>    2 classes: 'No', 'Yes' 
#> 
#> No pre-processing
#> Resampling: Cross-Validated (5 fold, repeated 3 times) 
#> Summary of sample sizes: 6609, 6610, 6609, 6610, 6610, 6610, ... 
#> Resampling results across tuning parameters:
#> 
#>   mtry  Accuracy   Kappa    
#>    2    0.7760017  0.5520022
#>   16    0.8911472  0.7822945
#>   30    0.8875167  0.7750336
#> 
#> Accuracy was used to select the optimal model using the largest value.
#> The final value used for the model was mtry = 16.

Setelah dilakukan upsampling, terlihat bahwa nilai akurasi di data train meningkat menjadi 89,11% dengan nilai mtry optimum sebanyak 16.

Model Evaluation

Terakhir, mari kita uji model random forest yang telah kita buat ke data test. Pada kasus ini, kita ingin memperoleh nilai recall atau sensitivitas yang sebesar mungkin agar model kita dapat mendeteksi pelanggan yang sebenarnya Churn sebanyak-banyaknya.

pred <- predict(model_forest_up, newdata = data_test, type = "prob")
pred$result <- as.factor(ifelse(pred$Yes > 0.45, "Yes", "No"))
confusionMatrix(pred$result, data_test$Churn, positive = "Yes")

#> Confusion Matrix and Statistics
#> 
#>           Reference
#> Prediction  No Yes
#>        No  849 109
#>        Yes 183 264
#>                                         
#>                Accuracy : 0.7922        
#>                  95% CI : (0.77, 0.8131)
#>     No Information Rate : 0.7345        
#>     P-Value [Acc > NIR] : 0.00000031    
#>                                         
#>                   Kappa : 0.4989        
#>                                         
#>  Mcnemar's Test P-Value : 0.00001937    
#>                                         
#>             Sensitivity : 0.7078        
#>             Specificity : 0.8227        
#>          Pos Pred Value : 0.5906        
#>          Neg Pred Value : 0.8862        
#>              Prevalence : 0.2655        
#>          Detection Rate : 0.1879        
#>    Detection Prevalence : 0.3181        
#>       Balanced Accuracy : 0.7652        
#>                                         
#>        'Positive' Class : Yes           
#> 

Dengan menggunakan threshold 0.45, diperoleh recall sebesar 70,78% dengan akurasi sebesar 79,22%.

Selain menggunakan confusion matrix, kita dapat membentuk kurva ROC beserta nilai AUC dengan menggunakan package ROCR sebagai berikut:

pred_prob <- predict(object = model_forest_up, newdata = data_test, type = "prob")
pred <-  prediction(pred_prob[,2], labels = data_test$Churn)
perf <- performance(prediction.obj = pred, measure = "tpr", x.measure = "fpr")
plot(perf)

auc <- performance(pred, measure = "auc")
auc@y.values[[1]]

#> [1] 0.8513259

Nilai AUC di atas menyatakan bahwa performa model kita sebesar 85,13% dalam memisahkan distribusi kelas positif Churn dengan negatif pada data test.