1. Pendahuluan

Latar Belakang

Berlian adalah salah satu ciptaan alam yang paling berharga, indah, dan langkah. Berlian merupakan zat alami terkeras yang tersusun dari karbon murni. Berlian juga merupakan batu permata paling populer. Berlian memiliki sejumlah aplikasi penting dalam dunia industri.Berlian 58 kali lebih keras daripada mineral paling keras kedua di bumi.

Diamond berasal dari bahasa Yunani adamas yang berarti ‘tidak terkalahkan’. Ini sesuai dengan sifat aslinya yaitu batu alam paling keras di bumi. Dalam bahasa Latin disebut juga adamare yang artinya mencintai dengan sepenuh hati.

Berlian ditemukan di tiga jenis endapan: kerikil aluvial, ladang glasial, dan pipa kimberlite. Pipa kimberlite (seperti yang ada di Kimberley, Afrika Selatan) terbentuk dari intrusi magma ke kerak bumi yang membawa berlian serta batuan dan mineral lainnya dari mantel bumi. Pipa-pipa itu sendiri seringkali berumur kurang dari 100 juta tahun. Namun, berlian yang mereka bawa terbentuk 1 hingga 3,3 miliar tahun dengan kedalaman lebih dari sekitar 75 mil (120 km). Sedangkan berlian jenis kerikil fluvial dan glasial berasal dari lepasan erosi fluvial dan glasial yang terjadi di matriks kimberlite. Kemudian disimpan kembali di sungai atau di glasial till (merupakan endapan yang bergerak menjauh dari inti atau gletser).

Memprediksi harga suatu berlian dengan beberapa variabel indikator menjadi sangat penting agar setiap orang yang membeli berlian tidak tertipu dengan kualitas berlian yang dibelinya. Dengan demikian pada project ini kami akan membangun model untuk memprediksi harga berlian berdasarkan ciri/karakteristik tertentu.

Goal

Tujuan dari project ini adalah membangun model yang dapat memprediksi harga suatu diamond berdasarkan karakteristik diamond tersebut.

Algoritma

Adapun algoritma yang akan digunakan dalam pemodelan ini adalah sebagai berikut :

2. Package

Berikut ini adalah package yang dipakai dalam project ini.

library(gbm)

## Loaded gbm 2.1.8

library(mlr3tuning)

## Loading required package: mlr3

## Loading required package: paradox

library(gam)

## Loading required package: splines

## Loading required package: foreach

## Loaded gam 1.20

library(rpart)
library(pROC)

## Type 'citation("pROC")' for a citation.

## 
## Attaching package: 'pROC'

## The following objects are masked from 'package:stats':
## 
##     cov, smooth, var

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.3     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x purrr::accumulate() masks foreach::accumulate()
## x dplyr::filter()     masks stats::filter()
## x dplyr::lag()        masks stats::lag()
## x purrr::when()       masks foreach::when()

library(mlr3verse)
library(mlr3extralearners)

## 
## Attaching package: 'mlr3extralearners'

## The following objects are masked from 'package:mlr3':
## 
##     lrn, lrns

library(precrec)

## 
## Attaching package: 'precrec'

## The following object is masked from 'package:pROC':
## 
##     auc

library(adabag)

## Loading required package: caret

## Loading required package: lattice

## 
## Attaching package: 'caret'

## The following object is masked from 'package:purrr':
## 
##     lift

## Loading required package: doParallel

## Loading required package: iterators

## Loading required package: parallel

library(ROCR)
library(ROCit)
library(magrittr)

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

library(visdat)
library(naniar)
library(UpSetR)

## 
## Attaching package: 'UpSetR'

## The following object is masked from 'package:lattice':
## 
##     histogram

library(laeken)
library(vcd)

## Loading required package: grid

library(VIM)

## Loading required package: colorspace

## 
## Attaching package: 'colorspace'

## The following object is masked from 'package:pROC':
## 
##     coords

## VIM is ready to use.

## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues

## 
## Attaching package: 'VIM'

## The following object is masked from 'package:datasets':
## 
##     sleep

library(sm)

## Package 'sm', version 2.2-5.7: type help(sm) for summary information

library(ggplot2)
library(dplyr)
library(mlbench)
library(caret)
library(mlr3verse)
library(mlr3fselect)
library(DataExplorer)
library(skimr)

