Package

Beberapa package yang akan digunakan dalam laporan ini antara lain sebagai berikut:

library(readxl) # untuk membaca file excel
library(sjmisc) # untuk deskripsi statistik data
library(grid) # untuk menyusun >1 plot dalam 1 frame
library(gridExtra) # untuk menyusun >1 plot dalam 1 frame
library(ggplot2) # membuat plot
library(lattice) # membuat barchart
library(smotefamily) # melakukan SMOTE
library(randomForest) # membentuk model Random Forest
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
## 
##     margin
## The following object is masked from 'package:gridExtra':
## 
##     combine
library(randomForestExplainer) # memperoleh visualisasi hasil pemodelan RF
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(caret) # membuat matriks konfusi
library(pROC) # membuat ROC curve
## Type 'citation("pROC")' for a citation.
## 
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

Import Data

Sebuah operator telekomunikasi seluler akan menjajagi kerja sama dengan sebuah bank untuk meningkatkan angka recharge pelanggan pra-bayar melalui channel perbankan (ATM, internet banking, dsb). Mereka ingin memperoleh model yang dapat digunakan untuk memprediksi pelanggan mana saja yang potensial untuk diajak menggunakan channel tersebut.

data <- read_excel('Tugas_STA581.xlsx')
data <- data[-1] # menghapus kolom ID dalam data

Digunakan package readxl untuk membaca data excel ke dalam bentuk dataframe di R.

head(data)
## # A tibble: 6 × 12
##      X1    X2    X3    X4    X5    X6    X7    X8    X9   X10   X11     Y
##   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1   257 13643    27 0.780 0.963  5217 0.769   768 30006   6.3   4.7     1
## 2    24 13897    15 0.739 0.906 18025 0.291  2368 39954   7.3   5.3     1
## 3     1   760     1 0.5   1      3275 0.962  2542 43544   1.3   1.3     0
## 4   627 41160    97 0.707 0.998  2864 0.723     0     0   9.3   5       0
## 5    42 19002    38 0.750 0.964  2794 0.748     0     0   4     1       0
## 6    65 38253    43 0.779 1      1122 1         5  1719  17    14.7     1

Sebelumnya dilakukan pembagian data ke dalam data latih dan data uji yang selanjutnya data yang sama ini akan dipakai menjadi data latih dan data uji di semua metode.

Eksplorasi Data

Statistik Deskriptif

eksp_data <- data
eksp_data$Y <- as.factor(eksp_data$Y)
as.data.frame(descr(eksp_data))
##    var        type label     n NA.prc         mean           sd           se
## 1   X1     numeric    X1 50000      0 1.106676e+02 1.806347e+02 8.078230e-01
## 4   X2     numeric    X2 50000      0 2.790754e+04 2.736161e+04 1.223649e+02
## 5   X3     numeric    X3 50000      0 4.023486e+01 5.120954e+01 2.290160e-01
## 6   X4     numeric    X4 50000      0 7.342076e-01 1.395345e-01 6.240174e-04
## 7   X5     numeric    X5 50000      0 8.824384e-01 1.950748e-01 8.724011e-04
## 8   X6     numeric    X6 50000      0 6.745711e+03 7.915107e+03 3.539743e+01
## 9   X7     numeric    X7 50000      0 7.718046e-01 2.549293e-01 1.140078e-03
## 10  X8     numeric    X8 50000      0 9.468719e+02 2.260855e+03 1.011085e+01
## 11  X9     numeric    X9 50000      0 1.614417e+04 2.682608e+04 1.199699e+02
## 2  X10     numeric   X10 50000      0 7.129130e+00 4.841823e+00 2.165329e-02
## 3  X11     numeric   X11 50000      0 3.676492e+00 2.414185e+00 1.079656e-02
## 12   Y categorical     Y 50000      0 1.867200e-01 3.896905e-01 1.742749e-03
##              md      trimmed             range          iqr      skew
## 1  4.700000e+01 7.218118e+01     4130 (0-4130) 1.170000e+02  4.657092
## 4  2.209200e+04 2.385257e+04 638172 (0-638172) 2.600850e+04  3.929080
## 5  2.400000e+01 3.022100e+01     1007 (0-1007) 4.000000e+01  3.717217
## 6  7.498324e-01 7.485643e-01           1 (0-1) 1.035349e-01 -2.745926
## 7  9.702921e-01 9.286175e-01           1 (0-1) 1.437488e-01 -2.588062
## 8  3.975000e+03 5.284058e+03 244444 (0-244444) 7.839000e+03  3.215393
## 9  8.626017e-01 8.155916e-01           1 (0-1) 3.416490e-01 -1.314261
## 10 7.300000e+01 4.329732e+02   59611 (0-59611) 9.250000e+02  5.749438
## 11 2.514500e+03 1.042832e+04 504086 (0-504086) 2.443125e+04  3.386872
## 2  6.000000e+00 6.523237e+00     30.7 (0-30.7) 6.400000e+00  1.207519
## 3  3.000000e+00 3.349972e+00     35.7 (0-35.7) 2.700000e+00  1.991754
## 12 0.000000e+00 1.084000e-01           1 (0-1) 0.000000e+00  1.607904

