Model Klasifikasi Berbasis Pohon

Pengantar

Materi ini dimaksudkan untuk memberikan ilustrasi pendekatan pemodelan prediktif berbasis Pohon Klasifikasi. Ada tiga pendekatan yang akan dicover dalam materi ini yaitu:

• pohon klasifikasi

• random forest untuk klasifikasi

• gradient boosting untuk pohon klasifikasi

Sebagaimana telah diketahui, pohon klasifikasi merupakan salah satu metode yang dapat menghasilkan model prediktif dengan variabel respon berupa variabel kategorik (kelas). Hal ini yang membedakan dengan istilah pohon regresi yang melibatkan variabel respon numerik/kontinu.

Ilustrasi yang akan digunakan pada materi ini adalah data yang digunakan untuk memprediksi apakah suatu email terkategori sebagai spam atau bukan. Dataset yang digunakan bersumber pada Hopkins,Mark, Reeber,Erik, Forman,George, and Suermondt,Jaap. (1999). Spambase. UCI Machine Learning Repository. https://doi.org/10.24432/C53G6X.

Dataset telah didownload dan disimpan sebagai suatu file CSV. Perintah di bawah ini berguna untuk mengimport data pada file CSV menjadi R-dataframe.

spam_data <- read.csv("spambase.csv",
                      stringsAsFactors =TRUE )
dim(spam_data)

[1] 4601   58

colnames(spam_data)

 [1] "word_freq_make"             "word_freq_address"         
 [3] "word_freq_all"              "word_freq_3d"              
 [5] "word_freq_our"              "word_freq_over"            
 [7] "word_freq_remove"           "word_freq_internet"        
 [9] "word_freq_order"            "word_freq_mail"            
[11] "word_freq_receive"          "word_freq_will"            
[13] "word_freq_people"           "word_freq_report"          
[15] "word_freq_addresses"        "word_freq_free"            
[17] "word_freq_business"         "word_freq_email"           
[19] "word_freq_you"              "word_freq_credit"          
[21] "word_freq_your"             "word_freq_font"            
[23] "word_freq_000"              "word_freq_money"           
[25] "word_freq_hp"               "word_freq_hpl"             
[27] "word_freq_george"           "word_freq_650"             
[29] "word_freq_lab"              "word_freq_labs"            
[31] "word_freq_telnet"           "word_freq_857"             
[33] "word_freq_data"             "word_freq_415"             
[35] "word_freq_85"               "word_freq_technology"      
[37] "word_freq_1999"             "word_freq_parts"           
[39] "word_freq_pm"               "word_freq_direct"          
[41] "word_freq_cs"               "word_freq_meeting"         
[43] "word_freq_original"         "word_freq_project"         
[45] "word_freq_re"               "word_freq_edu"             
[47] "word_freq_table"            "word_freq_conference"      
[49] "char_freq_."                "char_freq_..1"             
[51] "char_freq_..2"              "char_freq_..3"             
[53] "char_freq_..4"              "char_freq_..5"             
[55] "capital_run_length_average" "capital_run_length_longest"
[57] "capital_run_length_total"   "spam"                      

Dataset ini memuat 4601 amatan dimana satuan amatannya adalah sebuah email. Setiap amatan telah teridentifikasi apakah ini adalah email spam atau email biasa (bukan spam). Terdapat 57 variabel yang akan menjadi prediktor yang terdiri atas:

• 48 variabel yang menggambarkan frekuensi relatif banyaknya kata ‘make’, ‘addres’, ‘all’, dan lain-lain

