Pendahuluan

Tujuan

  • Melakukan cross validation untuk menghasilkan pemodelan antara tingkat efisiensi penggunaan bahan bakar (direpresentasikan oleh mpg atau miles per gallon) versus kekuatan mesin (horsepower) pada data ISLR::Auto dengan metode:

    1. Polynomial Regression
    2. Piecewise Constant
    3. Natural Cubic Splines

  • Pemodelan dilakukan untuk masing-masing subset data berdasarkan asal negara produsennya (Amerika, Eropa, Jepang).

Output

Mendapatkan model non-linear yang optimal berdasarkan masing-masing data.

Metodologi

Metode Cross Validaion

10-fold cross validation dengan meminimumkan RMSE (Root Mean Squared Error).

Pilihan k biasanya 5 atau 10, tetapi tidak ada aturan formal. Ketika k semakin besar, perbedaan ukuran antara himpunan pelatihan dan himpunan sampel ulang semakin kecil. Ketika perbedaan ini berkurang, bias dari teknik ini menjadi lebih kecil. (Page 70, Applied Predictive Modeling, 2013)

Skenario Tunning Parameter

Polynomial Regression

Derajat/degree polinomial awal: 2:4

d (derajat polinomial) yang digunakan biasanya hingga 4 karena d yang terlalu besar akan membuat kurva berbentuk tidak teratur dan tampak “liar” terutama pada bagian ujungnya.

Piecewise Constant

Jumlah breaks awal: 2:10data Auto, 2:6 untuk Data AutoEropa, dan 2:5 untuk Data AutoJepang

Natural Cubic Spline

Jumlah dan lokasi knots awal ditentukan berdasarkan basis yang terpilih 2:10

Tools (Package)

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'tibble' was built under R version 4.2.3
## Warning: package 'readr' was built under R version 4.2.3
## Warning: package 'dplyr' was built under R version 4.2.3
## Warning: package 'stringr' was built under R version 4.2.3
## Warning: package 'forcats' was built under R version 4.2.3
## Warning: package 'lubridate' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.3     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(splines)
library(rsample)
## Warning: package 'rsample' was built under R version 4.2.3
library(ISLR)
## Warning: package 'ISLR' was built under R version 4.2.3
library(cowplot)
## Warning: package 'cowplot' was built under R version 4.2.3
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
library(caret) # cross-validation
## Warning: package 'caret' was built under R version 4.2.3
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(gridExtra) # combining graphs
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(ggplot2)
library(dplyr)
library(purrr)
library(dplyr)
library(knitr)
## Warning: package 'knitr' was built under R version 4.2.3
library(DT)
## Warning: package 'DT' was built under R version 4.2.3
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.2.3
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(gam) # generalized additive models
## Warning: package 'gam' was built under R version 4.2.3
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 4.2.3
## 
## Attaching package: 'foreach'
## 
## The following objects are masked from 'package:purrr':
## 
##     accumulate, when
## 
## Loaded gam 1.22-2
library(rstatix) #Deteksi Outlier
## Warning: package 'rstatix' was built under R version 4.2.3
## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter

Hasil dan Pembahasan

Eksplorasi Data

Deksripsi Data

Auto MPG merupaka dataset mengenai penggunaan bahan bakar yang dihitung dengan satuan miles per gallon (mpg). mpg merupakan satuan nilai penggunaan bahan bakar sekaligus sebagai peubah respon pada dataset tersebut. Terdapat tujuh peubah prediktor yang memengaruhi penggunaan bahan bakar yakni sebagai berikut.

No Variabel Keterangan
1 Cylinder Silinder yang digunakan pada kendaraan
2 Displacement Besarnya perpindahan
3 Horsepower Tenaga kendaraan (Kekuatan Mesin)
4 Weight Berat kendaraan
5 Acceleration Percepatan kendaraan
6 Year Model Tahun Keluaran Kendaraan
7 Origin Asal Kendaraan 
head(Auto)
##   mpg cylinders displacement horsepower weight acceleration year origin
## 1  18         8          307        130   3504         12.0   70      1
## 2  15         8          350        165   3693         11.5   70      1
## 3  18         8          318        150   3436         11.0   70      1
## 4  16         8          304        150   3433         12.0   70      1
## 5  17         8          302        140   3449         10.5   70      1
## 6  15         8          429        198   4341         10.0   70      1
##                        name
## 1 chevrolet chevelle malibu
## 2         buick skylark 320
## 3        plymouth satellite
## 4             amc rebel sst
## 5               ford torino
## 6          ford galaxie 500

Data yang digunakan adalah Auto dari package ISLR yang terdiri dari 392 observasi dan 9 variabel. Namun, variabel yang digunakan kali ini hanya 3, yaitu:

  1. mpg : miles per gallon
  2. horsepower: Engine horsepower
  3. origin: Origin of car (1. American, 2. European, 3. Japanese)

Data Auto

Auto = Auto %>%  select(mpg, horsepower, origin)
datatable(Auto)
summary(Auto)
##       mpg          horsepower        origin     
##  Min.   : 9.00   Min.   : 46.0   Min.   :1.000  
##  1st Qu.:17.00   1st Qu.: 75.0   1st Qu.:1.000  
##  Median :22.75   Median : 93.5   Median :1.000  
##  Mean   :23.45   Mean   :104.5   Mean   :1.577  
##  3rd Qu.:29.00   3rd Qu.:126.0   3rd Qu.:2.000  
##  Max.   :46.60   Max.   :230.0   Max.   :3.000