## 
## Attaching package: 'skimr'

## The following object is masked from 'package:naniar':
## 
##     n_complete

## The following object is masked from 'package:mlr3':
## 
##     partition

library(corrplot)

## corrplot 0.92 loaded

library(leaps)
library(olsrr)

## 
## Attaching package: 'olsrr'

## The following object is masked from 'package:datasets':
## 
##     rivers

library(kableExtra) #Tampilan Tabel

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

library(agricolae)  #Pemeriksaan Asumsi
library(lmtest)     #Untuk pengecekan asumsi

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(car)        #Untuk pengecekan asumsi

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:ROCit':
## 
##     logit

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

library(tseries)    #Untuk pengecekan asumsi

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(glmnet)

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## Loaded glmnet 4.1-2

library(glmnetUtils)

## 
## Attaching package: 'glmnetUtils'

## The following objects are masked from 'package:glmnet':
## 
##     cv.glmnet, glmnet

library(broom)
library(ggpubr)
library(modelr)

## 
## Attaching package: 'modelr'

## The following object is masked from 'package:broom':
## 
##     bootstrap

## The following object is masked from 'package:mlr3':
## 
##     resample

library(precrec)
library(adabag)
library(rpart.plot)
library(mice)

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(caret)
library(randomForest)

## randomForest 4.6-14

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:ggplot2':
## 
##     margin

library(ROSE)

## Loaded ROSE 0.0-4

3. Pre-processing Data

Sebelum data yang ada kita modelkan, terlebih dahulu kita siapkan data yang akan digunakan (training & Testing) dan cleaning data (jika diperlukan).

data3 <- read.csv2("D:/Magister IPB/Kuliah/Semester 2/STA582_Pembelajaran_Mesin_Statistika/Kuliah/Pertemuan 5/diamonds.csv",stringsAsFactors = TRUE)
data3<-data3[,-1]
dim(data3)

## [1] 53940    10

str(data3)

## 'data.frame':    53940 obs. of  10 variables:
##  $ carat  : num  0.23 0.21 0.23 0.29 0.31 0.24 0.24 0.26 0.22 0.23 ...
##  $ cut    : Factor w/ 5 levels "Fair","Good",..: 3 4 2 4 2 5 5 5 1 5 ...
##  $ color  : Factor w/ 7 levels "D","E","F","G",..: 2 2 2 6 7 7 6 5 2 5 ...
##  $ clarity: Factor w/ 8 levels "I1","IF","SI1",..: 4 3 5 6 4 8 7 3 6 5 ...
##  $ depth  : num  61.5 59.8 56.9 62.4 63.3 62.8 62.3 61.9 65.1 59.4 ...
##  $ table  : num  55 61 65 58 58 57 57 55 61 61 ...
##  $ price  : int  326 326 327 334 335 336 336 337 337 338 ...
##  $ x      : num  3.95 3.89 4.05 4.2 4.34 3.94 3.95 4.07 3.87 4 ...
##  $ y      : num  3.98 3.84 4.07 4.23 4.35 3.96 3.98 4.11 3.78 4.05 ...
##  $ z      : num  2.43 2.31 2.31 2.63 2.75 2.48 2.47 2.53 2.49 2.39 ...

glimpse(data3)

## Rows: 53,940
## Columns: 10
## $ carat   <dbl> 0.23, 0.21, 0.23, 0.29, 0.31, 0.24, 0.24, 0.26, 0.22, 0.23, 0.~
## $ cut     <fct> Ideal, Premium, Good, Premium, Good, Very Good, Very Good, Ver~
## $ color   <fct> E, E, E, I, J, J, I, H, E, H, J, J, F, J, E, E, I, J, J, J, I,~
## $ clarity <fct> SI2, SI1, VS1, VS2, SI2, VVS2, VVS1, SI1, VS2, VS1, SI1, VS1, ~
## $ depth   <dbl> 61.5, 59.8, 56.9, 62.4, 63.3, 62.8, 62.3, 61.9, 65.1, 59.4, 64~
## $ table   <dbl> 55, 61, 65, 58, 58, 57, 57, 55, 61, 61, 55, 56, 61, 54, 62, 58~
## $ price   <int> 326, 326, 327, 334, 335, 336, 336, 337, 337, 338, 339, 340, 34~
## $ x       <dbl> 3.95, 3.89, 4.05, 4.20, 4.34, 3.94, 3.95, 4.07, 3.87, 4.00, 4.~
## $ y       <dbl> 3.98, 3.84, 4.07, 4.23, 4.35, 3.96, 3.98, 4.11, 3.78, 4.05, 4.~
## $ z       <dbl> 2.43, 2.31, 2.31, 2.63, 2.75, 2.48, 2.47, 2.53, 2.49, 2.39, 2.~

