Dalam dokumen ini diberikan ilustrasi penggunaan R untuk menjalankan algoritma recursive partition untuk menghasilkan pohon regresi (regression tree). Ilustrasi akan diberikan mulai dari proses mengimport data menjadi dataframe, melakukan eksplorasi sederhana, membuat model pohon regresi berdasarkan data training, hingga mengevaluasi ketepatan prediksinya menggunakan data testing.
Data yang digunakan adalah data banyaknya kontainer yang dilayani Terminal JICT dan rata-rata Dwelling Time-nya setiap hari, yang tersedia pada URL: (kontak penulis).
Proses mengimport data dalam format file csv dapat dilakukan menggunakan fungsi read.csv() dan menyebutkan nama file beserta folder tempat menyimpannya.
data <- read.csv("D:/jict.csv")
Dengan perintah di atas, pada workspace-R kita memiliki data frame dengan mana data. Selanjutnya kita akan bekerja dengan dataframe ini.
Sebelum lebih jauh bekerja dengan data, kita perlu mengetahui informasi dasar dari data yang kita miliki. Informasi dasar tersebut meliputi ukuran/banyaknya data dan nama-nama kolom atau variabel yang ada di dalamnya.
FUngsi dim() dapat digunakan untuk menampilkan berapa banyak baris dan kolom pada dataframe kita, dan fungsi colnames() dapat digunakan untuk menampilkan nama-nama kolom pada data.
dim(data)
## [1] 872 21
colnames(data)
## [1] "terminal_id" "tanggal" "jumlah_kontainer"
## [4] "NON.MITA" "MITA" "HIJAU"
## [7] "MERAH" "KUNING" "NON"
## [10] "BORDER" "POST.BORDER" "BORDER_POST.BORDER"
## [13] "lartas.N" "lartas.Y" "lartas.NA"
## [16] "KARANTINA.TUMBUHAN" "KARANTINA.HEWAN" "KARANTINA.IKAN"
## [19] "PRENO" "NON_PRENO" "avg_dt"
Tampak bahwa ada sebanyak 872 baris pada data dengan kolom sebanyak 21 buah. Nama-nama kolom sudah diberikan oleh output di atas.
Sebelum melakukan pemodelan, ada beberapa variabel baru yang dibuat yaitu persentase kontainer yang dilayani yang tergolong pada kategori MITA, PRENO, LARTAS, KARANTINA, JALUR HIJAU, dan JALUR MERAH. Berikut cara memperoleh variabel baru tersebut.
data$p.mita <- data$MITA /data$jumlah_kontainer
data$p.preno <- data$PRENO/data$jumlah_kontainer
data$p.merah <- data$MERAH/data$jumlah_kontainer
data$p.hijau <- data$HIJAU/data$jumlah_kontainer
data$p.lartas <- data$lartas.Y/data$jumlah_kontainer
data$p.karantina <-(data$KARANTINA.HEWAN+data$KARANTINA.IKAN+
data$KARANTINA.TUMBUHAN)/data$jumlah_kontainer
Kemudian, eksplorasi data kita awali dengan melihat hubungan antara peubah prediktor dengan peubah respon. Yang akan menjadi peubah respon adalah rata-rata dwelling time ("avg_dt"), sedangkan variabel persentase berbagai jenis kontainer akan dijadikan variabel prediktor.
Di bawah ini cara memperoleh scatter plot antara prediktor dan respon. Sebagai ilustrasi hanya Persentase Kontainer PRENO dan JALUR HIJAU yang dibuat scatter plotnya.
plot(data$p.preno, data$avg_dt, ylim=c(0,8))
plot(data$p.hijau, data$avg_dt, ylim=c(0,8))
Cara lain untuk melihat hubungan di atas adalah dengan menggunakan boxplot, namun prediktornya perlu di-kategorikan terlebih dahulu. Berikut ini prosesnya.
k.preno <- cut(data$p.preno, breaks=c(0, 0.2, 0.3, 0.4, 0.5,
0.6, 0.7, 1))
boxplot(data$avg_dt ~ k.preno, ylim=c(0,8))
k.hijau <- cut(data$p.hijau, breaks=c(0.6, 0.8, 0.85,
0.9, 0.95, 1))
boxplot(data$avg_dt ~ k.hijau, ylim=c(0,8))
Nampak jelas bahwa semakin besar persentase kontainer PRENO dan JALUR HIJAU, rata-rata DT cenderung semakin kecil.
Proses pemodelan akan diawali dengan membagi data menjadi dua bagian secara acak, yaitu data latih (training set) dan data uji (testing set). Data latih adalah data yang akan digunakan untuk melakukan pemodelan, sedangkan data uji adalah data untuk melihat kinerja prediksinya.
Di bawwah ini adalah perintah untuk mengambil secara acak 600 baris untuk disimpan pada data training dan sisanya disimpan pada data testing.
set.seed(100)
acak <- sample(1:nrow(data), 600)
data.training <- data[acak,]
data.testing <- data[-acak,]
Ada dua dataframe baru yang dihasilkan dari perintah di atas yaitu data.training dan data.testing, yang banyaknya baris untuk masing-masing dapat diperoleh sebagai berikut.
dim(data.training)
## [1] 600 27
dim(data.testing)
## [1] 272 27
Selanjutnya kita dapat menjalankan algoritma pohon regresi menggunakan fungsi rpart() yang tersedia pada package dengan nama yang sama yaitu rpart.
Pada fungsi ini kita tinggal menyebutkan dataframe yang digunakan sebagai sumber data yaitu data.training. Ekspresi avg_dt ~ p.mita+p.preno+ p.merah + p.hijau + p.lartas+ p.karantina+jumlah_kontainer memiliki makna bahwa variabel target-nya adalah "avg_dt" dan variabel yang disebutkan digunakan sebagai prediktor. Antar variabel prediktor dipisahkan dengan tanda "+".
library(rpart)
pohon <- rpart(data=data.training, avg_dt ~ p.mita+p.preno+
p.merah + p.hijau + p.lartas+
p.karantina+jumlah_kontainer,
control=rpart.control(cp=0, minsplit=100))
Hasil dari pemodelan di atas tersimpan pada objek dengan nama pohon yang berisi model pohon regresi yang dihasilkan.
Perintah di bawah ini dapat digunakan untuk menampilkan pohon yang terbentuk. Terlihat bahwa partisi pertama adalah menggunakan variabel "p.preno".
pohon
## n= 600
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 600 2778.64100 3.149017
## 2) p.preno>=0.2034935 533 1281.19700 2.921501
## 4) p.preno>=0.4713983 248 154.58760 2.555363
## 8) p.hijau>=0.8814856 201 116.18250 2.422040
## 16) p.preno>=0.5527603 105 80.42338 2.260381
## 32) p.karantina>=0.06734603 63 12.83620 2.180000 *
## 33) p.karantina< 0.06734603 42 66.56956 2.380952 *
## 17) p.preno< 0.5527603 96 30.01377 2.598854 *
## 9) p.hijau< 0.8814856 47 19.55296 3.125532 *
## 5) p.preno< 0.4713983 285 1064.43400 3.240105
## 10) p.karantina< 0.1866407 252 166.17790 3.089087
## 20) p.hijau>=0.9087217 125 43.62584 2.746960
## 40) p.preno>=0.3471257 64 20.20406 2.589219 *
## 41) p.preno< 0.3471257 61 20.15853 2.912459 *
## 21) p.hijau< 0.9087217 127 93.51969 3.425827
## 42) p.mita>=0.1742588 88 52.79086 3.273750 *
## 43) p.mita< 0.1742588 39 34.10136 3.768974 *
## 11) p.karantina>=0.1866407 33 848.62070 4.393333 *
## 3) p.preno< 0.2034935 67 1250.37000 4.958955 *
Guna memudahkan melihat pohon klasifikasi yang terbentuk, kita dapat menggunakan fungsi rpart.plot() seperti di bawah ini.
library(rpart.plot)
rpart.plot(pohon, box.palette="RdBu", shadow.col="gray", nn=TRUE)
Bagaimana kualitas prediksinya? Berikut ini adalah proses untuk menghasilkan beberapa ukuran kinerja berdasarkan prediksinya di data testing, yaitu dengan menggunakan ukuran MAPE (mean absolute percentage error), yang diperoleh sebesar 28.5%.
prediksi <- predict(pohon, data.testing)
MAPE <- mean(abs(prediksi-data.testing$avg_dt)/data.testing$avg_dt)*100
MAPE
## [1] 28.53151
Di bawah ini proses menjalankan algoritma Random Forest di R dan menghitung kinerja prediksinya, dan diperoleh MAPE sebesar 25.8%.
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
forest <- randomForest(data=data.training, avg_dt ~ p.mita+p.preno+
p.merah + p.hijau + p.lartas+
p.karantina+jumlah_kontainer,
ntree=500)
prediksi.forest <- predict(forest, data.testing)
MAPE <- mean(abs(prediksi.forest-data.testing$avg_dt)/data.testing$avg_dt)*100
MAPE
## [1] 25.79675