Moving Beyond Linearity

Pendahuluan

Tujuan

• Melakukan cross validation untuk menghasilkan permodelan mpg vs horsepower yang optimal pada data ISLR :: Auto menggunakan metode

  • Polynomial regression

  • Piecewise constant

  • Natural cubic spline

• Melakukan cross validation untuk masing-masing subset data berdasarkan asal negara produsennya (Amerika, Eropa dan Jepang)

Output

Mendapatkan model non-linier yang optimal berdasarkan masing-masing data negara produsennya

Metodologi

Tahapan

Tuning Parameter

• Cross Validation : Cross validation dengan meminimumkan RMSE (root mean squared error)

• Polynomial regression : Derajat/degree polinomial awal: c(2:10)

• Piecewise constant Jumlah bins/breaks awal: c(2:6)

• Natural Cubic Spline : Jumlah dan lokasi knots awal ditentukan berdasarkan quantile

  Jumlah Knots	Lokasi knots (quantile)
  1	.5,
  2	.33 * c(1:2)
  3	.25 * c(1:3)
  3*	.25 * c(.10 , .5 , .90)
  4	.20 * c(1:4)
  4*	c(.05, .35, .65 , .95)
  5	.167 * c(1:5)
  5*	c(.05, .275, .5, .725, .95)
  6	.143 * c(1:6)
  6*	c(.05, .23, .41, .59, .77, .95)
  7	.125 * c(1:7)
  7*	c(.025, .1833, .3417, .5, .6583, .8167, .975)

Sumber : Frank E. Harrell, Jr. (2015). Regression Modeling Strategies

Tools

• R

• splines dan stats untuk pemodelan

• caret untuk cross validation

• ggplot2 dan cowplot untuk grafik/visualisasi

• Tidyverse untuk data wrangling

Hasil Pembahasan

Fungsi - Fungsi

Membuat Formula Model

library(splines)
library(caret)

## Loading required package: ggplot2

## Loading required package: lattice

library(ggplot2)
library(cowplot)

# Fungsi untuk membuat formula model berdasarkan jenis model dan parameternya
createFormula <- function(y, x, model, param) {
  if (model == "polynomial") {
    f <- sprintf("%s ~ poly(%s, degree = %s)", y, x, param)
  } else if (model == "piecewise") {
    f <- sprintf("%s ~ cut(%s, breaks = %s)", y, x, param)
  } else if (model == "ncubicspline") {
    f <- sprintf("%s ~ ns(%s, knots = c(%s))", y, x, paste(param, collapse = ","))
  }
  return(f)
}

K-Fold Cross Validation

# Fungsi untuk melatih model dengan k-fold cross-validation
trainModel <- function(formula, df) {
  trCtl <- trainControl(method='cv', number = 5)
  models <- train(as.formula(formula), data = df, trControl = trCtl, method = 'glm')
  return(models)
}

Fungsi Performa Model

getModelPerformance <- function(models) {
  RMSE <- models$results$RMSE
  fModel <- models$finalModel
  return(list(RMSE = RMSE, fModel = fModel))
}

Modeling dan Cross Validation

beyondLinear <- function(df, y, x, model, param){
  
  cvResult <- mapply(function(param, model){
    # Buat formula berdasarkan model dan parameter
    formula <- createFormula(y, x, model, param)
    
    # Latih model dengan k-fold cross-validation
    models <- trainModel(formula, df)
    
    # Dapatkan RMSE dan model terbaik
    performance <- getModelPerformance(models)
    
    return(list(models, performance$RMSE, performance$fModel))
    
  }, model = model, param = param, USE.NAMES = FALSE, SIMPLIFY = FALSE )
  if(!is.null(names(param))){
    paramIdx <- names(param)
  } else {
    paramIdx <- param
  }
  
  cvScore <- data.frame(paramIdx = as.character(unlist(paramIdx)),
                        score = unlist(sapply(cvResult, "[", 2, USE.NAMES = F)))
  cvScore$param <- param
  
  bestPos <- which.min(cvScore$score)[1]
  bestParam <- param[[bestPos]]
  bestScore <- cvResult[[bestPos]][[2]]
  bestModel <- cvResult[[bestPos]][[1]]
  
  predValue <- data.frame(df[c(y,x)], 
                          pred = predict(bestModel, df, interval = "prediction"))
  
  return(list(method = model,
              cvScore = cvScore, 
              whichBest = bestPos, 
              bestParamIdx = paramIdx[bestPos],
              bestParam = bestParam, 
              bestScore = bestScore, 
              bestModel = bestModel,
              predValue = predValue))
}