Visualisasi Hubungan mpg terhadap horsepower

plotmpg <- function(x){
  ggplot(Auto) +
  geom_point(aes_string(x = tolower(x), y = "mpg"), color="dodgerblue2", alpha = 0.8) +
  ggtitle(paste0("Miles per Gallon vs ", x)) +
  xlab(x) + ylab("Miles per Gallon") +
  theme_bw() 
}
p <- plotmpg("horsepower")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
p

Ini merupakan scatter plot dari 392 kendaraan dimana dari masing-masing kendaraan diukur dua peubah/variabel, yaitu kekuatan mesin (horse power pada sumbu X) dan berapa jarak yang ditempuh kendaraan tersebut untuk setiap galon bahan bakarnya (mile per gallon pada sumbu y).

Berdasarkan scatter plot ini, secara umum terlihat bahwa semakin besar kekuatan mesin (horsepower) suatu kendaraan, semakin pendek jarak yang ditempuh kendaraan itu per galon bahan bakarnya. Artinya semakin kuat mesinnya, maka semakin boros. Dapat dilihat juga bahwa hubungan antara kekuatan mesin kendaraan (pada sumbu X) dan jarak yang dapat ditempuh oleh kendaraan per galon bahan bakarnya (pada sumbu Y) cenderung bersifat tidak linear.

Beberapa hal yang kurang menyenangkan dari model linier pada kasus ini, yaitu:

  • Model linier yang kita gunakan tidak mengikuti pola hubungan pada data. Sebelah kiri curam, sebelah kanan cenderung landai sehingga prediksinya tidak cukup baik.

  • Dapat diperoleh nilai penduga peubah respon yang “tidak masuk akal”. Jika garis linier diteruskan, untuk kekuatan mesin tertentu, akan diperoleh respon yang negatif sehingga tidak mungkin terjadi.

  • Sifat dan perilaku sisaan cenderung tidak sesuai harapan. Jika wilayah dibagi menjadi tiga bagian, bagian yang tengah, eror cenderung negatif karena garis cenderung ada di atas titik-titik yang sebenarnya. Sedangkan wilayah yang kanan sisaan cenderung positif karena garis cenderung ada di bawah titik-titik yang sebenarnya. Padahal biasanya lebih disukai error dari model ada yang positif dan negatif dengan nilai harapan nol. Di kasus ini, jika disekat, ada wilayah yang nilai harapan positif dan negatif sehingga tidak sesuai harapan.

Data Menurut Negara

Asal Negara Produsen Amerika

AutoAmerika <- Auto[Auto$origin==1,] %>%  select(mpg, horsepower, origin)
datatable(AutoAmerika)

Visualisasi Hubungan mpg terhadap horsepower

ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="dodgerblue2") +
  theme_bw() +
   labs(y="Mile per Gallon", x="Horse Power", title = "Scatter Plot `mpg` dan `horsepower`",
        subtitle = "Amerika")

Asal Negara Produsen Eropa

AutoEropa <- Auto[Auto$origin==2,] %>%  select(mpg, horsepower, origin)
datatable(AutoEropa)

Visualisasi Hubungan mpg terhadap horsepower

ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="dodgerblue2") +
  theme_bw() +
   labs(y="Mile per Gallon", x="Horse Power", title = "Scatter Plot `mpg` dan `horsepower`",
        subtitle = "Eropa")

Deteksi Outlier

# Mengidentifikasi outlier pada kolom mpg dan horsepower
outliers_mpg_er <- AutoEropa %>% identify_outliers(mpg)
outliers_hp_er <- AutoEropa %>% identify_outliers(horsepower)

# Menggabungkan outlier dari kedua kolom
outliers <- rbind(outliers_mpg_er, outliers_hp_er)

# Membuat data frame baru tanpa outlier
AutoEropa_clean <- AutoEropa[!(AutoEropa$mpg %in% outliers$mpg & AutoEropa$horsepower %in% outliers$horsepower), ]
ggplot(AutoEropa_clean,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="dodgerblue2") +
  theme_bw() +
   labs(y="Mile per Gallon", x="Horse Power", title = "Scatter Plot `mpg` dan `horsepower`",
        subtitle = "Eropa")

Berdasarkan plot di atas, terlihat bahwa tidak terdapat perubahan yang signifikan setelah dilakukan penanganan outlier sehingga data yang digunakan untuk analisis tetap berdasarkan data awal.

Asal Negara Produsen Jepang

AutoJepang <- Auto[Auto$origin==3,] %>%  select(mpg, horsepower, origin)
datatable(AutoJepang)

Visualisasi Hubungan mpg terhadap horsepower

ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="dodgerblue2") +
  theme_bw() +
   labs(y="Mile per Gallon", x="Horse Power", title = "Scatter Plot `mpg` dan `horsepower`",
        subtitle = "Jepang")

Deteksi Outlier

# Mengidentifikasi outlier pada kolom mpg dan horsepower
outliers_mpg_jp <- AutoJepang %>% identify_outliers(mpg)
outliers_hp_jp <- AutoJepang %>% identify_outliers(horsepower)

# Menggabungkan outlier dari kedua kolom
outliers <- rbind(outliers_mpg_jp, outliers_hp_jp)

# Membuat data frame baru tanpa outlier
AutoJepang_clean <- AutoJepang[!(AutoJepang$mpg %in% outliers$mpg & AutoJepang$horsepower %in% outliers$horsepower), ]
ggplot(AutoJepang_clean,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="dodgerblue2") +
  theme_bw() +
   labs(y="Mile per Gallon", x="Horse Power", title = "Scatter Plot `mpg` dan `horsepower`",
        subtitle = "Eropa")