# membagi data training dan data testing
set.seed(134)
Sample <- sample (1:53940, 43152) 
testing <- data3 [Sample, ]
training <- data3 [-Sample, ]

dim(testing)

## [1] 43152    10

dim(training)

## [1] 10788    10

sum_mis3<-miss_var_summary(training)
sum_mis_plot3<-head(sum_mis3,11)
sum_mis_plot3

#Sebaran data training (price) sebelum transformasi
ggpubr::gghistogram(data = training,x = "price",fill = "darkgreen")+scale_y_continuous(expand = c(0,0))

## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.

plot_intro(training)

#Eksplorasi data
DataExplorer::plot_scatterplot(data = training,
                               by = "price",nrow = 3,ncol = 3,geom_point_args = list(color="Steelblue"))

#Boxplot price sebelum transformasi
boxplot(training$price)

#Transformasi variabel price
training<-training %>% 
  mutate(price=log(price))

#Sebaran data training (price) setelah transformasi
ggpubr::gghistogram(data = training,x = "price",fill = "darkgreen")+scale_y_continuous(expand = c(0,0))

## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.

#Boxplot price setelah transformasi
boxplot(training$price)

testing<-testing %>% 
  mutate(price=log(price))

ggpubr::gghistogram(data = testing,x = "price",fill = "darkgreen")+scale_y_continuous(expand = c(0,0))

## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.

4. Pemodelan

Pada pertemuan ini metode yang akan digunakan dalam pomodelan adalah regresi gradient boost.

Gradient Boosting

Pertama kita akan bangun dengan menggunakan model gradient boosting, untuk kemudian akan kita uji akurasinya dengan menggunakan fungsi predict yang mengacu pada nilai RMSE dan MAE.

Membuat Model Awal

Pada tahapan ini kita akan membangun model gradient boosting dengan menggunakan hyperparameter awal.

gbm <- gbm(formula = price ~ .,
             distribution = "gaussian",
             data =training,
             n.trees = 500,  
            interaction.depth = 9,
            cv.folds = 5,
            n.cores = NULL, # will use all cores by default
            verbose = FALSE)

print(gbm)

## gbm(formula = price ~ ., distribution = "gaussian", data = training, 
##     n.trees = 500, interaction.depth = 9, cv.folds = 5, verbose = FALSE, 
##     n.cores = NULL)
## A gradient boosted model with gaussian loss function.
## 500 iterations were performed.
## The best cross-validation iteration was 498.
## There were 9 predictors of which 9 had non-zero influence.

Variabel Importance

Berikut ini akan ditampilkan urutan variabel yang paling berpengaruh terhadap variabel respon secara terurut tersusun menurun.

sqrt(min(gbm$cv.error))

## [1] 0.09739018

summary(gbm)

Secara berturut-turut lima variabel yang paling berpengaruh terhadap variabel respon adalah y, carat, x, z, dan clarity.

Tuning Hyperparameter

Untuk mendapatkan model terbaik, perlu dilakukan tuning hyperparameter. Beberapa variasi dalam hyperparameter ditentukan secara manual, dan dalam hal ini setting tunning hyperparameter mengacu pada kedalaman pohon, jumlah pohon, dan jumlah sampel minimum pada setiap node.

hyper_grid <- expand.grid(
  interaction.depth = c(2:9), #Maximum depth of each tree
  optimal_trees = 0, 
  min_RMSE = 0,
  n.minobsinnode = c(2000, 1000, 500, 100,50,25)# minimum number of observations in the terminal nodes of the trees
)

nrow(hyper_grid) # Número total de combinaciones

## [1] 48