• 6 variabel yang menggambarkan banyaknya karakter ;, (, [, !, $, dan #

• 3 variabel yang menggambarkan keberadaan huruf kapital dalam badan email

Beberapa package yang digunakan dalam materi ini kita load terlebih dahulu sehingga fungsi-fungsi di dalamnya siap untuk kita gunakan.

library(ggstatsplot)
library(ggpubr)
library(dplyr)
library(rsample)
library(rpart)
library(rpart.plot)
library(yardstick)
library(ranger)
library(xgboost)
library(Ckmeans.1d.dp)

Pengenalan Lebih Jauh Kondisi Data

Sudah disampaikan sebelumnya bahwa dataset memuat 4601 amatan. Sebaran banyaknya email spam dan bukan spam dapat diperoleh pada pie-chart di bawah ini. Sekitar 39% amatan adalah email spam, dan 61% sisanya bukan spam.

ggpiestats(data = spam_data,
           x = "spam",
           label = "both",
           results.subtitle = FALSE )

Nilai rataan variabel prediktor dapat dilihat di bawah ini.

rataanX <- colMeans(select(spam_data,-spam))
data.rataan <- data.frame(
  x = names(rataanX),
  rataan = rataanX
)
ggbarplot(data = data.rataan[1:54,],
          x = "x",
          y="rataan",
          fill="steelblue",
          color = "steelblue",
          orientation="horizontal",
          ggtheme = theme_pubr(base_size = 5))

ggbarplot(data = data.rataan[55:57,],
          x = "x",
          y="rataan",
          fill="steelblue",
          color = "steelblue",
          orientation="horizontal",
          ggtheme = theme_pubr(base_size = 9))

Selanjutnya, kita barangkali tertarik untuk melihat bagaimana keterkaitan antara variabel prediktor dengan variabel respon. Di bawah ini kita coba untuk melihat apakah ada keterkaitan keberadaan kata “MONEY” pada email dengan kejadian email spam menggunakan stacked-bar chart. Sebelum membuat grafik, kita siakan dulu sebuah variabel baru yang nilainya “tidak ada” apabila kata money tidak ada pada email (dicirikan dengan kondisi word_freq_money == 0) serta bernilai ‘ada’ jika sebaliknya.

spam_data0 <- spam_data
spam_data0$word_money <- ifelse(spam_data0$word_freq_money==0,"Tidak Ada","Ada")
ggbarstats(data = spam_data0,
            x="spam",
            y="word_money",
            results.subtitle = FALSE)

Diagram di atas menginformasikan bahwa pada kelompok email yang memuat kata ‘money’, sebesar 93%-nya merupakan email spam. Sedangkan email-email yang tidak memuat kata tersebut, yang merupakan spam hanya sebesar 29%. Ketimpangan persentase/proporsi email spam pada kedua kelompok email berdasarkan keberadaan kata ‘money’ ini menegaskan bahwa keberadaan kata money penting untuk mendeteksi email spam.

Prosedur pemeriksaan yang sama kita coba jalankan untuk melihat efek dari keberadaan kata “RE”. Di bawah ini adalah perintah yang sama dengan sebelumnya untuk menghasilkan stacked-bar chart variabel keberadaan kata ‘re’ dengan status email spam.

spam_data0$word_re <- ifelse(spam_data0$word_freq_re==0,"Tidak Ada","Ada")
ggbarstats(data = spam_data0,
            x="spam",
            y="word_re",
            results.subtitle = FALSE)

Bar chart di atas menunjukkan bahwa proporsi email spam pada himpunan email yang memuat kata ‘re’ dan yang tidak memuat kata tersebut relatif tidak jauh berbeda. Karenanya kita dapat mengatakan bahwa variabel prediktor ini bukanlah variabel yang penting untuk digunakan membedakan email spam dari yang bukan.

Penyusunan Model Prediktif

Pada bagian ini akan diberikan ilustrasi penyusunan model prediktif dengan tiga pendekatan, yaitu: (1) Pohon Klasifikasi, (2) Random Forest, dan (3) Gradient Boosting. Untuk setiap pendekatan, juga diberikan ilustrasi evaluasi kinerja model yang dihasilkan dalam melakukan prediksi. Agar dapat daengan objektif membandingkan kinerjanya, pada ilustrasi ini dataset awal kita pisah secara acak menjadi dua bagian dengan perbandingan 80:20. Bagian yang 80% akan digunakan sebagai data latih (training dataset) untuk memperoleh model prediktif, dan bagian yang 20% akan digunakan sebagai data uji (testing set).

Berikut ini adalah perintah/program untuk melakukan pemisahan data secara acak berstrata untuk menjamin persentase email spam di data latih dan data uji relatif sama.

#levels(spam_data$spam)
spam_data$spam <- relevel(spam_data$spam,
                          ref = "yes" )
#levels(spam_data$spam)

set.seed(10)
define_split <- initial_split(data = spam_data,
                              prop = 0.8,
                              strata = "spam")

latih <- training(define_split)
uji <- testing(define_split)

dim(latih)

[1] 3680   58

prop.table(table(latih$spam))

      yes        no 
0.3940217 0.6059783 

dim(uji)

[1] 921  58

prop.table(table(uji$spam))

      yes        no 
0.3941368 0.6058632 

Berdasarkan output di atas kita peroleh bahwa data latih berisi 3680 amatan, sedangkan data uji berisi 921 amatan. Proporsi email spam di kedua dataset tersebut relatif sama yaitu 39.4%.

Pemodelan Pohon Klasifikasi

Berikut ini diberikan program pembentukan pohon klasifikasi menggunakan fungsi rpart. Ada dua pohon yang dibentuk. Masing-masing menggunakan hyperparameter minsplit yang berbeda. Pohon pertama (diberi nama pohon1) menggunakan minsplit=50, sedangkan pohon kedua menggunakan minsplit=100.

pohon1 <- rpart(spam~.,
                data = latih,
                control = rpart.control(cp=0, minsplit=50), 
                method = "class") 
rpart.plot(pohon1,
           type = 2,
            faclen = -1,
           extra=104,
            tweak = 1)

pohon2 <- rpart(spam~.,
                data = latih,
                control = rpart.control(cp=0, minsplit=100), 
                method = "class") 
rpart.plot(pohon2,
           type = 2,
            faclen = -1,
           extra=104,
            tweak = 1)

Minsplit adalah banyaknya amatan minimum dalam suatu simpul yang dipersyaratkan agar proses splitting dapat dilakukan. Mengatur minsplit yang besar akan membuat algoritma splitting lebih cepat berhenti sehingga pohon yang dihasilkan cenderung lebih kecil dan sederhana. Sebaliknya, minsplit yang kecil akan menyebabkan model pohon menjadi lebih kompleks. Berdasarkan gambar di atas, jelas bahwa pohon1 berukuran lebih besar karena dihasilkan dari algoritma yang menggunakan minsplit lebih kecil.

Dari kedua pohon tersebut, mana yang memberikan ketepatan prediksi yang lebih baik? Kita akan coba melihat perbandingan kinerjanya dengan perintah-perintah di bawah ini.

Proses pengukuran kinerja dilakukan dengan tahapan sebagai berikut. Pertama, model yang ada digunakan untuk melakukan prediksi bagi amatan-amatan yang ada pada data uji. Kemudian, hasil prediksi tersebut dibandingkan dengan data aktualnya. Beberapa statistik kinerja selanjutnya dihitung untuk melihat ketepatan prediksi tersebut seperti akurasi, sensitiviti, spesifisitas, dan sebagainya.

PRogram di bawah ini bekerja untuk memperoleh prediksi amatan data uji, dan kemudian membuat dataframe yang berisi kelas aktual, kelas prediksi, dan peluang suatu email masuk kategori spam berdasarkan model pohon1.

prob_pohon1 <- predict(pohon1, newdata=uji,
                      type = "prob")[,1]
pred_pohon1 <- ifelse(prob_pohon1>=0.5,"yes","no")
eval_pohon1 <- data.frame(aktual=uji$spam,
                          prediksi=pred_pohon1,
                          peluang=prob_pohon1)
eval_pohon1$prediksi <- factor(eval_pohon1$prediksi,
                               levels=c("yes","no"))

Kolom kelas aktual dan kelas prediksi selanjutnya dijadikan dasar untuk memperoleh confusion matrix atau tabel klasifikasi. Matriks konfusi itu yang kemudian dijadikan dasar perhitungan akurasi, sensitivitas, dan lain-lain.

conf_pohon1 <- conf_mat(eval_pohon1,
         truth = aktual,
         estimate = prediksi)
autoplot(conf_pohon1,type = "heatmap")+
  scale_fill_gradient(low = "#F4AFAB",high = "#EE847E")

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

summary(conf_pohon1)

# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.903
 2 kap                  binary         0.796
 3 sens                 binary         0.851
 4 spec                 binary         0.937
 5 ppv                  binary         0.898
 6 npv                  binary         0.906
 7 mcc                  binary         0.797
 8 j_index              binary         0.789
 9 bal_accuracy         binary         0.894
10 detection_prevalence binary         0.374
11 precision            binary         0.898
12 recall               binary         0.851
13 f_meas               binary         0.874

Berdasarkan output di atas, tampak bahwa nilai akurasi yang diperoleh berdasarkan prediksi pada data uji adalah 90.33%.

Selain menggunakan statistik berdasarkan matriks konfusi di atas, kita juga dapat menilai kinerja prediksi pohon klasifikasi menggunakan Kurva ROC. Luas di bawah kurva (AUC = area under the curve) dapat dijadikan salah satu ukuran tambahan mengenai kinerja model. Pada kasus kita, pohon1 memiliki nilai AUC sebesar 0.950.

roc_pohon1 <- roc_curve(data = eval_pohon1,
          truth = aktual,
          peluang)
autoplot(roc_pohon1)

roc_auc(data = eval_pohon1,
        truth = aktual,
        peluang)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.950

**Evaluasi kinerja pohon2*

Sekarang mari kita lakukan hal yang sama untuk mengukur kinerja prediksi dari pohon2. Meskipun memiliki bentuk model yang lebih sederhana, pohon ini memiliki akurasi sedikit lebih tinggi dibandingkan pohon1. Akurasi yang didapatkan adalah 91.3%. Nilai AUC-nya pun juga lebih besar yaitu 0.953.

prob_pohon2 <- predict(pohon2,
                      newdata=uji,
                      type = "prob")[,1]
pred_pohon2 <- ifelse(prob_pohon2>=0.5,"yes","no")
eval_pohon2 <- data.frame(aktual=uji$spam,
                          prediksi=pred_pohon2,
                          peluang=prob_pohon2)
eval_pohon2$prediksi <- factor(eval_pohon2$prediksi,
                                  levels=c("yes","no"))

conf_pohon2 <- conf_mat(eval_pohon2,
         truth = aktual,
         estimate = prediksi)
autoplot(conf_pohon2,type = "heatmap")+
  scale_fill_gradient(low = "#F4AFAB",high = "#EE847E")

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

summary(conf_pohon2)

# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.913
 2 kap                  binary         0.818
 3 sens                 binary         0.893
 4 spec                 binary         0.927
 5 ppv                  binary         0.888
 6 npv                  binary         0.930
 7 mcc                  binary         0.818
 8 j_index              binary         0.819
 9 bal_accuracy         binary         0.910
10 detection_prevalence binary         0.396
11 precision            binary         0.888
12 recall               binary         0.893
13 f_meas               binary         0.890

roc_pohon2 <- roc_curve(data = eval_pohon2,
          truth = aktual,
          peluang)

autoplot(roc_pohon2)

roc_auc(data = eval_pohon2,
        truth = aktual,
        peluang)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.953

Pemodelan Random Forest

Ilustrasi berikutnya adalah implementasi algoritma random forest menggunakan package ranger. Di bawah ini adalah contoh perintah program untuk menjalankan algoritma random forest yang berisi 500 buah pohon klasifikasi.

rf1 <- ranger(spam~.,data = latih, num.trees=500,
              importance = "impurity",
              probability = TRUE)

Kinerja prediksi model hasil algoritma random forest dapat dilihat pada hasil-hasil di bawah ini dimana akurasinya mencapai 95.4% dan AUC-nya sebesar 0.988.

prob_rf1 <- predict(rf1,
                      data=uji,
                      type = "response")
prob_rf1 <- prob_rf1$predictions[,1]
pred_rf1 <- ifelse(prob_rf1>=0.5,"yes","no")
eval_rf1 <- data.frame(aktual=uji$spam,
                          prediksi=pred_rf1,
                          peluang=prob_rf1)
eval_rf1$prediksi <- factor(eval_rf1$prediksi,
                               levels=c("yes","no"))

Confussion Matrix

conf_rf1 <- conf_mat(eval_rf1,
         truth = aktual,
         estimate = prediksi)
autoplot(conf_rf1,type = "heatmap")+
  scale_fill_gradient(low = "#F4AFAB",high = "#EE847E")

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

summary(conf_rf1)

# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.954
 2 kap                  binary         0.904
 3 sens                 binary         0.931
 4 spec                 binary         0.970
 5 ppv                  binary         0.952
 6 npv                  binary         0.956
 7 mcc                  binary         0.904
 8 j_index              binary         0.901
 9 bal_accuracy         binary         0.950
10 detection_prevalence binary         0.385
11 precision            binary         0.952
12 recall               binary         0.931
13 f_meas               binary         0.942

roc_rf1 <- roc_curve(data = eval_rf1,
          truth = aktual,
          peluang)

autoplot(roc_rf1)

roc_auc(data = eval_rf1,
        truth = aktual,
        peluang)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.988

Pemodelan Gradient Boosting

Algortima ketiga yang diimplementasikan pada materi ini adalah gradient boosting.

mat_train <- as.matrix(select(latih, -spam))
label_train <- pull(latih, spam)
label_train <- ifelse(label_train=="yes",1,0)
dtrain <- xgb.DMatrix(data = mat_train, label=label_train)

gb1 <- xgb.train(params = list(objective = "binary:logistic"),
                 data=dtrain,
                 nrounds = 2)

Menggunakan program di atas, model hasil algoritma gradient boosting disimpan dalam objek dengan nama gb1. Selanjutnya program di bawah ini dapat digunakan untuk memperoleh ukuran kinerja model gb1 sebagaimana yang juga digunakan pada model pohon klasifikasi dan random forest sebelumnya.

Kinerja prediksi model hasil algoritma gradient boosting dapat dilihat pada hasil-hasil di bawah ini dimana akurasinya adalah 92.4% dan AUC-nya sebesar 0.952.

mat_test <- as.matrix(select(uji, -spam))
dtest <- xgb.DMatrix(data = mat_test)

prob_gb1 <- predict(gb1,
                    newdata=dtest)
pred_gb1 <- ifelse(prob_gb1>=0.5,"yes","no")
eval_gb1 <- data.frame(aktual = uji$spam,
                          prediksi = pred_gb1,
                          peluang = prob_gb1)
eval_gb1$prediksi <- factor(eval_gb1$prediksi,
                               levels=c("yes","no"))

conf_gb1 <- conf_mat(eval_gb1,
         truth = aktual,
         estimate = prediksi)
autoplot(conf_gb1,type = "heatmap")+
  scale_fill_gradient(low = "#F4AFAB",high = "#EE847E")

Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.

summary(conf_gb1)

# A tibble: 13 × 3
   .metric              .estimator .estimate
   <chr>                <chr>          <dbl>
 1 accuracy             binary         0.924
 2 kap                  binary         0.839
 3 sens                 binary         0.879
 4 spec                 binary         0.953
 5 ppv                  binary         0.925
 6 npv                  binary         0.924
 7 mcc                  binary         0.840
 8 j_index              binary         0.832
 9 bal_accuracy         binary         0.916
10 detection_prevalence binary         0.375
11 precision            binary         0.925
12 recall               binary         0.879
13 f_meas               binary         0.901

roc_gb1 <- roc_curve(data = eval_gb1,
          truth = aktual,
          peluang)

autoplot(roc_gb1)

roc_auc(data = eval_gb1,
        truth = aktual,
        peluang)

# A tibble: 1 × 3
  .metric .estimator .estimate
  <chr>   <chr>          <dbl>
1 roc_auc binary         0.952

Ringkasan Perbandingan Kinerja

Ringkasan perbandingan kinerja model pohon klasifikasi, random forest, dan gradient boosting dapat disajikan pada tabel di bawah ini.

akurasi = as.numeric(c(summary(conf_pohon2)[1,3], 
            summary(conf_rf1)[1,3], 
            summary(conf_gb1)[1,3]))
AUC = as.numeric(c(
  roc_auc(data = eval_pohon2, truth = aktual,peluang)[1,3],
  roc_auc(data = eval_rf1, truth = aktual,peluang)[1,3],
  roc_auc(data = eval_gb1, truth = aktual,peluang)[1,3]))
kinerja = cbind(akurasi, AUC)
rownames(kinerja) = c("pohon", "random forest", "gradient boosting")
kinerja

                    akurasi       AUC
pohon             0.9131379 0.9534075
random forest     0.9543974 0.9877267
gradient boosting 0.9239957 0.9523189

Analisis Tingkat Kepentingan Variabel

imp_rf1 <-data.frame(
            variable=names(importance(rf1)),
            impurity=importance(rf1)
            )
rownames(imp_rf1) <- NULL

imp_rf1 <- slice_max(imp_rf1,
           n = 10,
           order_by=impurity)

ggbarplot(data = arrange(imp_rf1,impurity),
           x = "variable",
           y = "impurity",
          orientation="horizontal",
          fill="steelblue")

imp_gb1 <- xgb.importance(model = gb1)
xgb.ggplot.importance(importance_matrix = imp_gb1)