Berdasarkan plot di atas, terlihat bahwa tidak terdapat perubahan yang signifikan setelah dilakukan penanganan outlier sehingga data yang digunakan untuk analisis tetap berdasarkan data awal.

Pemodelan Data Auto

Polynomial Regression

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Polynomial Regression.

set.seed(123)

cv.poly1 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly1 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,1,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "RMSE"
      return(metric)
      }
  )
  RegLin<- colMeans(metric_poly1) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(RegLin) 
}

cv.poly2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly2 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,2,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse 
      names(metric) <- "RMSE"
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly3 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,3,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "RMSE"
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly4 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,4,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "RMSE"
      return(metric)
      }
  )
  poly4<- colMeans(metric_poly4) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly4) 
}
cvpoly <- rbind (cv.poly1(Auto),cv.poly2(Auto),cv.poly3(Auto),cv.poly4(Auto))
cvpoly
##          RMSE
## [1,] 4.889086
## [2,] 4.351129
## [3,] 4.355430
## [4,] 4.367009

Secara empiris jika dilihat dari RMSE terlihat bahwa polinomial derajat 2 memberikan nilai RMSE yang terkecil sehingga dipilih Polynomial Regression derajat 2.

Secara Visual

poly2 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  ggtitle("Polynomial Regression Derajat 2 (RMSE=4.35)") +
  theme_bw()
poly2

Piecewise Constant

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Piecewise Constant.

set.seed(123)

cv.pc <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 2:10

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[x$in_id,]
training$horsepower <- cut(training$horsepower,i)

mod <- lm(mpg ~ horsepower,
         data=training)

labs_x <- levels(mod$model[,2])

labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
                  upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))


testing <- dataset[-x$in_id,]

horsepower_new <- cut(testing$horsepower,c(labs_x_breaks[1,1],labs_x_breaks[,2]))

pred <- predict(mod,
               newdata=list(horsepower=horsepower_new))
truth <- testing$mpg

data_eval <- na.omit(data.frame(truth,pred))

rmse <- mlr3measures::rmse(truth = data_eval$truth,
                           response = data_eval$pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
return(best_tangga)
}
cv.pc(Auto)
##   breaks     rmse
## 1      2 6.147626
## 2      3 5.775792
## 3      4 4.985270
## 4      5 4.712623
## 5      6 4.644383
## 6      7 4.554116
## 7      8 4.437597
## 8      9 4.549668
## 9     10 4.589172

Secara empiris berdasarkan nilai RMSE, banyaknya interval yang terbaik adalah 8 interval sehingga dipilih Piecewise Constant dengan 8 interval.

Secara Visual

pc8 <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,8), 
               lty = 1, col = "blue",se = F)+
   ggtitle("Piecewise Constant 8 Interval (RMSE=4.44)") +
  theme_bw()
pc8

Natural Cubic Splines

Secara Empiris dengan Cross Validation

Menentukan Banyaknya Fungsi Basis, sebagai berikut:

set.seed(123)
cross.val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 2:10
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- map_dfr(cross.val$splits,
                  function(x){
                  mod <- lm(mpg ~ ns(horsepower,df=i),data=Auto[x$in_id,])
                  pred <- predict(mod,newdata=Auto[-x$in_id,])
                  truth <- Auto[-x$in_id,]$mpg
                  rmse <- mlr3measures::rmse(truth = truth,response = pred)
                  metric <- rmse
                  names(metric) <- "rmse"
                  return(metric)
                }
      )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-best.spline3 %>% slice_min(rmse) #berdasarkan rmse

basis
##   df     rmse
## 1 10 4.260316

Banyaknya fungsi basis yang terpilih adalah 10. Selanjutnya adalah menentukan knot berdasarkan banyaknya fungsi basis yang terpilih.

df_10<-attr(ns(Auto$horsepower, df=10),"knots")
df_10
##   10%   20%   30%   40%   50%   60%   70%   80%   90% 
##  67.0  72.0  80.0  88.0  93.5 100.0 110.0 140.0 157.7

Cross Validation

set.seed(123)

cross.val <- vfold_cv(Auto,v=10,strata = "mpg")

metric.spline3.ns4 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(80, 120, 160, 200)),
         data=Auto[x$in_id,])
pred <- predict(mod,
               newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns4 <- colMeans(metric.spline3.ns4)# menghitung rata-rata untuk 10 folds



metric.spline3.ns9 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(67,72,80,88,93.5,100,110,140,157.7)),
         data=Auto[x$in_id,])
pred <- predict(mod,
               newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns9 <- colMeans(metric.spline3.ns9)# menghitung rata-rata untuk 10 folds



nilai.cv.ns <- rbind("NCS4" = mean.metric.spline3.ns4,
                     "NCS9" = mean.metric.spline3.ns9)

nama.model.ncs <- c("4 - NCS Knots 80, 120, 160, 200",
                    "9 - NCS (df=10)")

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs
##                                Model     rmse
## NCS4 4 - NCS Knots 80, 120, 160, 200 4.359156
## NCS9                 9 - NCS (df=10) 4.261893

Secara empiris berdasarkan nilai RMSE, banyaknya knot terbaik adalah 9 knots, yaitu knot 67,72,80,88,93.5,100,110,140,157.7 sehingga dipilih Natural Cubic Splines dengan 9 knots.

Secara Visual

nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(67,72,80,88,93.5,100,110,140,157.7)), col = "blue", se = F)+
  ggtitle("Natural Cubic Splines 9 Knots (RMSE=4.26)") +
  theme_bw()