Keterangan:

  • X1 : rata-rata durasi melakukan panggilan per bulan (menit)
  • X2 : rata-rata penghasilan dari transaksi panggilan per bulan (Rp)
  • X3 : rata-rata frekuensi melakukan panggilan per bulan
  • X4 : persentase panggilan di jam kerja
  • X5 : persentase panggilan sesama operator
  • X6 : rata-rata penghasilan dari transaksi SMS per bulan (Rp)
  • X7 : persentase transaksi SMS sesama operator
  • X8 : rata-rata pemakaian data per bulan (MB)
  • X9 : rata-rata penghasilan dari transaksi pemakaian data per bulan (Rp)
  • X10 : rata-rata frekuensi melakukan recharge pulsa per bulan
  • X11 : rata-rata banyaknya wilayah yang dikunjungi per bulan
  • Y : status kepotensialan pelanggan untuk diajak menggunakan channel perbankan. (1 = potensial, 0 = tidak potensial)

Proporsi Pelanggan Potensial dan Tidak Potensial

ggplot(eksp_data, aes(x = factor(1), fill = Y)) + geom_bar(width = 1) + scale_fill_brewer(palette="Set2") + coord_polar("y") + stat_count(geom = "text", aes(label = stat(count)/50000), position=position_stack(vjust=0.5), colour="white")
## Warning: `stat(count)` was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Dari plot di atas diperoleh bahwa dari 50000 pelanggan yang diamati, sebanyak 18.672% pelanggan merupakan pelanggan potensial dan sebanyak 81.328% sisanya bukan merupakan pelanggan yang potensial untuk diajak menggunakan channel perbankan.

Visualisasi Persebaran Masing-masing Peubah prediktor

