Bahan Kuliah Pemodelan Klasifikasi - Prodi Statistika Terapan - IPB University

Pada dokumen ini akan dipaparkan bagaimana penggunaan R untuk melakukan pemodelan prediktif menggunakan pohon klasifikasi (classification tree).

Yang akan didiskusikan antara lain adalah:

Data yang akan digunakan sebagai ilustrasi adalah data lowbwt.csv yang diperoleh dari buku Applied Logistic Regression karya Hosmer dan Lemeshow.

Data ini memuat variabel-variabel:

Kita mulai dengan mengimport data menjadi dataframe di R, dengan cara berikut:

data <- read.csv('D:/lowbwt.csv', header = TRUE)

memperoleh model pohon klasifikasi

Andaikan saja kita ingin melakukan pemodelan pohon klasifikasi dengan variabel target/respon adalah “low” dan variabel prediktornya: * lwt * smoke * ht * ui

Fungsi yang dapat digunakan untuk memperoleh model pohon klasifikasi adalah rpart() yang ada pada package rpart. Pada fungsi itu kita harus menambahkan opsi method = ‘class’ sehingga fungsi akan mengenali variabel respon sebagai variabel kategorik dan algoritma rpart (recursive partition) akan menghasilkan pohon klasifikasi. Sebagai catatan, untuk variabel respon yang numerik, rpart juga dapat menjalankan algoritma untuk menghasilkan pohon regresi.

Berikut ini program untuk menghasilkan pohon klasifikasinya. Perhatikan cara penulisan modelnya, yang serupa dengan cara penulisan model regresi di fungsi lm() maupun glm() dan fungsi-fungsi lainnya.

Opsi minsplit=30 mengandung pengertian bahwa jika ada node yang berukuran kurang dari 30 amatan maka algoritma dihentikan. Kita pending dulu pembahasan mengenai cp dan akan diskusikan selanjutnya. Sementara ini kita gunakan cp=0.

library(rpart)
pohon <- rpart(low ~ lwt + smoke + ht + ui, data=data,
               method='class',
               control=rpart.control(minsplit = 30, cp=0))

Model yang diperoleh, disimpan dalam objek dengan nama pohon yang kalau dilihat isinya adalah sebagai berikut:

pohon
## n= 189 
## 
## node), split, n, loss, yval, (yprob)
##       * denotes terminal node
## 
##  1) root 189 59 >= 2500 g (0.3121693 0.6878307)  
##    2) lwt< 106 37 17 < 2500 g (0.5405405 0.4594595)  
##      4) lwt>=100.5 13  4 < 2500 g (0.6923077 0.3076923) *
##      5) lwt< 100.5 24 11 >= 2500 g (0.4583333 0.5416667) *
##    3) lwt>=106 152 39 >= 2500 g (0.2565789 0.7434211)  
##      6) ui=Yes 20 10 < 2500 g (0.5000000 0.5000000) *
##      7) ui=No 132 29 >= 2500 g (0.2196970 0.7803030)  
##       14) ht=Yes 10  4 < 2500 g (0.6000000 0.4000000) *
##       15) ht=No 122 23 >= 2500 g (0.1885246 0.8114754) *

Jika kita perhatikan dengan seksama hasil tersebut kita dapat membayangkan bentuknya kira-kira sebagai berikut. Rootnode berisi 189 amatan terbagi menjadi dua dengan aturan apakah lwt<106 atau lwt>106. Node dengan lwt<106 berisi 37 amatan, dan node lwt >=106 berisi 152 amatan.

Selanjutnya node lwt<36 terbagi menjadi dua yaitu apakah lwt>=100 dan lwt<100 yang masing-masing berisi 13 dan 24 amatan. KEduanya tidak lagi terpisah, yang salah satu alasannya karena kita stop dengan opsi minsplit=30. Keduanya berisi kurang dari 30 amatan sehingga tidak dilanjutkan lagi pemisahan oleh rpart dan menjadi terminal node.

Perhatikan angka lain yang diberikan… misal pada root node. Seperti yang diberikan pada baris awal output:

node), split, n, loss, yval, (yprob)

nilai 189 adalah n, nilai 59 adalah frekuensi amatan yang salah klasifikasi, >= 2500 atau yval adalah prediksi pada node tersebut yaitu kelas yang dominan, dan (0.3121693 0.6878307) masing-masing adalah proporsi dari yang “< 2500g” dan “>= 2500 g”. Yang lebih dominan adalah yang “>= 2500 g”, sehingga itu yang muncul pada yval.

Tentu saja akan lebih mudah mengetahui bentuk pohon klasifikasi yang terbentuk dengan cara visualisasi yang lebih baik. Berikut ini akan diberikan salah satu cara menggambar pohon klasifikasi yang telah diperoleh

Menggambar pohon klasifikasi

Fungsi yang dapat digunakan untuk visualisasi pohon klasifikasi adalah rpart.plot() yang ada pada package rpart.plot. Kita tinggal menyebutkan nama modelnya dan opsi nomor pilihan apa saja yang ingin ditampilkan di setiap node.

Opsi extra=4 meminta fungsi menampilkan proporsi dari masing-masing kelas variabel target dan kelas apa yang dominan. Berikut ini ilustrasinya.

library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 3.5.3
rpart.plot(pohon, extra=4)

Melakukan Prediksi menggunakan Pohon Klasifikasi