nc9

Perbandingan Model Terbaik

Secara Empiris

poly2.a <- cv.poly2(Auto)
pc.a <- cv.pc(Auto)[1,-1]
ncs.a <- model.ncs[2,-1]

# Membuat dataframe
df <- data.frame("Metode" = c("Polinomial Derajat 2", "Piecewise Constant 8", "NCS 9 Knots"), "RMSE" = c(poly2.a, pc.a, ncs.a))
df
##                 Metode     RMSE
## 1 Polinomial Derajat 2 4.351129
## 2 Piecewise Constant 8 6.147626
## 3          NCS 9 Knots 4.261893

Berdasarkan hasil analisis, model terbaik untuk data Auto adalah Natural Cubic Splines dengan 9 knots.

Secara Visual

# Membuat plot gabungan
combined_plot <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  ggtitle("Comparison of Best Model") +
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", 
              formula = y~poly(x,2,raw=T), aes(col = "Polinomial Derajat 2"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~cut(x,8), aes(col = "Piecewise Constant 8"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = c(67,72,80,88,93.5,100,110,140,157.7)), 
              aes(col = "NCS 9 Knots"), se = F, linetype = "solid") +
  labs(color="Model")+
  scale_color_manual(values=c("Polinomial Derajat 2"="blue", "Piecewise Constant 8"="red","NCS 9 Knots"= "green")) +
  theme_bw()
combined_plot

Secara umum, model Natural Cubic Splines 9 Knots tampaknya paling akurat dalam memodelkan hubungan ini. Model ini tampaknya paling akurat dalam memprediksi mpg berdasarkan horsepower. Garis biru mewakili model ini dan tampaknya paling dekat dengan titik data, menunjukkan bahwa model ini memiliki kesalahan prediksi yang paling rendah.

Pemodelan Data Amerika

Polynomial Regression

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan polynomial regression.

cvpoly <- rbind (cv.poly1(AutoAmerika),cv.poly2(AutoAmerika),cv.poly3(AutoAmerika),cv.poly4(AutoAmerika))
cvpoly
##          RMSE
## [1,] 4.237175
## [2,] 3.799433
## [3,] 3.804542
## [4,] 3.823584

Secara empiris jika dilihat dari nilai RMSE terlihat bahwa polinomial derajat 2 memberikan nilai RMSE yang terkecil sehingga dipilih Polynomial Regression derajat 2.

Secara Visual

poly2 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "coral2",se = F)+
  theme_bw()+
  labs(title = "Polynomial Derajat 2 (RMSE=3.80)",
        subtitle = "Amerika")
poly2

Piecewise Constant

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Piecewise Constant.

cv.pc(AutoAmerika)
##   breaks     rmse
## 1      2 4.997476
## 2      3 4.759908
## 3      4 4.139362
## 4      5 3.966889
## 5      6 4.222039
## 6      7 4.150320
## 7      8 3.945534
## 8      9 3.737522
## 9     10 3.816464

Secara empiris berdasarkan RMSE, banyaknya interval yang terbaik adalah 9 sehingga model yang dipilih adalah Piecewise Constant dengan 9 Interval

pc9 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,9), 
               lty = 1, col = "coral2",se = F)+
  theme_bw()+
  labs(title = "Piecewise Constant 9 Interval (RMSE=3.74)",
        subtitle = "Amerika")
pc9

Natural Cubic Spline

Secara Empiris dengan Cross Validation

Menentukan Banyaknya Fungsi Basis, sebagai berikut

set.seed(123)
cross.val <- vfold_cv(AutoAmerika,v=10,strata = "mpg")
df <- 2:10
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- map_dfr(cross.val$splits,
                  function(x){
                  mod <- lm(mpg ~ ns(horsepower,df=i),data=AutoAmerika[x$in_id,])
                  pred <- predict(mod,newdata=AutoAmerika[-x$in_id,])
                  truth <- AutoAmerika[-x$in_id,]$mpg
                  rmse <- mlr3measures::rmse(truth = truth,response = pred)
                  metric <- rmse
                  names(metric) <- "rmse"
                  return(metric)
                }
      )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-best.spline3 %>% slice_min(rmse) #berdasarkan rmse

basis
##   df     rmse
## 1  7 3.769022

Diperoleh df atau banyaknya fungsi basis terpilih sebesar 7. Selanjutnya menentukan knot berdasarkan banyaknya fungsi basis yang terpilih.

df_7 <-attr(ns(AutoAmerika$horsepower, df=7),"knots")
df_7
## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429% 
##  82.85714  90.00000 100.00000 110.00000 145.00000 165.00000

Cross Validation

set.seed(123)

cross.val <- vfold_cv(AutoAmerika,v=10,strata = "mpg")

metric.spline3.ns4 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(80, 120, 160, 200)),
         data=AutoAmerika[x$in_id,])
pred <- predict(mod,
               newdata=AutoAmerika[-x$in_id,])
truth <- AutoAmerika[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns4 <- colMeans(metric.spline3.ns4)# menghitung rata-rata untuk 10 folds



metric.spline3.ns7 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(82.85714,90,100,110,145,165)),
         data=AutoAmerika[x$in_id,])
pred <- predict(mod,
               newdata=AutoAmerika[-x$in_id,])