Grafik & Plot

library(ggplot2)
library(cowplot)

cvPlot <- function(df, title, xlab){
  
  #' cvPlot()
  #' 
  #' Create cross-validation plot for Polynomial, Piecewise, Natural Cubic Spline
  #' based on given tuning parameters
  
  scaleFUN <- function(x) as.character(round(x,2))
  df$paramIdx <- factor(df$paramIdx, levels = df$paramIdx)
  
  ggplot(df, aes(x=paramIdx, y=score)) + 
    geom_bar(stat = "identity", width=0.1, fill = "black") +
    coord_cartesian(ylim = c(0.99*min(df$score), NA)) +
    geom_hline(yintercept = min(df$score), linetype="dashed", color = "red") +
    ggtitle(title) +
    xlab(xlab) + ylab("CV Root Mean Squared Error") +
    theme_classic()
}

predPlot <- function(df, y, laby, x, labx, title){
  
  #' predPlot()
  #' 
  #' Create plot of best model's fitted value vs response
  #' Along with its observed values
  
  ggplot(df) +
    geom_point(aes_string(x = x, y = y, fill = as.factor(1))) +
    geom_line(aes_string(x = x, y = "pred", color = as.factor(1)), size = 1) +
    scale_fill_manual(labels = c("Observed"), values = c("black")) +
    scale_color_manual(labels = c("Predicted"), values = c("blue"))  +
    theme_classic() +
    ggtitle(title) +
    xlab(labx) + ylab(laby) +
    theme(legend.title = element_blank(), 
          legend.position = "bottom", 
          legend.box = "horizontal")
}

predPlotAll <- function(df, y, x, laby, labx) {
  
  #' predPlotAll()
  #' 
  #' Create plot of all best model's fitted value vs response 
  #' In one frame
  
  ggplot(df) +
    geom_point(aes_string(x = x, y = y, fill = as.factor(1))) +
    geom_line(aes_string(x = x, y = "polynomial", color = as.factor(2)), size = 1) +
    geom_line(aes_string(x = x, y = "piecewise", color = as.factor(3)), size = 1) +
    geom_line(aes_string(x = x, y = "ncubicspline", color = as.factor(4)), size = 1) +
    scale_fill_manual(labels = c("Observed"), values = c("black")) +
    scale_color_manual(labels = c("Polynomial", "Piecewise", "Natural Cubic Splines"), 
                       values = c("blue", "red", "green"))  +
    theme_classic() +
    ggtitle("Comparison of Best Predictions") +
    xlab(labx) + ylab(laby) +
    theme(legend.title = element_blank(), 
          legend.position = "bottom", 
          legend.box = "horizontal")
}


plotBestModel <- function(df, bestModel, y, x, laby, labx, method, bestParam, titleSuffix) {
  
  title <- paste(titleSuffix, 
                 sprintf(" (%s: %s)", method, bestParam))
  
  p <- ggplot(df) +
    geom_point(aes_string(x = x, y = y), color = "black", size = 1) +
    geom_line(aes_string(x = x, y = "pred"), data = bestModel$predValue, color = "blue", size = 1) +
    scale_fill_manual(labels = c("Observed"), values = c("black")) +
    scale_color_manual(labels = c("Predicted"), values = c("blue")) +
    theme_classic() +
    ggtitle(title) +
    xlab(labx) + 
    ylab(laby) +
    theme(legend.title = element_blank(), 
          legend.position = "bottom", 
          legend.box = "horizontal")
  
  return(p)  # Mengembalikan objek plot
}