Pohon yang diperoleh selanjutnya dapat digunakan untuk melakukan prediksi kejadian apakah bayi lahir dengan berat badan rendah (< 2500 g) ataukah normal (>= 2500 g).

Misalnya saja ada seorang ibu dengan informasi berikut:

  • lwt = 140
  • smoke = “Yes”
  • ht = “No”
  • ui = “No”

maka mengikuti pohon yang ada pertama akan dicek apakah lwt<106? Jika ya, ambil jalur kiri. Jika tidak, ambil jalur kanan. Karena lwt = 140 maka berarti ambil jalur kanan.

Selanjutnya ketemu dengan pertanyaan apakah ui=“Yes”? Jika ya, ambil jalur kiri. Jika tidak, ambil jalur kanan. Karena untuk data kita ui=“No”, tidak “Yes”, maka ambil jalur kanan.

Kemudian bertemu dengan pertanyaan apakah ht=“Yes”? Jika ya, ambil jalur kiri. Jika tidak, ambil jalur kanan. Data kita meniliki nilai ht=“No”, maka ambil jalur kanan.

Sekarang kita ketemu dengan terminal node, dan diprediksi memiliki Peluang(< 2500 g) = 0.19 dan Peluang (>= 2500 g) = 0.81. Dengan menggunakan threshold peluang 0.5, amatan ini akan diprediksi memiliki bayi dengan berat tergolong “>= 2500 g” seperti yang ada pada label terminal node-nya.

Terminal node yang berwarna hijau menunjukkan kalau prediksinya adalah “>= 2500 g”. Semakin gelap hijaunya, semakin tinggi peluangnya. Sedangkan untuk kelas sebaliknya, “< 2500 g”, terminal node diberikan warna biru.

Menggunakan R, kita dapat melakukan prediksi menggunakan fungsi predict dengan input nama model dan data yang diprediksi. Pastikan data yang diprediksi berupa dataframe dengan nama-nama variabel sama dengan nama yang ada pada model.

Berikut contoh melakukan prediksi pada individu dengan data lwt = 140, smoke = “Yes”, ht = “No”, ui=“No”.

bayi.baru <- data.frame(lwt = 140, smoke = "Yes", ht = "No", ui="No")
predict(pohon, bayi.baru)
##    < 2500 g >= 2500 g
## 1 0.1885246 0.8114754

Beberapa hyperparameter dari pohon klasifikasi

Kita telah melihat bagaimana pohon yang dihasilkan pada saat kita gunakan opsi minsplit=30. Minsplit adalah salah satu hyperparameter pada algoritma pohon klasifikasi dan perlu kita tentukan di awal.

Seperti yang telah dijelaskan, minsplit adalah opsi untuk menentukan berapa ukuran node minimal yang diperbolehkan untuk melakukan splitting/pemisahan. Dengan opsi minsplit=30 maka node-node yang berisi kurang dari 30 amatan tidak akan dilanjutkan proses splitting-nya dan algoritma berhenti.

Seandainya nilai minsplit kita perkecil, maka kita akan memperoleh pohon yang lebih kompleks, karena node yang sedikit isinya masih boleh untuk di-split. Berikut tampilan pohon dengan minsplit=10.

pohon10 <- rpart(low ~ lwt + smoke + ht + ui, data=data,
               method='class',
               control=rpart.control(minsplit = 10, cp=0))
rpart.plot(pohon10, extra=4)

Selain minsplit, ada juga beberapa opsi/hyperparameter lain, diantaranya:

opsi keterangan
minsplit the minimum number of observations that must exist in a node in order for a split to be attempted. Default = 20.
minbucket the minimum number of observations in any terminal node. If only one of minbucket or minsplit is specified, the code either sets minsplit to minbucket*3 or minbucket to minsplit/3, as appropriate.
maxdepth Set the maximum depth of any node of the final tree, with the root node counted as depth 0. Default = 30.

Optimasi Hyperparameter

Berapa nilai hyperparameter minsplit yang sebaiknya kita gunakan?

Untuk menjawabnya kita bisa menerapkan proses validasi atau validasi silang dengan cara melakukan percobaan beberapa kali pembuatan pohon klasifikasi dengan minsplit yang berbeda dan kita tentukan mana yang optimal.

library(caret)
## Loading required package: lattice
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 3.5.3
set.seed(100)
akurasi.semua <- NULL

for(ulangan in 1:30){
acak <- createDataPartition(data$low, p=0.7, list=FALSE)
data.training <- data[acak,]
data.testing <- data[-acak,]

for (k in 1:40){
pohon <- rpart(low ~ lwt + smoke + ht + ui 
               + age +ptl + ftv, 
               data=data.training,
               method='class',
               control=rpart.control(minsplit = k, cp=0))
prediksi.prob <- predict(pohon, data.testing)
prediksi <- ifelse(prediksi.prob > 0.5, ">= 2500 g", "< 2500 g")[,2]
akurasi <- mean(prediksi == data.testing$low)
akurasi.semua <- rbind(akurasi.semua, c(k, akurasi))
}
}
mean.akurasi <- tapply(akurasi.semua[,2], akurasi.semua[,1], mean)
plot(names(mean.akurasi),mean.akurasi, type="b", xlab="minsplit", ylab="rata-rata akurasi data testing")

Menggunakan rata-rata dari 30 kali ulangan, tampaknya minsplit yang optimal adalah 25.

Karenanya, pada untuk memperoleh model akhir selanjutnya pada data ini gunakan seluruh data dengan minsplit=25.