truth <- AutoAmerika[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns7 <- colMeans(metric.spline3.ns7)# menghitung rata-rata untuk 10folds



nilai.cv.ns <- rbind("NCS4" = mean.metric.spline3.ns4,
                     "NCS6" = mean.metric.spline3.ns7)

nama.model.ncs <- c("4 - NCS Knots 80, 120, 160, 200",
                    "6 - NCS (df=7)")

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs
##                                Model     rmse
## NCS4 4 - NCS Knots 80, 120, 160, 200 3.823310
## NCS6                  6 - NCS (df=7) 3.761205

Secara empiris berdasarkan RMSE, banyaknya knot terbaik adalah 6 knots sehingga model yang dipilih adalah Natural Cubic Splines dengan 6 Knots.

Secara Visual

nc6 <- ggplot(AutoAmerika, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(82.85714,90,100,110,145,165)), col = "coral2", se = F)+
  theme_bw()+
   labs(title = "Natural Cubic Splines 6 Knots (RMSE=3.76)",
        subtitle = "Amerika")
nc6

Perbandingan Model Terbaik

Secara Empiris

poly2.b <- cv.poly2(AutoAmerika)
pc.b <- cv.pc(AutoAmerika)[1,-1]
ncs.b <- model.ncs[2,-1]

# Membuat dataframe
df <- data.frame("Metode" = c("Polinomial Derajat 2", "Piecewise Constant 9", "NCS 6 Knots"), "RMSE" = c(poly2.b, pc.b, ncs.b))
df
##                 Metode     RMSE
## 1 Polinomial Derajat 2 3.799433
## 2 Piecewise Constant 9 4.997476
## 3          NCS 6 Knots 3.761205

Berdasarkan hasil analisis, model terbaik untuk data Auto Amerika adalah Piecewise Constant dengan 9 interval.

Secara Visual

# Membuat plot gabungan
combined_plot <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  ggtitle("Comparison of Best Model") +
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", 
              formula = y~poly(x,2,raw=T), aes(col = "Polinomial Derajat 2"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~cut(x,9), aes(col = "Piecewise Constant 9"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = c(82.85714,90,100,110,145,165)), 
              aes(col = "NCS 6 Knots"), se = F, linetype = "solid") +
  labs(color="Model")+
  scale_color_manual(values=c("Polinomial Derajat 2"="blue", "Piecewise Constant 9"="red","NCS 6 Knots"= "green")) +
  theme_bw()
combined_plot

Secara keseluruhan, plot ini menunjukkan bahwa ada hubungan non-linear antara horsepower dan mpg. Model Piecewise Constant 9, sebagai model terbaik berdasarkan hasil Cross Validation, membagi data menjadi 9 segmen dan mengasumsikan fungsi konstan di setiap segmen. Ini memungkinkan model untuk menangkap pola non-linear dalam data tanpa harus memperkirakan fungsi yang rumit.

Pemodelan Data Eropa

Polynomial Regression

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Polynomial Regression.

set.seed(123)


cv.poly2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly2 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,2,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse 
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly3 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,3,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly4 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,4,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse 
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly4<- colMeans(metric_poly4) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly4) 
}
cvpoly <- rbind (cv.poly2(AutoEropa),cv.poly3(AutoEropa),cv.poly4(AutoEropa))
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
## The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
## The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
cvpoly
##          rmse
## [1,] 4.805676
## [2,] 4.832309
## [3,] 4.878846

Secara empiris jika dilihat dari nilai RMSE terlihat bahwa polinomial derajat 2 memberikan nilai RMSE yang terkecil sehingga dipilih Polynomial Regression derajat 2.

Secara Visual

poly2 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Polynomial Derajat 2 (RMSE=4.805)",
        subtitle = "Eropa")
poly2

Piecewise Constant

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Piecewise Constant.

set.seed(123)

cv.pc2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 2:6

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[x$in_id,]
training$horsepower <- cut(training$horsepower,i)

mod <- lm(mpg ~ horsepower,
         data=training)

labs_x <- levels(mod$model[,2])

labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
                  upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))


testing <- dataset[-x$in_id,]

horsepower_new <- cut(testing$horsepower,c(labs_x_breaks[1,1],labs_x_breaks[,2]))

pred <- predict(mod,
               newdata=list(horsepower=horsepower_new))
truth <- testing$mpg

data_eval <- na.omit(data.frame(truth,pred))

rmse <- mlr3measures::rmse(truth = data_eval$truth,
                           response = data_eval$pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
return(best_tangga)
}
cv.pc2(AutoEropa)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
##   breaks     rmse
## 1      2 5.462277
## 2      3 5.395409
## 3      4 5.226264
## 4      5 4.963259
## 5      6 5.133934

Secara empiris berdasarkan nilai RMSE, banyaknya interval yang terbaik adalah 5 interval sehingga model yang dipilih adalah Piecewise Contants dengan 5 interval.

pc5 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,5), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Piecewise Constant 5 (RMSE=4.963)",
        subtitle = "Eropa")
pc5

Natural Cubic Spline

Secara Empiris dengan Cross Validation

Menentukan Banyaknya Fungsi Basis

set.seed(123)
cross.val <- vfold_cv(AutoEropa,v=10,strata = "mpg")
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
df <- 2:10
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- map_dfr(cross.val$splits,
                  function(x){
                  mod <- lm(mpg ~ ns(horsepower,df=i),data=AutoEropa[x$in_id,])
                  pred <- predict(mod,newdata=AutoEropa[-x$in_id,])
                  truth <- AutoEropa[-x$in_id,]$mpg
                  rmse <- mlr3measures::rmse(truth = truth,response = pred)
                  metric <- rmse
                  names(metric) <- "rmse"
                  return(metric)
                }
      )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-data.frame(best.spline3 %>% slice_min(rmse)) #berdasarkan rmse