x1 <- ggplot(eksp_data, aes(X1)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3F2B44") +
  geom_density(color = "#000000", fill = "#E3ADB5", alpha = 0.6) + ggtitle('Rata-rata Durasi Melakukan Panggilan per Bulan (menit)')
x2 <- ggplot(eksp_data, aes(X2)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#231F20") +
  geom_density(color = "#000000", fill = "#95DFE3", alpha = 0.6) + ggtitle('Rata-rata Penghasilan dari Transaksi Panggilan per Bulan (Rp)')
x3 <- ggplot(eksp_data, aes(X3)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3F2B44") +
  geom_density(color = "#000000", fill = "#E3ADB5", alpha = 0.6) + ggtitle('Rata-rata Frekuensi Melakukan Panggilan per Bulan')
x4 <- ggplot(eksp_data, aes(X4)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3F2B44") +
  geom_density(color = "#000000", fill = "#E3ADB5", alpha = 0.6) + ggtitle('Persentase Panggilan di Jam Kerja')
x5 <- ggplot(eksp_data, aes(X5)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#454727") +
  geom_density(color = "#000000", fill = "#95B8E3", alpha = 0.6) +ggtitle('Persentase Panggilan Sesama Operator')
x6 <- ggplot(eksp_data, aes(X6)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#231F20") +
  geom_density(color = "#000000", fill = "#95DFE3", alpha = 0.6) +ggtitle('Rata-rata Penghasilan dari Transaksi SMS per Bulan (Rp)')
x7 <- ggplot(eksp_data, aes(X7)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#454727") +
  geom_density(color = "#000000", fill = "#95B8E3", alpha = 0.6) + ggtitle("Persentase Transaksi SMS Sesama Operator")
x8 <- ggplot(eksp_data, aes(X8)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3E2A1F") +
  geom_density(color = "#000000", fill = "#A99887", alpha = 0.6) + ggtitle("Rata-rata Pemakaian Data per Bulan (MB)")
x9 <- ggplot(eksp_data, aes(X9)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#231F20") +
  geom_density(color = "#000000", fill = "#95DFE3", alpha = 0.6) + ggtitle("Rata-rata Penghasilan dari Transaksi Pemakaian Data per Bulan (Rp)")
x10 <- ggplot(eksp_data, aes(X10)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3E2A1F") +
  geom_density(color = "#000000", fill = "#A99887", alpha = 0.6) + ggtitle("Rata-rata Frekuensi Melakukan Recharge Pulsa per Bulan")
x11 <- ggplot(eksp_data, aes(X11)) +
  geom_histogram(aes(y = ..density..), color = "#000000", fill = "#3E2A1F") +
  geom_density(color = "#000000", fill = "#A99887", alpha = 0.6) + ggtitle("Rata-rata Banyaknya Wilayah yang Dikunjungi per Bulan")
blank<-rectGrob(gp=gpar(col="white"))
figure1 <- grid.arrange(blank, x1, blank,x3, blank, x4, blank, ncol=1, heights = c(1/40, 3/10, 1/40, 3/10, 1/40, 3/10, 1/40))
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

figure2 <- grid.arrange(blank, x2, blank,x6, blank, x9, blank, ncol=1, heights = c(1/40, 3/10, 1/40, 3/10, 1/40, 3/10, 1/40))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

figure3 <- grid.arrange(blank, x8, blank,x10, blank, x11, blank, ncol=1, heights = c(1/40, 3/10, 1/40, 3/10, 1/40, 3/10, 1/40))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

figure4 <- grid.arrange(blank, x5, blank,x7, blank, ncol=1, heights = c(1/30, 4.5/10, 1/30, 4.5/10, 1/30))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Menangani Imbalance dalam Data

Metode Synthetic Minority Over-sampling Technique (SMOTE) merupakan metode yang populer diterapkan dalam rangka menangani ketidak seimbangan kelas. Teknik ini mensintesis sampel baru dari kelas minoritas untuk menyeimbangkan dataset dengan cara sampling ulang sampel kelas minoritas.

# data tanpa SMOTE
set.seed(2022)
sample <- sample(c(TRUE, FALSE), nrow(data), replace=TRUE, prob=c(0.8,0.2))
train  <- data[sample, ]
test   <- data[!sample, ]

Dalam hal ini, ketidakseimbangan kelas dalam data train dilihat dalam grafik berikut:

barchart(as.factor(train$Y), col='maroon')

Terlihat bahwa proporsi pelanggan yang tidak potensial jauh lebih besar daripada pelanggan potensial. Karena data latih digunakan untuk membentuk model prediksi padahal datanya tidak seimbang, maka perlu dilakukan SMOTE.

smote_data <- SMOTE(data[,-12], data$Y)
newdata <- smote_data$data
set.seed(2022)
newdata$class <- as.numeric(newdata$class)
sample <- sample(c(TRUE, FALSE), nrow(newdata), replace=TRUE, prob=c(0.8,0.2))
newtrain  <- newdata[sample, ]
newtest   <- newdata[!sample, ]
barchart(as.factor(newtrain$class), col='navy')

Terlihat bahwa data latih sudah balance dan siap digunakan dalam klasifikasi.

Klasifikasi dengan Random Forest Classifier

Random Forest adalah kombinasi dari masing – masing tree yang baik kemudian dikombinasikan ke dalam satu model. Random Forest bergantung pada sebuah nilai vector random dengan distribusi yang sama pada semua pohon yang masing masing decision tree memiliki kedalaman yang maksimal. Random Forest adalah classifier yang terdiri dari classifier yang berbentuk pohon {\(h(x, \theta_k), k=1. ...\)} dimana \(\theta_k\) adalah random vector yang diditribusikan secara independen dan masing masing tree pada sebuah unit kan memilih class yang paling popular pada input x.

rf_data <- newdata
rf_train <- newtrain
rf_test <- test

Model RF dibentuk dengan menggunakan fungsi randomForest() dan kemudian diperoleh hasil akurasi sebagai berikut:

set.seed(2022)
rf_model <- randomForest(data=rf_train,
               as.factor(class)~.,
               ntree=500)
rf_model
## 
## Call:
##  randomForest(formula = as.factor(class) ~ ., data = rf_train,      ntree = 500) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 14.17%
## Confusion matrix:
##       0     1 class.error
## 0 28377  4159   0.1278276
## 1  4678 25158   0.1567905
plot(rf_model)

Plot di atas merupakan error yang dihasilkan dari banyak pohon yang berbeda-beda. Model yang dibentuk sebelumnya memiliki 500 pohon dengan error rate terkecil.

varImpPlot(rf_model)

Plot di atas menjelaskan kontribusi masing-masing peubah di mana X11 menjadi peubah prediktor yang paling penting dalam memprediksi Y.

rf_pred <- predict(rf_model, rf_test, type="prob")
rf_predict <- ifelse(rf_pred[,2] > 0.5, 1, 0)
rf_conf <- confusionMatrix(as.factor(rf_predict), as.factor(rf_test$Y), positive = "1", mode='everything')
rf_conf
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 7983  170
##          1  203 1656
##                                           
##                Accuracy : 0.9627          
##                  95% CI : (0.9588, 0.9664)
##     No Information Rate : 0.8176          
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.876           
##                                           
##  Mcnemar's Test P-Value : 0.09754         
##                                           
##             Sensitivity : 0.9069          
##             Specificity : 0.9752          
##          Pos Pred Value : 0.8908          
##          Neg Pred Value : 0.9791          
##               Precision : 0.8908          
##                  Recall : 0.9069          
##                      F1 : 0.8988          
##              Prevalence : 0.1824          
##          Detection Rate : 0.1654          
##    Detection Prevalence : 0.1857          
##       Balanced Accuracy : 0.9411          
##                                           
##        'Positive' Class : 1               
##

Diperoleh nilai-nilai akurasi sebagai berikut:

rf_roc <- roc(rf_test$Y, rf_pred[,2], direction = "<", plot = T)
## Setting levels: control = 0, case = 1

rf_roc$auc
## Area under the curve: 0.9854

Diperoleh nilai AUC sebesar 98.54% yang artinya model mampu mengklasifikasikan 98.54% kelas dengan benar. Sehingga, berdasarkan klasifikasi dengan random forest diperoleh ukuran-ukuran akurasi sebagai berikut:

rf_df <- data.frame(Metode = 'Rand Forest',
                     Akurasi = rf_conf$overall[1],
                     Sensitivitas = rf_conf$byClass[1],
                     Spesifisitas = rf_conf$byClass[2],
                     AUC = rf_roc$auc)
rf_df
##               Metode   Akurasi Sensitivitas Spesifisitas       AUC
## Accuracy Rand Forest 0.9627447    0.9069003    0.9752016 0.9854496

Hasil Prediksi

rf_result <- data.frame(actual = rf_test$Y,
                        predicted = rf_predict)
ggplot(rf_result, aes(x = factor(1), fill = as.factor(predicted))) + geom_bar(width = 1) + scale_fill_brewer(palette="Set1") + coord_polar("y") + stat_count(geom = "text", aes(label = stat(count)/10012), position=position_stack(vjust=0.5), colour="white")

Berdasarkan hasil prediksi yang terbentuk, diperoleh bahwa sebanyak 18.57% pelanggan diprediksi akan menjadi pelanggan potensial dan sisanya sebanyak 81.43 menjadi pelanggan tidak potensial. Selanjutnya akan dilihat bagaimana hubungan masing-masing peubah prediktor dengan peubah respon.

plot_predict_interaction(rf_model, rf_test, 'X2', 'X10')

Pelanggan yang menghasilkan transaksi panggilan kecil dan sering melakukan recharge pulsa per bulan cenderung diklasifikasikan dalam kelas pelanggan tidak potensial. Sebaliknya, pelanggan dengan nominal transaksi panggilan besar (> Rp50000) dengan frekuensi recharge pulsa per bulan < 10, cenderung diklasifikasikan ke dalam pelanggan potensial.

plot_predict_interaction(rf_model, rf_test, 'X4', 'X6')

Pelanggan yang melakukan panggilan di jam kerja secara intens (presentase > 75%) dengan nominal transaksi SMS besar ( > Rp50000) cenderung diklasifikasikan ke dalam pelanggan potensial.

plot_predict_interaction(rf_model, rf_test, 'X6', 'X10')

Pelanggan yang memiliki nominal transaksi SMS bulanan kecil (< Rp50000) dengan rata-rata frekuensi recharge pulsa per bulan banyak (> 10 kali) cenderung diklasifikasikan ke dalam pelanggan tidak potensial.

plot_predict_interaction(rf_model, rf_test, 'X8', 'X9')

Pelanggan yang memiliki rata-rata pemakaian data kecil dalam sebulan serta rata-rata mengeluarkan nominal yang kecil untuk data bulanan cenderung diklasifikasikan ke dalam kelompok pelanggan tidak potensial.

plot_predict_interaction(rf_model, rf_test, 'X10', 'X11')

Pelanggan yang tidak mengunjungi banyak wilayah selama sebulan (< 5 wilayah) dengan rata-rata recharge pulsa di bawah 10 kali selama sebulan cenderung diklasifikasikan ke dalam pelanggan potensial.