Sekarang kita mulai proses tunning hyperparameter untuk mendapatkan hyperparameter dengan model terbaik yang mengacu pada nilai RMSE nya.

for(i in 1:nrow(hyper_grid)) {
  set.seed(123)
  gbm.tune <- gbm(
    formula = price ~ .,
    distribution = "gaussian",
    data = training,
    n.trees = 1000,
    interaction.depth = hyper_grid$interaction.depth[i],
    n.minobsinnode = hyper_grid$n.minobsinnode[i],
    train.fraction = 0.75,
    n.cores = NULL,
    verbose = FALSE
  )
  
# Agregamos a hyper_grid la información que nos interesa que nos interesa
  hyper_grid$optimal_trees[i] <- which.min(gbm.tune$valid.error)
  hyper_grid$min_RMSE[i] <- sqrt(min(gbm.tune$valid.error))

}

hyper_grid %>% 
  dplyr::arrange(min_RMSE) %>%
  head(10)

min(hyper_grid$min_RMSE)

## [1] 0.1271677

Berdasarkan output di atas, dapat disimpulkan bahwa model ini akan menghasilkan nilai RMSE minimum saat kita menggunakan interaction.depth = 3, n-trees = 975, dan n.minobsinnode = 50. Hal ini terlihat dari nilai RMSE nya adalah yang paling kecil.

Final Model

Sekarang kita bangun ulang model dengan menggunakan model yang terbaik di atas.

# train GBM model
gbm.fit.final <- gbm(
  formula = price ~ .,
  distribution = "gaussian",
  data = training,
  n.trees = 975,
  interaction.depth = 3,
  n.minobsinnode = 50,
  n.cores = NULL,
  verbose = FALSE
  ) 

summary(gbm.fit.final, cBars = 10,
        method = relative.influence, las = 2)

Didapatkan secara terurut variabel importance dari model tsb.

5. Predict Data

Selanjutnya kita gunakan model terbaik tersebut untuk memprediksi data testing.

gbm_pred <- predict(object=gbm.fit.final, newdata=testing)

## Using 975 trees...

gbm_pred<-as.data.frame(gbm_pred)
a<-gbm_pred$gbm_pred
head(gbm_pred)

#Mengembalikan fungsi log natural
price_predict<-  2.718282^(a)
price_real<- 2.718282^(testing$price)

hasil_akhir<- as.data.table(cbind(price_predict,price_real))
hasil_akhir

MAE_Testing<-mean(abs(price_predict-price_real))
MAE_Testing

## [1] 294.3114

6. Kesimpulan

Adapun yang dapat disimpulkan dari simulasi di atas adalah sbb:

1. Model terbaik adalah regresi gradient boosting dengan hyperparameter tertentu.

2. Metode ini dapat sangat baik dalam menjelaskan hubungan antara variabel respon dan predictor.

3. Tuning hyperparameter adalah hal yang sangat penting untuk dilakukan.

4. Semakin presisi model pada data training, maka berpeluang akan overfit pada data test.

5. Model regresi gradient boosting ini lebih baik dibanding model random forest.

Referensi

Breheny P. (n.d.). Getting started with grpreg. GitHub Pages. https://pbreheny.github.io/grpreg/articles/getting-started.html

Hastie T. (2013, May 9). glmnet: Lasso and elastic-net regularization in R. Revolutions. https://blog.revolutionanalytics.com/2013/05/hastie-glmnet.html

Post J. (2014, September 29). LASSO, Ridge, and Elastic Net. https://www4.stat.ncsu.edu/~post/josh/LASSO_Ridge_Elastic_Net_-_Examples.html

Tibshirani, R., Saunders, M., Rosset, S., Zhu, J., & Knight, K. (2005). Sparsity and smoothness via the fused lasso. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 67(1), 91-108. https://doi.org/10.1111/j.1467-9868.2005.00490.x

Tripathy A. (2013, July 14). Regularization – Predictive Modeling Beyond Ordinary Least Squares Fit. ShatterLine Blog. https://shatterline.com/blog/2013/07/

Xiaotong C., Chen G., & Chong W. (n.d.). Statistical Learning and Data Mining Codes. Biostatistics - Academic Divisions - School of Public Health - University of Minnesota. https://www.biostat.umn.edu/~weip/course/dm/examples/