basis
##   df     rmse
## 1  2 4.807144

Diperoleh basis yang terbaik adalah 2, selanjutnya menentukan knot berdasarkan banyaknya fungsi basis yang terpilih

df_2<-attr(ns(AutoEropa$horsepower, df=2),"knots")
df_2
##  50% 
## 76.5

Cross Validation

metric.spline3.ns1 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(76.5)),
         data=AutoEropa[x$in_id,])
pred <- predict(mod,
               newdata=AutoEropa[-x$in_id,])
truth <- AutoEropa[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns1 <- colMeans(metric.spline3.ns1)# menghitung rata-rata untuk 10 folds



metric.spline3.ns2 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(88, 105)),
         data=AutoEropa[x$in_id,])
pred <- predict(mod,
               newdata=AutoEropa[-x$in_id,])
truth <- AutoEropa[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns2 <- colMeans(metric.spline3.ns2)# menghitung rata-rata untuk 10 folds



nilai.cv.ns <- rbind("NCS1" = mean.metric.spline3.ns1,
                     "NCS2" = mean.metric.spline3.ns2)

nama.model.ncs <- c("1 - NCS Knots 76.5",
                    "2 - NCS Knots 88, 105")

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs
##                      Model     rmse
## NCS1    1 - NCS Knots 76.5 4.806982
## NCS2 2 - NCS Knots 88, 105 4.825955

Secara empiris berdasarkan nilai RMSE, banyaknya knot terbaik adalah 1 knot sehingga model yang dipilih adalah Natural Cubic Splines dengan 1 knot.

Secara Visual

nc1 <- ggplot(AutoEropa, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(76.5)), col = "blue", se = F) +
  theme_bw()+
  labs(title = "Natural Cubic Splines 1 Knot (RMSE=4.806)",
        subtitle = "Eropa")
nc1

Perbandingan Model Terbaik

Secara Empiris

poly2.c <- cv.poly2(AutoEropa)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
pc.c <- cv.pc2(AutoEropa)[1,-1]
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
ncs.c <- model.ncs[1,-1]

# Membuat dataframe
df <- data.frame("Metode" = c("Polinomial Derajat 2", "Piecewise Constant 5", "NCS 1 Knot"), "RMSE" = c(poly2.c, pc.c, ncs.c))
df
##                 Metode     RMSE
## 1 Polinomial Derajat 2 4.805676
## 2 Piecewise Constant 5 5.462277
## 3           NCS 1 Knot 4.806982

Berdasarkan hasil analisis, model terbaik untuk data Auto Eropa adalah Polynomial Regression derajat 2.

Secara Visual

# Membuat plot gabungan
combined_plot <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  ggtitle("Comparison of Best Model") +
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", 
              formula = y~poly(x,2,raw=T), aes(col = "Polinomial Derajat 2"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~cut(x,5), aes(col = "Piecewise Constant 5"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = c(76.5)), 
              aes(col = "NCS 1 Knot"), se = F, linetype = "solid") +
  labs(color="Model")+
  scale_color_manual(values=c("Polinomial Derajat 2"="blue", "Piecewise Constant 5"="red","NCS 1 Knot"= "green")) +
  theme_bw()
combined_plot

Secara keseluruhan, plot ini menunjukkan bahwa ada hubungan non-linear antara tenaga kuda dan efisiensi bahan bakar, dan model Polynomial Derajat 2 tampaknya paling akurat dalam memodelkan hubungan ini.

Berdasarkan model Polynomial Derajat 2, dapat disimpulkan bahwa ada hubungan non-linear (horsepower) dan efisiensi bahan bakar (mpg). Model ini tampaknya memberikan perkiraan yang baik tentang hubungan antara tenaga kuda dan efisiensi bahan bakar, meskipun ada beberapa penyimpangan dari titik data.

Pemodelan Data Jepang

Polynomial Regression

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Polynomial Regression.

set.seed(123)

cv.poly2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly2 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,2,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse 
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly3 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,3,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  metric_poly4 <- map_dfr(cross_val$splits,
    function(x){  
      mod <- lm(mpg ~ poly(horsepower,4,raw = T),
         data=dataset[x$in_id,]) 
      pred <- predict(mod,
               newdata=dataset[-x$in_id,]) 
      truth <- dataset[-x$in_id,]$mpg
      rmse <- mlr3measures::rmse(truth = truth,response = pred)
      metric <- rmse
      names(metric) <- "rmse"
      return(metric)
      }
  )
  poly4<- colMeans(metric_poly4) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly4) 
}
cvpoly <- rbind (cv.poly2(AutoJepang),cv.poly3(AutoJepang),cv.poly4(AutoJepang))
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
## The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
## The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
cvpoly
##          rmse
## [1,] 4.522705
## [2,] 3.962162
## [3,] 4.380986

Secara empiris jika dilihat dari nilai RMSE terlihat bahwa polinomial derajat 2 memberikan nilai rmse yang terkecil sehingga dipilih Polynomial Regression derajat 3.

Secara Visual

poly3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Polinomial Derajat 3 (RMSE=3.96)",
        subtitle = "Jepang")
poly3

Piecewise Constant

Secara Empiris dengan Cross Validation

Cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan Piecewise Constant.

set.seed(123)

cv.pc2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")

breaks <- 2:5