Fungsi Gabungan

generateAllREg <- function(df, y, x, laby, labx, degree, cut, knot, best=NULL){
  
  methods <- c("polynomial", "piecewise", "ncubicspline")
  method_alias <- c("Polynomial Regression", "Piecewise Constant", "Natural Cubic Splines")
  params <- list(Degree = degree, `N Bins` = cut, `N Knots` = knot)
  
  outModel <- list()
  pltCV <- list()
  pltPred <- list()
  
  for( i in seq(length(methods)) ){
    
    model <- beyondLinear(df = df, y = y, x = x, model = methods[i], param = params[[i]])
    
    outModel[[i]] <- model
    
    
    p <- cvPlot(model$cvScore,
                title = method_alias[i],
                xlab =  names(params)[i])
    
    pltCV[[i]] <- p
    
    p <- predPlot(df = model$predValue, 
                  y = y, laby = laby,
                  x = x, labx = labx,
                  title = sprintf("%s with %s=%s (RMSE=%s)", 
                                  method_alias[i], names(params)[i], 
                                  model$bestParamIdx, round(model$bestScore,2)))
    
    if(methods[i] == "ncubicspline"){
      p <- p + geom_vline(xintercept = model$bestParam, linetype="dashed", color = "gray")
    }
    
    pltPred[[i]] <- p
  }
  
  bestSummary <- data.frame(t(sapply(outModel, "[", c("method", "bestParamIdx", "bestScore"))))
  names(bestSummary) <- c("method", "param", "RMSE")
  
  predAllBest <- data.frame(outModel[[1]]$predValue, 
                            outModel[[2]]$predValue$pred, 
                            outModel[[3]]$predValue$pred)
  names(predAllBest) <- c(y, x, methods)
  
  pltPredAll <- predPlotAll(df = predAllBest, 
                          y = y, x = x, 
                          laby = laby, labx = labx)
    
    if(is.null(best)){
      bIdx <- which.min(bestSummary$RMSE)
    } else {
      bIdx <- best
    }
  
  bestModel <- bestSummary[bIdx, ]
  
  predicted <- outModel[[bIdx]]$predValue
  
  return(list(bestSummary = bestSummary, 
              bestModel = bestModel,
              bestModelDetail = outModel[[bIdx]]$bestModel,
              bestParamDetail = outModel[[bIdx]]$bestParam,
              model = outModel, 
              pltCV = pltCV, 
              pltPred = pltPred,
              pltPredAll = pltPredAll,
              pltCVGrid = plot_grid(pltCV[[1]], pltCV[[2]], pltCV[[3]]),
              pltPredGrid = plot_grid(pltPred[[1]], pltPred[[2]], pltPred[[3]], pltPredAll),
              predBest3 = predAllBest,
              predBestFinal = predicted))
  
}

Knots Natural Cubic Spline

pctKnots <- function(x){
  
  #' pctKnots()
  #' 
  #' Calculate knots based on quartiles
  
  pct <- list(
    .5,
    .33*c(1:2),
    .25*c(1:3),
    c(.10,.5,.90),
    .20*c(1:4),
    c(.05,.35,.65,.95),
    .167*c(1:5),
    c(.05,.275,.5,.725,.95),
    .143*c(1:6),
    c(.05,.23,.41,.59,.77,.95),
    .125*c(1:7),
    c(.025,.1833,.3417,.5,.6583,.8167,.975))
  
  knt <- sapply(pct, function(p, i){round(quantile(i, p),2)}, i = x)
  names(knt) <- c(1,2,3,"3*",4,"4*",5,"5*",6,"6*",7,"7*")
  return(knt)
  
}

Dataset

Data Auto

library(ISLR)
library(ggplot2)
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ purrr::lift()      masks caret::lift()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(splines)
library(rsample)
library(cowplot)
library(caret)
library(gridExtra)

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

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

Auto2 <- Auto %>% select(mpg, horsepower, origin)

summary(Auto)

##       mpg          cylinders      displacement     horsepower        weight    
##  Min.   : 9.00   Min.   :3.000   Min.   : 68.0   Min.   : 46.0   Min.   :1613  
##  1st Qu.:17.00   1st Qu.:4.000   1st Qu.:105.0   1st Qu.: 75.0   1st Qu.:2225  
##  Median :22.75   Median :4.000   Median :151.0   Median : 93.5   Median :2804  
##  Mean   :23.45   Mean   :5.472   Mean   :194.4   Mean   :104.5   Mean   :2978  
##  3rd Qu.:29.00   3rd Qu.:8.000   3rd Qu.:275.8   3rd Qu.:126.0   3rd Qu.:3615  
##  Max.   :46.60   Max.   :8.000   Max.   :455.0   Max.   :230.0   Max.   :5140  
##                                                                                
##   acceleration        year           origin                      name    
##  Min.   : 8.00   Min.   :70.00   Min.   :1.000   amc matador       :  5  
##  1st Qu.:13.78   1st Qu.:73.00   1st Qu.:1.000   ford pinto        :  5  
##  Median :15.50   Median :76.00   Median :1.000   toyota corolla    :  5  
##  Mean   :15.54   Mean   :75.98   Mean   :1.577   amc gremlin       :  4  
##  3rd Qu.:17.02   3rd Qu.:79.00   3rd Qu.:2.000   amc hornet        :  4  
##  Max.   :24.80   Max.   :82.00   Max.   :3.000   chevrolet chevette:  4  
##                                                  (Other)           :365

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)

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() 
}
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() 
## }

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 Berdasarkan Negara

Asal Negara Produsen Amerika

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

##   mpg horsepower origin
## 1  18        130      1
## 2  15        165      1
## 3  18        150      1
## 4  16        150      1
## 5  17        140      1
## 6  15        198      1

Visualisasi Hubungan mpg terhadap horsepower pada data Amerika

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)
head(AutoEropa)

##    mpg horsepower origin
## 20  26         46      2
## 21  25         87      2
## 22  24         90      2
## 23  25         95      2
## 24  26        113      2
## 51  28         90      2

Visualisasi Hubungan mpg terhadap horsepower terhadap data eropa

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")

Asal Negara Produsen Jepang

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

##    mpg horsepower origin
## 15  24         95      3
## 19  27         88      3
## 30  27         88      3
## 32  25         95      3
## 54  31         65      3
## 55  35         69      3

Visualisasi Hubungan mpg terhadap horsepower terhadap data jepang

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")

Permodelan data Auto

AutoKnots <- pctKnots(Auto$horsepower)

set.seed(1019)
mpg.Auto <- generateAllREg(df = Auto2,
                         y = "mpg", x = "horsepower", 
                         laby = "Miles per Gallon", labx = "Horsepower", 
                         degree = c(2:10), cut = c(2:6), knot = AutoKnots)

## 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.

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Cross Validation

mpg.Auto$pltCVGrid

Tuning parameter terbaik

mpg.Auto$bestSummary

##         method param     RMSE
## 1   polynomial     8 4.324576
## 2    piecewise     6 4.667371
## 3 ncubicspline     5 4.309537

mpg.Auto$bestModel

##         method param     RMSE
## 3 ncubicspline     5 4.309537

mpg.Auto$bestParamDetail

##  16.7%  33.4%  50.1%  66.8%  83.5% 
##  70.00  84.00  93.89 110.00 150.00

mpg.Auto$pltPredGrid

Plot Model Terbaik

bestModelAuto <- mpg.Auto$bestModel
bestParamIdxAuto <- bestModelAuto$param
bestMethodAuto <- bestModelAuto$method

BestPlotAuto <- plotBestModel(df = mpg.Auto$predBestFinal, 
              bestModel = bestModelAuto, 
              y = "mpg", 
              x = "horsepower", 
              laby = "Miles per Gallon", 
              labx = "Horsepower", 
              method = bestMethodAuto, 
              bestParam = bestParamIdxAuto, 
              titleSuffix = "Prediksi Model Terbaik Auto \n")
BestPlotAuto

BestPlotGridAuto<- mpg.Auto$pltPred[3]
BestPlotGridAuto

## [[1]]

Permodelan Data Amerika

AutoKnotsAmerika <- pctKnots(AutoAmerika$horsepower)

set.seed(1019)
mpg.Amerika <- generateAllREg(df = AutoAmerika,
                         y = "mpg", x = "horsepower", 
                         laby = "Miles per Gallon", labx = "Horsepower", 
                         degree = c(2:10), cut = c(2:6), knot = AutoKnotsAmerika)

Cross Validation

mpg.Amerika$pltCVGrid

Tuning Parameter Terbaik

mpg.Amerika$bestSummary

##         method param     RMSE
## 1   polynomial     7 3.786994
## 2    piecewise     5 4.088909
## 3 ncubicspline     1 3.732735

mpg.Amerika$bestModel

##         method param     RMSE
## 3 ncubicspline     1 3.732735

mpg.Amerika$bestParamDetail

## 50% 
## 105

mpg.Amerika$pltPred[3]

## [[1]]

Plot Model Terbaik

bestModelAmerika <- mpg.Amerika$bestModel
bestParamIdxAmerika <- bestModelAmerika$param
bestMethodAmerika <- bestModelAmerika$method

BestPlotAmerika <- plotBestModel(df = mpg.Amerika$predBestFinal, 
              bestModel = bestModelAmerika, 
              y = "mpg", 
              x = "horsepower", 
              laby = "Miles per Gallon", 
              labx = "Horsepower", 
              method = bestMethodAmerika, 
              bestParam = bestParamIdxAmerika, 
              titleSuffix = "Prediksi Model Terbaik (Amerika) \n")
BestPlotAmerika

BestPlotGridAmerika<- mpg.Amerika$pltPred[3]
BestPlotGridAmerika

## [[1]]

Permodelan Data Eropa

AutoKnotsEropa <- pctKnots(AutoEropa$horsepower)

set.seed(1019)
mpg.Eropa <- generateAllREg(df = AutoEropa,
                         y = "mpg", x = "horsepower", 
                         laby = "Miles per Gallon", labx = "Horsepower", 
                         degree = c(2:10), cut = c(2:6), knot = AutoKnotsEropa)

Cross Validation

mpg.Eropa$pltCVGrid

Tuning Parameter Terbaik

mpg.Eropa$bestSummary

##         method param     RMSE
## 1   polynomial     2 4.906158
## 2    piecewise     5 5.084536
## 3 ncubicspline    4* 4.832777

mpg.Eropa$bestModel

##         method param     RMSE
## 3 ncubicspline    4* 4.832777

mpg.Eropa$bestParamDetail

##     5%    35%    65%    95% 
##  48.00  71.00  86.55 115.00

mpg.Eropa$pltPredGrid

Plot Model Terbaik

bestModelEropa <- mpg.Eropa$bestModel
bestParamIdxEropa <- bestModelEropa$param
bestMethodEropa <- bestModelEropa$method

BestPlotEropa <- plotBestModel(df = mpg.Eropa$predBestFinal, 
              bestModel = bestModelEropa, 
              y = "mpg", 
              x = "horsepower", 
              laby = "Miles per Gallon", 
              labx = "Horsepower", 
              method = bestMethodEropa, 
              bestParam = bestParamIdxEropa, 
              titleSuffix = "Prediksi Model Terbaik (Eropa) \n")
BestPlotEropa

BestPlotGridEropa<- mpg.Eropa$pltPred[3]
BestPlotGridEropa

## [[1]]

Permodelan Data Jepang

AutoKnotsJepang <- pctKnots(AutoJepang$horsepower)

set.seed(1019)
mpg.Jepang <- generateAllREg(df = AutoJepang,
                         y = "mpg", x = "horsepower", 
                         laby = "Miles per Gallon", labx = "Horsepower", 
                         degree = c(2:10), cut = c(2:6), knot = AutoKnotsJepang)

Cross Validation

mpg.Jepang$pltCVGrid

mpg.Jepang$bestSummary

##         method param     RMSE
## 1   polynomial     5 4.069017
## 2    piecewise     6 4.252005
## 3 ncubicspline     2 4.002123

mpg.Jepang$bestModel

##         method param     RMSE
## 3 ncubicspline     2 4.002123

mpg.Jepang$bestParamDetail

##   33%   66% 
## 67.74 90.00

mpg.Jepang$pltPredGrid

Plot Model Terbaik

bestModelJepang <- mpg.Jepang$bestModel
bestParamIdxJepang <- bestModelJepang$param
bestMethodJepang <- bestModelJepang$method

BestPlotJepang <- plotBestModel(df = mpg.Jepang$predBestFinal, 
              bestModel = bestModelJepang, 
              y = "mpg", 
              x = "horsepower", 
              laby = "Miles per Gallon", 
              labx = "Horsepower", 
              method = bestMethodJepang, 
              bestParam = bestParamIdxJepang, 
              titleSuffix = "Prediksi Model Terbaik (Jepang) \n")
BestPlotJepang

BestPlotGridJepang<- mpg.Jepang$pltPred[3]
BestPlotGridJepang

## [[1]]

Model Terbaik

grid.arrange(BestPlotAuto,BestPlotAmerika, BestPlotEropa, BestPlotJepang,  ncol = 2)

Interpretasi

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

Auto (Seluruh Data):

Model: Natural cubic spline dengan Knots 5

• Untuk seluruh dataset (Auto), tren umum menunjukkan bahwa saat Horsepower meningkat, mpg cenderung menurun. Pada Horsepower rendah, ada variasi besar dalam mpg, tetapi pada Horsepower yang lebih tinggi, hubungan menjadi lebih stabil dengan mpg cenderung mencapai nilai minimum dan tetap konstan.

• Hal ini masuk akal karena mobil dengan tenaga kuda lebih tinggi umumnya membutuhkan lebih banyak bahan bakar, sehingga menghasilkan mpg yang lebih rendah. Stabilitas pada tenaga kuda yang tinggi mungkin disebabkan oleh batas efisiensi mesin.

Amerika:

Model: Natural cubic spline dengan Knot 1

Interpretasi

• Data Amerika menunjukkan penurunan mpg yang cepat saat Horsepower meningkat, terutama pada kisaran tenaga rendah hingga menengah.

• Spline dengan Knot 1 dmenghasilkan model yang mendekati linear, mengindikasikan bahwa mobil buatan Amerika cenderung memiliki hubungan yang lebih sederhana antara Horsepower dan mpg. Pada kisaran Horsepower yang lebih tinggi, hubungan ini menjadi lebih stabil dengan sedikit perubahan pada mpg.

Eropa:

Model: Natural cubic spline Knot 4*

Interpretasi

• Untuk mobil Eropa, spline dengan Knot 4* menunjukkan lebih banyak variasi di antara hubungan Horsepower dan mpg dibandingkan dengan data Amerika. Pada kisaran tenaga rendah hingga menengah, mpg menurun dengan beberapa fluktuasi yang lebih halus.

• ini mungkin mencerminkan keanekaragaman desain mobil Eropa yang umumnya lebih efisien dalam hal tenaga dan konsumsi bahan bakar. Pada kisaran tenaga rendah, mobil Eropa menunjukkan kinerja mpg yang lebih baik dibandingkan mobil Amerika.

Jepang:

Model: Natural cubic spline dengan Knot 2.

Interpretasi:

• Data mobil Jepang menunjukkan tren non-linear yang menarik. Pada kisaran Horsepower rendah hingga menengah, ada penurunan mpg yang diikuti oleh kenaikan kembali pada kisaran Horsepower menengah.

• Ini mungkin menunjukkan bahwa mobil Jepang memiliki desain yang lebih bervariasi dengan fokus pada efisiensi bahan bakar, terutama pada kisaran Horsepower menengah.

Kesimpulan

1. 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.

2. 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.

3. 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.