Pola Hubungan Tak Linear (Beyond Linearity)

A. Pendahuluan

Soal

Lakukan cross-validation untuk menghasilkan pemodelan optimal data Horse Power terhadap Mile per Gallon pada dataset Auto menggunakan 3 metode berikut, selanjutnya melakukan komparasi antar model.
* Regresi polynomial

  • Fungsi tangga (Piecewise constant)

  • Natural cubic splines

Lakukan hal di atas untuk masing-masing subset data berdasarkan asal negara produsennya (Amerika, Japan, dan Eropa).

Plot model terbaik dalam satu frame kemudian berikan ulasannya.

A.1 Eksplorasi Data

Deskripsi 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 Skala Pengukuran
1 Cylinder Silinder yang digunakan pada kendaraan Nominal
2 Displacement Besarnya perpindahan Interval
3 Horsepower Tenaga kendaraan Interval
4 Weight Berat kendaraan Interval
5 Acceleration Percepatan kendaraan Interval
6 Year Model Tahun keluaran kendaraan Nominal
7 Origin Asal kendaraan Nominal 
library(pacman)
p_load("tidyverse","ISLR","ggplot2","splines", "caret", 
       "rsample", "mlr3measures", "cowplot", "rmarkdown",
       "DT")

# Import data
data(Auto)
dataall <- Auto
glimpse(dataall)
## Rows: 392
## Columns: 9
## $ mpg          <dbl> 18, 15, 18, 16, 17, 15, 14, 14, 14, 15, 15, 14, 15, 14, 2~
## $ cylinders    <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 6, 6, 6, 4, ~
## $ displacement <dbl> 307, 350, 318, 304, 302, 429, 454, 440, 455, 390, 383, 34~
## $ horsepower   <dbl> 130, 165, 150, 150, 140, 198, 220, 215, 225, 190, 170, 16~
## $ weight       <dbl> 3504, 3693, 3436, 3433, 3449, 4341, 4354, 4312, 4425, 385~
## $ acceleration <dbl> 12.0, 11.5, 11.0, 12.0, 10.5, 10.0, 9.0, 8.5, 10.0, 8.5, ~
## $ year         <dbl> 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 7~
## $ origin       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 3, ~
## $ name         <fct> chevrolet chevelle malibu, buick skylark 320, plymouth sa~
datatable(dataall,class = 'cell-border stripe', caption = htmltools::tags$caption(
    style = 'caption-side: bottom; text-align: center;',
    'Table 1: ', htmltools::em('Tabel keseluruhan dataset Auto')
    )
    )

A.2 Melihat dugaan hubungan non-linear antara peubah prediktor terhadap peubah respon

Pada analisis ini dipilih 3 hubungan antara peubah prediktor terhadap peubah respon yakni:

  1. weight vs mpg

  2. acceleration vs mpg

  3. displacement vs mpg

Menyimpan data yang terdiri dari peubah mpg, weight, acceleration, displacement, origin ke variabel baru bernama data3