best_tangga <- map_dfr(breaks, function(i){

metric_tangga <- map_dfr(cross_val$splits,
    function(x){

training <- dataset[x$in_id,]
training$horsepower <- cut(training$horsepower,i)

mod <- lm(mpg ~ horsepower,
         data=training)

labs_x <- levels(mod$model[,2])

labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
                  upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))


testing <- dataset[-x$in_id,]

horsepower_new <- cut(testing$horsepower,c(labs_x_breaks[1,1],labs_x_breaks[,2]))

pred <- predict(mod,
               newdata=list(horsepower=horsepower_new))
truth <- testing$mpg

data_eval <- na.omit(data.frame(truth,pred))

rmse <- mlr3measures::rmse(truth = data_eval$truth,
                           response = data_eval$pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
}
)
best_tangga <- cbind(breaks=breaks,best_tangga) # menampilkan hasil all breaks
return(best_tangga)
}
cv.pc2(AutoJepang)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
##   breaks     rmse
## 1      2 4.306154
## 2      3 4.201949
## 3      4 4.072662
## 4      5 4.226607

Secara empiris berdasarkan nilai RMSE, banyaknya interval yang terbaik adalah 4 interval sehingga model yang dipilih adalah Piecewiseconstant dengan 4 interval,

Secara Visual

pc4 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,4), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Picewise Constant 4 (RMSE=4.07)",
        subtitle = "Jepang")
pc4

Natural Cubic Spline

Secara Empiris dengan Cross Validation

Menentukan Banyaknya Fungsi Basis

set.seed(123)
cross.val <- vfold_cv(AutoJepang,v=10,strata = "mpg")
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
df <- 2:10
best.spline3 <- map_dfr(df, function(i){
metric.spline3 <- map_dfr(cross.val$splits,
                  function(x){
                  mod <- lm(mpg ~ ns(horsepower,df=i),data=AutoJepang[x$in_id,])
                  pred <- predict(mod,newdata=AutoJepang[-x$in_id,])
                  truth <- AutoJepang[-x$in_id,]$mpg
                  rmse <- mlr3measures::rmse(truth = truth,response = pred)
                  metric <- rmse
                  names(metric) <- "rmse"
                  return(metric)
                }
      )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
)


best.spline3 <- cbind(df=df,best.spline3)
basis<-best.spline3 %>% slice_min(rmse) #berdasarkan rmse

basis
##   df     rmse
## 1  3 4.143445

Diperoleh basis yang terbaik adalah 3, selanjutnya menentukan knot berdasarkan banyaknya fungsi basis yang terpilih.

df_3<-attr(ns(AutoJepang$horsepower, df=3),"knots")
df_3
## 33.33333% 66.66667% 
##        68        90

Cross Validation

set.seed(123)

cross.val <- vfold_cv(AutoJepang,v=10,strata = "mpg")
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
metric.spline3.ns2 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(68,90)),
         data=AutoJepang[x$in_id,])
pred <- predict(mod,
               newdata=AutoJepang[-x$in_id,])
truth <- AutoJepang[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns2 <- colMeans(metric.spline3.ns2) # menghitung rata-rata untuk 10 folds



metric.spline3.ns3 <- map_dfr(cross.val$splits,
    function(x){
mod <- lm(mpg ~ ns(horsepower, knots = c(60, 88, 105)),
         data=AutoJepang[x$in_id,])
pred <- predict(mod,
               newdata=AutoJepang[-x$in_id,])
truth <- AutoJepang[-x$in_id,]$mpg

rmse <- mlr3measures::rmse(truth = truth,
                           response = pred
                           )
metric <- rmse
names(metric) <- "rmse"
return(metric)
}
)
mean.metric.spline3.ns3 <- colMeans(metric.spline3.ns3) # menghitung rata-rata untuk 10 folds



nilai.cv.ns <- rbind("NCS2" = mean.metric.spline3.ns2,
                     "NCS3" = mean.metric.spline3.ns3)

nama.model.ncs <- c("2 - NCS Knots 68,90",
                    "3 - NCS Knots 60, 88, 105")

model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)

model.ncs
##                          Model     rmse
## NCS2       2 - NCS Knots 68,90 4.140871
## NCS3 3 - NCS Knots 60, 88, 105 4.217066

Secara empiris berdasarkan nilai RMSE, banyaknya knots terbaik adalah 2 knots sehingga model yang terpilih adalah Natural Cubic Splines dengan 2 knots.

Secara Visual

nc2 <- ggplot(AutoJepang, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(68,90)), col = "blue", se = F) +
  theme_bw()+
    labs(title = "Natural Cubic Splines 2 Knots (RMSE=4.14)",
        subtitle = "Jepang")
nc2

Perbandingan Model Terbaik

Secara Empiris

poly2.d <- cv.poly3(AutoJepang)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
pc.d <- cv.pc2(AutoJepang)[1,-1]
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
ncs.d <- model.ncs[1,-1]

# Membuat dataframe
df <- data.frame("Metode" = c("Polinomial Derajat 2", "Piecewise Constant 4", "NCS 2 Knots"), "RMSE" = c(poly2.d, pc.d, ncs.d))
df
##                 Metode     RMSE
## 1 Polinomial Derajat 2 3.962162
## 2 Piecewise Constant 4 4.306154
## 3          NCS 2 Knots 4.140871

Berdasarkan hasil analisis, model terbaik untuk data Auto Eropa adalah Polynomial Regression derajat 2.

Secara Visual

# Membuat plot gabungan
combined_plot <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  ggtitle("Comparison of Best Model") +
  geom_point(alpha = 0.8, color="black") + 
  stat_smooth(method = "lm", 
              formula = y~poly(x,3,raw=T), aes(col = "Polinomial Derajat 3"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~cut(x,4), aes(col = "Piecewise Constant 4"), se = F, 
              linetype = "solid") +
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = c(68,90)), 
              aes(col = "NCS 2 Knots"), se = F, linetype = "solid") +
  labs(color="Model")+
  scale_color_manual(values=c("Polinomial Derajat 3"="blue", "Piecewise Constant 4"="red","NCS 2 Knots"= "green")) +
  theme_bw()
combined_plot

Polynomial Derajat 3: Model ini, yang diwakili oleh garis biru, tampaknya paling akurat dalam memprediksi efisiensi bahan bakar berdasarkan tenaga kuda. Garis ini tampaknya paling dekat dengan titik data, menunjukkan bahwa model ini memiliki kesalahan prediksi yang paling rendah.

Secara keseluruhan, plot ini menunjukkan bahwa ada hubungan non-linear antara horsepower dan mpg.

Kesimpulan

Berdasarkan hasil modeling non linier, pemilihan tuning parameter menggunakan 10-fold cross validation, dan evaluasi diperoleh model terbaik untuk masing-masing data sebagai berikut:

mod.auto <- model.ncs[2,-1]
mod.Amerika <- cv.pc(AutoAmerika)[1,-1]
mod.Eropa <- cv.poly2(AutoEropa)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
mod.Jepang <- cv.poly3(AutoJepang)
## Warning: The number of observations in each quantile is below the recommended threshold of 20.
## • Stratification will use 3 breaks instead.
df <- data.frame("Metode" = c("NCS 9 Knots", "Piecewise Constant 9", "Polynomial Derajat 2", "Polynomial Derajat 3"), "RMSE" = c(mod.auto, mod.Amerika, mod.Eropa, mod.Jepang))
df$Negara <- c("Global", "Amerika", "Eropa", "Jepang")
df
##                 Metode     RMSE  Negara
## 1          NCS 9 Knots 4.217066  Global
## 2 Piecewise Constant 9 4.997476 Amerika
## 3 Polynomial Derajat 2 4.805676   Eropa
## 4 Polynomial Derajat 3 3.962162  Jepang

Plot Model Terbaik Berdasarkan Asal Negara

Hasil modeling dapat disajikan dalam diagram pencar dan kurva prediksi sebagai berikut:

pc9 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,9), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Piecewise Constant 9 Interval",
        subtitle = "Amerika")

poly2 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Polynomial Derajat 2",
        subtitle = "Eropa")

poly3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.8, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Polinomial Derajat 3",
        subtitle = "Jepang")

plot_grid(pc9, poly2, poly3)

Berdasarkan model terbaik, ada beberapa kesimpulan yang diperoleh yaitu sebagai berikut:

  1. Model non linear terbaik yang terpilih untuk data Auto adalah Natural Cubic Splines dengan 9 knots, sedangkan untuk data berdasarkan asal negara model non linear terbaik yang terpilih adalah: Amerika (Piecewise Constant dengan 9 Interval), Eropa (Polynomial Regression derajat 2), dan Jepang (Polynomial Regression derajat 2).

  2. Secara umum, model terbaik yang terpilih untuk semua data baik itu secara keseluruhan ataupun berdasarkan asal negara menunjukkan bahwa peningkatan tenaga kuda tidak selalu berarti penurunan efisiensi bahan bakar yang signifikan. Sebaliknya, ada titik di mana peningkatan tenaga kuda mulai memiliki dampak negatif yang lebih besar pada efisiensi bahan bakar. Namun, perlu diingat bahwa model ini hanyalah perkiraan dan mungkin tidak selalu berlaku untuk setiap kendaraan atau situasi.

Saran

Saran yang dapat diberikan berdasarkan analisis adalah:

  1. Terdapat hubungan non-linear antara kekuatan mesin horsepower dan efisiensi bahan bakar mpg. Secara umum, peningkatan kekuatan mesin dapat mengakibatkan penurunan efisiensi bahan bakar. Oleh karena itu, jika efisiensi bahan bakar adalah prioritas, pertimbangkan untuk memilih kendaraan dengan tenaga kuda yang lebih rendah.

  2. Model yang digunakan dalam analisis ini memberikan perkiraan yang baik tentang hubungan antara kekuatan mesin horsepower dan efisiensi bahan bakar mpg. Namun, model ini hanyalah perkiraan dan mungkin tidak selalu berlaku untuk setiap kendaraan atau situasi. Oleh karena itu, selalu baik untuk memvalidasi model dengan data tes atau data baru sebelum membuat keputusan berdasarkan prediksi model.

  3. Pertimbangkan Model yang Lebih Sederhana: Model yang lebih sederhana bisa lebih berguna jika mereka cukup akurat dan lebih mudah diinterpretasikan.

Referensi

Sartono, B. (November 1, 2020). Moving Beyond Linearity. Retrieved from: https://rpubs.com/bagusco/taklinear

Amalia, R. (April 23, 2021). Moving Beyond Linearity. Retrieved from: https://rpubs.com/reniamelia/MovingBeyondLinearity

Setiabudi, N.A. (November 3, 2021). Moving Beyond Linearity. Retrieved from: https://rpubs.com/nurandi/moving-beyond-linearity

Kuhn, M., & Johnson, K. (2013). Applied Predictive Modeling. New York, NY: Springer.

Mahasiswa Pascasarjana Statistika dan Sains Data, IPB University,