data3 <- dataall %>% select(mpg, weight, acceleration, displacement, origin)
# Plot weight vs mpg
A <- ggplot(dataall,aes(x=weight, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
              formula = y~x,lty = 1,
              col = "blue",se = F) +
  xlab("weight") +
  ylab("miles per gallon") +
  ggtitle("1. Regresi Linear Weight VS MPG") +
  theme_light(base_size = 8)

# Plot acceleration vs mpg
B <- ggplot(dataall,aes(x=acceleration, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
              formula = y~x,lty = 1,
              col = "orange",se = F) +
  xlab("acceleration") +
  ylab("miles per gallon") +
  ggtitle("2. Regresi Linear Acceleration VS MPG") +
  theme_light(base_size = 8)

# Plot displacement vs mpg
C <- ggplot(dataall,aes(x=displacement, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
              formula = y~x,lty = 1,
              col = 3,se = F) +
  xlab("displacement") +
  ylab("miles per gallon") +
  ggtitle("3. Regresi Linear Displacement VS MPG") +
  theme_light(base_size = 8)

# install.packages("ggpubr")
ggpubr::ggarrange(A, B, C, ncol = 2, nrow = 2)

Berdasarkan visualisasi di atas, dari ketiga hubungan antar peubah terlihat seperti hubungan linear yang digambarkan oleh garis lurus berwarna pada setiap gambar. Namun apabila diperhatikan lebih detail, garis lurus tersebut nampak kurang sesuai (mengikuti) kondisi data yang ada.

  • Kasus 1: weight vs mpg

    Bila hubungan tersebut mengikuti regresi linear, untuk nilai weight (berat kendaraan) yang besar akan diperoleh nilai mpg (jarak yang dapat ditempuh per gallon bahan bakar) negatif. Hal ini kurang relevan dengan konteks miles per gallon.

  • Kasus 2: acceleration vs mpg

    Terlihat bahwa banyak data ketika nilai acceleration [10,20] terlalu jauh dengan garis regresi linear. Hal ini ada dugaan nilai error yang cukup besar. Kemudian di sekitaran nilai acceleration yang berdekatan, nilai mpg-nya cenderung berlainan jauh (dengan nilai x yang hampir sama, nilai y ada yang cukup tinggi (jauh di atas garis regresi) ada yang cukup rendah (jauh di bawah garis regresi)). Berdasarkan visualisasi di atas menunjukkan bahwa hubungan accelartion vs mpg tidak linear.

  • Kasus 3: displacement vs mpg

    Serupa dengan kasus 1, bila hubungan tersebut mengikuti regresi linear, untuk nilai displacment kendaraan (perpindahan posisi) yang besar akan diperoleh nilai mpg (jarak yang dapat ditempuh per gallon bahan bakar) negatif. Hal ini tidak mungkin terjadi atau dapat dikatakan model tidak sesuai konteks

B. Menentukan Model Setiap Pasangan Peubah

1. Hubungan weight vs mpg

set.seed(12345)

# Regresi Polynomial
cv_reg.poly_A = function(derajat, dataset){
  hasil = NULL
  for (i in derajat) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ poly(weight,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "polynomial", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Tangga
cv_fg.tangga_A = function(breaks, dataset){
  hasil = NULL
  for (i in breaks) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ cut(weight,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model ="Fungsi tangga", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Natural Cubic Spline
cv_nc.spline_A = function(df, dataset){
  hasil = NULL
  for (i in df) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ ns(weight,df = ",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "Natural cubic spline", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Komparasi antar model
run_komparasi_A = function(dataset, derajat=2:20, breaks=2:20, df=2:50){
  hasil_poly = cv_reg.poly_A(derajat, dataset)
  hasil_fg.tangga = cv_fg.tangga_A(breaks, dataset)
  hasil_ncspline = cv_nc.spline_A(df,dataset)
  
  all_model = rbind(hasil_poly, hasil_fg.tangga, hasil_ncspline)
  all_model
}


#Fungsi plot
plotmodel_A <- function(model, dataset){
  ggplot(dataset, aes(x=weight, y=mpg)) +
    geom_point(alpha=0.55, color="black") +
    stat_smooth(formula = as.formula(model),
                method= "lm",
                lty=1, lwd=1.1, col= "darkorange2", se=T) +
    labs(title=paste("Final model:", model)) +
    theme_light(base_size = 10)
}

Perbandingan antar model

Pada kasus ini parameter model: derajat pada regresi polynomial, interval fungsi tangga, dan derajat kebebasan modelnaural cubic spline dibatasi dari 2 sampai 10.

comp.model_A <- run_komparasi_A(data3,2:10, 2:10, 2:10)

ggplot(comp.model_A, aes(x=hyper, y=RMSE, group=model)) +
  geom_line(aes(color = model), size=1) +
  labs(title = "Perbandingan RMSE antar model") + 
  theme_light(base_size = 10)

Model terbaik yang menggambarkan hubungan weight terhadap mpg berdasarakan nilai RMSE yang minimum

comp.model_A %>% slice_min(RMSE)
##        model hyper intercept     RMSE  Rsquared      MAE    RMSESD RsquaredSD
## 1 polynomial     2      TRUE 4.114084 0.7231533 3.062698 0.7105695 0.07934069
##       MAESD
## 1 0.4482645

Plot model terbaik

plotmodel_A(model="y ~ ns(x, 2)", dataset=data3)

NB: Pada kasus ini, model terbaik dipilih berdasarkan nilai RMSE

2. Hubungan acceleration vs mpg

set.seed(12345)

# Regresi Polynomial
cv_reg.poly_B = function(derajat, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in derajat) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ poly(acceleration,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "polynomial", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Tangga
cv_fg.tangga_B = function(breaks, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in breaks) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ cut(acceleration,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model ="Fungsi tangga", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Natural Cubic Spline
cv_nc.spline_B = function(df, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in df) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ ns(acceleration,df = ",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "Natural cubic spline", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Komparasi antar model
run_komparasi_B = function(dataset, derajat=2:20, breaks=2:20, df=2:50){
  hasil_poly = cv_reg.poly_B(derajat, dataset)
  hasil_fg.tangga = cv_fg.tangga_B(breaks, dataset)
  hasil_ncspline = cv_nc.spline_B(df,dataset)
  
  all_model = rbind(hasil_poly, hasil_fg.tangga, hasil_ncspline)
  all_model
}


#Fungsi plot
plotmodel_B <- function(model, dataset){
  ggplot(dataset, aes(x=acceleration, y=mpg)) +
    geom_point(alpha=0.55, color="black") +
    stat_smooth(formula = as.formula(model),
                method= "lm",
                lty=1, lwd=1.1, col= "darkorange2", se=T) +
    labs(title=paste("Final model:", model)) +
    theme_light(base_size = 10)
}

Perbandingan antar model

Pada kasus ini parameter model: derajat pada regresi polynomial, interval fungsi tangga, dan derajat kebebasan modelnaural cubic spline dibatasi dari 2 sampai 10.

comp.model_B <- run_komparasi_B(data3,2:10, 2:10, 2:10)

ggplot(comp.model_B, aes(x=hyper, y=RMSE, group=model)) +
  geom_line(aes(color = model), size=1) +
  labs(title = "Perbandingan RMSE antar model") + 
  theme_light(base_size = 10)

Model terbaik yang menggambarkan hubungan weight terhadap mpg berdasarakan nilai RMSE yang minimum

comp.model_B %>% slice_min(RMSE)
##                  model hyper intercept     RMSE  Rsquared      MAE    RMSESD
## 1 Natural cubic spline     7      TRUE 6.910232 0.2321443 5.622919 0.8617502
##   RsquaredSD     MAESD
## 1  0.1264537 0.6085781

Plot model terbaik

plotmodel_B(model="y ~ ns(x,7)", dataset=data3)

NB: Pada kasus ini, model terbaik dipilih berdasarkan nilai RMSE

3. Hubungan displacement vs mpg

set.seed(12345)

# Regresi Polynomial
cv_reg.poly_C = function(derajat, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in derajat) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ poly(displacement,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "polynomial", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Tangga
cv_fg.tangga_C = function(breaks, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in breaks) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ cut(displacement,",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model ="Fungsi tangga", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Natural Cubic Spline
cv_nc.spline_C = function(df, dataset){
  # set.seed(12345)
  hasil = NULL
  for (i in df) {
    train.control <- trainControl(method = "cv", number = 10)
    mdl = paste0("mpg ~ ns(displacement,df = ",i,")")
    model <- train(as.formula(mdl), data = dataset, method = "lm",
                   trControl = train.control)
    tampung = cbind(model = "Natural cubic spline", hyper = i , model$results )
    hasil = rbind(hasil, tampung)
    
    model$finalModel
  }
  hasil
}


# Fungsi Komparasi antar model
run_komparasi_C = function(dataset, derajat=2:20, breaks=2:20, df=2:50){
  hasil_poly = cv_reg.poly_B(derajat, dataset)
  hasil_fg.tangga = cv_fg.tangga_B(breaks, dataset)
  hasil_ncspline = cv_nc.spline_B(df,dataset)
  
  all_model = rbind(hasil_poly, hasil_fg.tangga, hasil_ncspline)
  all_model
}


#Fungsi plot
plotmodel_C <- function(model, dataset){
  ggplot(dataset, aes(x=displacement, y=mpg)) +
    geom_point(alpha=0.55, color="black") +
    stat_smooth(formula = as.formula(model),
                method= "lm",
                lty=1, lwd=1.1, col= "darkorange2", se=T) +
    labs(title=paste("Final model:", model)) +
    theme_light(base_size = 10)
}

Perbandingan antar model

Pada kasus ini parameter model: derajat pada regresi polynomial, interval fungsi tangga, dan derajat kebebasan modelnaural cubic spline dibatasi dari 2 sampai 10.

comp.model_C <- run_komparasi_C(data3,2:10, 2:10, 2:10)

ggplot(comp.model_C, aes(x=hyper, y=RMSE, group=model)) +
  geom_line(aes(color = model), size=1) +
  labs(title = "Perbandingan RMSE antar model") + 
  theme_light(base_size = 10)

Model terbaik yang menggambarkan hubungan weight terhadap mpg berdasarakan nilai RMSE yang minimum

comp.model_C %>% slice_min(RMSE)
##                  model hyper intercept     RMSE  Rsquared      MAE    RMSESD
## 1 Natural cubic spline     7      TRUE 6.910232 0.2321443 5.622919 0.8617502
##   RsquaredSD     MAESD
## 1  0.1264537 0.6085781

Plot model terbaik

plotmodel_C(model="y ~ ns(x,7)", dataset=data3)

NB: Pada kasus ini, model terbaik dipilih berdasarkan nilai RMSE

Referensi: