Nonlinear Regression

Soal

Libraries

library(ISLR)
library(tidyverse)
library(splines)    # Natural Cubic Spline
library(rsample)    # validasi silang k fold (vfold_cv)
library(cowplot)    # plot_grid
library(DT)         # datatable

Dataset

Auto Dataset

Tugas memodelkan hubungan antara variabel horsepower dan mpg, weight dan mpg, serta displacement dan mpg dari dataset Auto (data mobil) yang disediakan oleh package ISLR. Datanya terlihat seperti berikut.

data("Auto",package="ISLR")
attach(Auto)
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

Pemodelan horsepower dan mpg

Visualisasi Data

ggplot(Auto,aes(x=horsepower, y=mpg)) +
  geom_point(alpha=2.55, color="coral") + 
  theme_bw()

Dari 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 pula bahwa hubungan antara kekuatan mesin kendaraan (pada sumbu X) dan jarak yang dapat ditempuh oleh kendaraan per galon bahan bakarnya (pada sumbu Y) bersifat tidak linear.

Jika dipaksakan dengan regresi linear sederhana, maka akan terlihat seperti berikut.

ggplot(Auto,aes(x=horsepower, y=mpg)) +
  geom_point(alpha=2.05, color="coral") +
  stat_smooth(method = "lm", 
              formula = y~x,lty = 1,
              col = "blue",se = F)+
  theme_bw()

Tampak bahwa hubungan linear tidak masuk akal karena seolah jika mengikuti garis linearnya maka kendaraan dengan kekuatan mesin yang lebih besar lagi akan menempuh jarak nol atau bahkan negatif per galon bahan bakarnya. Ini jelas tidak mungkin.

Oleh karena regresi linear tidak cocok, maka untuk memodelkan masalah ini, dapat digunakan pendekatan Regresi Polinomial, Piecewise Constant (Fungsi Tangga), atau Natural Cubic Splines.

Regresi Polinomial

Cross Validation: untuk memilih model terbaik

polinomial <- function(dataset){
  bestmodel = c()
  derajat = 1:10
  for (i in derajat){
    set.seed(101)
    cross_val <- vfold_cv(dataset,v=10,strata = "mpg") #dari rsample
    metric_poly <- map_dfr(cross_val$splits, #dari tydiverse
                           function(x){  
                             mod <- lm(mpg ~ poly(horsepower,i,raw = T), #memodelkan data pada data training
                                       data=dataset[x$in_id,]) 
                             pred <- predict(mod, # nilai prediksi untuk data testing
                                             newdata=dataset[-x$in_id,]) 
                             truth <- dataset[-x$in_id,]$mpg
                             rmse <- mlr3measures::rmse(truth = truth,response = pred)
                             mae <- mlr3measures::mae(truth = truth,response = pred)
                             metric <- c(rmse,mae) 
                             names(metric) <- c("rmse","mae")
                             return(metric)
                           }
    )
    average = colMeans(metric_poly)
    bestmodel = c(bestmodel,average)
  }
  
  M = matrix(bestmodel,byrow = T,ncol = 2)
  colnames(M) <- c("rmse","mae")
  bestModel = as.data.frame(cbind(derajat,M))
  return(bestModel)
}
bestModel = polinomial(Auto); bestModel

##    derajat     rmse      mae
## 1        1 4.879566 3.848220
## 2        2 4.329478 3.266067
## 3        3 4.335502 3.269710
## 4        4 4.340408 3.276116
## 5        5 4.293865 3.237505
## 6        6 4.280648 3.236147
## 7        7 4.257757 3.205254
## 8        8 4.263997 3.212749
## 9        9 4.279362 3.213649
## 10      10 4.336515 3.260937

Dari tabel di atas, terlihat bahwa regresi polinomial dengan derajat 7 memiliki rmse maupun mae terkecil. Agar lebih mudah membandingkan, dapat dilihat pada plot berikut.

bestModel$derajat <- as.factor(bestModel$derajat)
ggplot(data=bestModel[-1,], aes(x=derajat, y=rmse, group=1)) +
  geom_line(linetype = "dashed")+
  geom_point()+
  ggtitle("K-Fold Cross Validation") +
  xlab("Degree of Polynomial") + 
  ylab("Root Mean Square Error")

Secara empiris regresi polinomial dengan derajat 7 memiliki rmse maupun mae terkecil. Namun, di sini dipilih regresi polinomial dengan derajat 2 karena menurut kebiasaan dalam regresi polinomial hanya digunakan sampai dejarat 4.

Hal ini karena meskipun memiliki rmse atau mae yang kecil, polinomial dengan derajat lebih besar dari 4 cenderung memiliki pola yang tidak teratur, khususnya pada ujung-ujung kurvanya. Untuk lebih jelas, dapat dilihat kurva berikut.

p2 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  ggtitle("Regresi Polinomial Derajat 2")

p7 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,7,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw() +
  ggtitle("Regresi Polinomial Derajat 7")

plot_grid(p2, p7)

Terlihat bahwa polinomial dengan derajat 7 memiliki ujung-ujung kurva yang tidak stabil.

Fungsi Tangga (Piecewise Constant)

Secara empiris

cv.pc <- function(dataset){
  set.seed(101)
  cross_val <- vfold_cv(dataset,v=10,strata = "mpg")
  interval <- 2:18
  best_tangga <- map_dfr(interval, 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)
                               mae <- mlr3measures::mae(truth = data_eval$truth,
                                                        response = data_eval$pred)
                               metric <- c(rmse,mae)
                               names(metric) <- c("rmse","mae")
                               return(metric)
                             }
    )
    colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
  }
  )
  best_tangga <- cbind(interval=interval,best_tangga) # menampilkan hasil all breaks
  basedonrmse <- as.data.frame(best_tangga %>% slice_min(rmse))
  basedonmae <- as.data.frame(best_tangga %>% slice_min(mae))
  return(list(tabel=best_tangga,best=rbind(basedonrmse,basedonmae)))
}
tangga = cv.pc(Auto); tangga

## $tabel
##    interval     rmse      mae
## 1         2 6.101201 4.766114
## 2         3 5.767713 4.520783
## 3         4 4.970451 3.782808
## 4         5 4.689346 3.556735
## 5         6 4.634534 3.524957
## 6         7 4.488716 3.366269
## 7         8 4.416584 3.409604
## 8         9 4.490070 3.467869
## 9        10 4.533882 3.416048
## 10       11 4.496908 3.337706
## 11       12 4.437961 3.309810
## 12       13 4.499759 3.362132
## 13       14 4.331626 3.231907
## 14       15 4.325459 3.291656
## 15       16 4.329302 3.343669
## 16       17 4.355959 3.329588
## 17       18 4.352910 3.234846
## 
## $best
##   interval     rmse      mae
## 1       15 4.325459 3.291656
## 2       14 4.331626 3.231907

Secara empiris berdasarkan rmse banyaknya interval yang terbaik adalah 15 sedangkan berdasarkan mae banyaknya interval terbaik adalah 14.

Secara visual

pc15 <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,15), 
               lty = 1, col = "blue",se = F)+
  theme_bw() +
  ggtitle("Regresi Fungsi Tangga 15 interval")

pc14 <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,14), 
               lty = 1, col = "blue",se = F)+
  theme_bw() +
  ggtitle("Regresi Fungsi Tangga 14 interval")


plot_grid(pc15, pc14)

Berdasarkan gambar, yang 15 interval lebih baik dibandingkan yang 14 interval karene memiliki pola yang lebih teratur, sehingga dipilih yang 15 interval.

Natural Cubic Spline

Menentukan Banyaknya Knot

nspline = function(dataset){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10,strata = "mpg")
  df <- 1:20
  best.spline3 <- map_dfr(df, function(i){
    metric.spline3 <- map_dfr(cross.val$splits,
                              function(x){
                                mod <- lm(mpg ~ ns(horsepower,df=i),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)
                                mae <- mlr3measures::mae(truth = truth,response = pred)
                                metric <- c(rmse,mae)
                                names(metric) <- c("rmse","mae")
                                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(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                          best.spline3 %>% slice_min(mae))) #berdasarkan mae

  df=basis[,"df"]
  knot1 = attr(ns(dataset$horsepower, df=df[1]),"knots")
  knot2 = attr(ns(dataset$horsepower, df=df[2]),"knots")
  length = c(paste0("length knot1 = ",length(knot1)),paste0("length knot2 = ",length(knot2)))
  return(list(basis=basis,knot1=round(knot1,2),knot2=round(knot2,2),length=length))
}

basisAuto = nspline(Auto);basisAuto

## $basis
##   df     rmse      mae
## 1 10 4.249102 3.184156
## 2 12 4.275312 3.182377
## 
## $knot1
##   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 
## 
## $knot2
## 8.333333% 16.66667%       25% 33.33333% 41.66667%       50% 58.33333% 66.66667% 
##     65.00     70.00     75.00     84.00     88.00     93.50    100.00    110.00 
##       75% 83.33333% 91.66667% 
##    126.00    150.00    165.83 
## 
## $length
## [1] "length knot1 = 9"  "length knot2 = 11"

Dari hasil di atas, berdasarkan rmse dipilih 9 knot (df=10), sedangkan berdasarkan mae dipilih 11 knot (df=12).

Cross Validation

k1 = basisAuto$knot1
k2 = basisAuto$knot2
nsplineval = function(dataset,k1,k2){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10,strata = "mpg")
  
  metric1 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(horsepower, knots = k1),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm1 <- colMeans(metric1)# menghitung rata-rata untuk 10 folds
  
  
  metric2 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(horsepower, knots = k2),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm2 <- colMeans(metric2)# menghitung rata-rata untuk 10 folds
  
  
  nilai.cv.ns <- rbind("df1" = meanm1,"df2" = meanm2)
  
  nama.model.ncs <- c(paste0("knot=",length(k1)),paste0("knot=",length(k2)))
  
  model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)
  
  return(model.ncs)
}
model.ncs = nsplineval(Auto,k1,k2); model.ncs

##       Model     rmse      mae
## df1  knot=9 4.247127 3.180381
## df2 knot=11 4.274068 3.179753

Secara empiris berdasarkan rmse banyaknya knot terbaik adalah 9 knot, berdasarkan mae banyaknya knot terbaik adalah 11 knot.

Visualisasi Grafik

nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = k1), col = "blue", se = F) +
  theme_bw()+
  ggtitle("Regresi Natural Cubic Spline 9 Knot")


nc11 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = k2), col = "blue", se = F)+
  theme_bw()+
  ggtitle("Regresi Natural Cubic Spline 11 Knot")

plot_grid(nc9, nc11)

Secara visual, grafik dengan 9 knot terlihat lebih rapi daripada grafik dengan 11 knot, sehingga dipilih grafik dengan 9 knot.

Perbandingan Model

Secara Empiris

a <- as.matrix(bestModel[2,])
b <- as.matrix(tangga$best[1,])
c <- as.matrix(model.ncs[1,])
perbandingan.model <- rbind(a,b,c)
rownames(perbandingan.model) <- c("Polinomial Derajat 2","Fungsi Tangga 15 Interval",
                                  "Natural Cubic Spline 9 Knot")
as.data.frame(perbandingan.model[,-1])

##                                         rmse              mae
## Polinomial Derajat 2                4.329478         3.266067
## Fungsi Tangga 15 Interval   4.32545924917122 3.29165646820982
## Natural Cubic Spline 9 Knot         4.247127         3.180381

Berdasarkan nilai rmse dan mae, model terbaik adalah Natural Cubic Spline dengan 9 Knot. Namun ini tidak dipilih karena dari segi grafik tidak mendukung. Berikut grafiknya.

plot_grid(p2,pc15,nc9)

Secara grafik Regresi Polinomial Derajat 2 terlihat lebih baik dibandingkan Natural Cubic Spline yang mana di ujungnya cenderung tidak stabil meskipun memiliki rmse dan mea yang lebih kecil. Sehingga model terbaik yang dipilih adalah Regresi Polinomial Derajat 2.

Pemodelan weight dan mpg

Visualisasi Data

ggplot(Auto,aes(x=weight, y=mpg)) +
  geom_point(alpha=2.55, color="coral") + 
  theme_bw()

Sebagaimana scatter plot antara horsepower dan mpg, scatter plot antara weight dan mpg juga terlihat tidak linear. Secara umum terlihat bahwa semakin berat suatu mobil, semakin banyak bahan bakar yang dikonsumsi, namun jarak yang ditempuh makin pendek.

Pemodelan hubungan antara weight dan mpg di sini juga menggunakan Regresi Polinomial, Piecewise Constant (Fungsi Tangga), atau Natural Cubic Splines.

Regresi Polinomial

Cross Validation: untuk memilih model terbaik

polinomial <- function(dataset){
  bestmodel = c()
  derajat = 1:10
  for (i in derajat){
    set.seed(101)
    cross_val <- vfold_cv(dataset,v=10) #dari rsample
    metric_poly <- map_dfr(cross_val$splits, #dari tydiverse
                           function(x){  
                             mod <- lm(mpg ~ poly(weight,i,raw = T), #memodelkan data pada data training
                                       data=dataset[x$in_id,]) 
                             pred <- predict(mod, # nilai prediksi untuk data testing
                                             newdata=dataset[-x$in_id,]) 
                             truth <- dataset[-x$in_id,]$mpg
                             rmse <- mlr3measures::rmse(truth = truth,response = pred)
                             mae <- mlr3measures::mae(truth = truth,response = pred)
                             metric <- c(rmse,mae) 
                             names(metric) <- c("rmse","mae")
                             return(metric)
                           }
    )
    average = colMeans(metric_poly)
    bestmodel = c(bestmodel,average)
  }
  
  M = matrix(bestmodel,byrow = T,ncol = 2)
  colnames(M) <- c("rmse","mae")
  bestModel = as.data.frame(cbind(derajat,M))
  brmse <- as.data.frame(bestModel %>% slice_min(rmse))
  bmae <- as.data.frame(bestModel %>% slice_min(mae))
  return(list(tabel=bestModel,best=rbind(brmse,bmae)))
}
bestModel = polinomial(Auto); bestModel

## $tabel
##    derajat     rmse      mae
## 1        1 4.319529 3.300646
## 2        2 4.142592 3.074852
## 3        3 4.146418 3.076157
## 4        4 4.153659 3.081786
## 5        5 4.156560 3.089703
## 6        6 4.159830 3.104884
## 7        7 4.150325 3.109720
## 8        8 4.164884 3.122320
## 9        9 4.286755 3.181968
## 10      10 4.256240 3.154162
## 
## $best
##   derajat     rmse      mae
## 1       2 4.142592 3.074852
## 2       2 4.142592 3.074852

Dari tabel di atas, terlihat bahwa regresi polinomial dengan derajat 2 memiliki rmse maupun mae terkecil. Agar lebih mudah membandingkan, dapat dilihat pada plot berikut.

bestModel$tabel$derajat <- as.factor(bestModel$tabel$derajat)
ggplot(data=bestModel$tabel[-1,], aes(x=derajat, y=rmse, group=1)) +
  geom_line(linetype = "dashed")+
  geom_point()+
  ggtitle("K-Fold Cross Validation") +
  xlab("Degree of Polynomial") + 
  ylab("Root Mean Square Error")

Terlihat bahwa regresi polinomial dengan derajat 2 memiliki rmse terkecil.

Berikut plot regresi polinomial dengan derajat 2.

p2 <- ggplot(Auto,aes(x=weight, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  ggtitle("Regresi Polinomial Derajat 2")

p2

Fungsi Tangga (Piecewise Constant)

Secara empiris

cv.pc <- function(dataset){
  set.seed(101)
  cross_val <- vfold_cv(dataset,v=10)
  interval <- 2:15
  best_tangga <- map_dfr(interval, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
                             function(x){
                               training <- dataset[x$in_id,]
                               training$weight <- cut(training$weight,i)
          
                               mod <- lm(mpg ~ weight,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,]
                               weight_new <- cut(testing$weight,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
                               pred <- predict(mod,
                                               newdata=list(weight=weight_new))
                               truth <- testing$mpg
                               
                               data_eval <- na.omit(data.frame(truth,pred))
                               
                               rmse <- mlr3measures::rmse(truth = data_eval$truth,
                                                          response = data_eval$pred)
                               mae <- mlr3measures::mae(truth = data_eval$truth,
                                                        response = data_eval$pred)
                               metric <- c(rmse,mae)
                               names(metric) <- c("rmse","mae")
                               return(metric)
                             }
    )
    colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
  }
  )
  best_tangga <- cbind(interval=interval,best_tangga) # menampilkan hasil all breaks
  brmse <- as.data.frame(best_tangga %>% slice_min(rmse))
  bmae <- as.data.frame(best_tangga %>% slice_min(mae))
  return(list(tabel=best_tangga,best=rbind(brmse,bmae)))
}
tangga = cv.pc(Auto); tangga

## $tabel
##    interval     rmse      mae
## 1         2 5.629421 4.387157
## 2         3 4.831505 3.642780
## 3         4 4.697535 3.613520
## 4         5 4.380339 3.250564
## 5         6 4.166984 3.098930
## 6         7 4.272382 3.186792
## 7         8 4.429601 3.277997
## 8         9 4.267901 3.138984
## 9        10 4.320372 3.187828
## 10       11 4.288085 3.196973
## 11       12 4.133035 3.034924
## 12       13 4.233133 3.149525
## 13       14 4.173105 3.092575
## 14       15 4.223529 3.110735
## 
## $best
##   interval     rmse      mae
## 1       12 4.133035 3.034924
## 2       12 4.133035 3.034924

Secara empiris berdasarkan rmse dan mae banyaknya interval terbaik adalah 12. Berikut visualisasi dari regresi fungsi tangga 12 interval.

pc12 <- ggplot(Auto,aes(x=weight, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,12), 
               lty = 1, col = "blue",se = F)+
  theme_bw() +
  ggtitle("Regresi Fungsi Tangga 12 interval")

pc12

Natural Cubic Spline

Menentukan Banyaknya Knot

nspline = function(dataset){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10)
  df <- 1:20
  best.spline3 <- map_dfr(df, function(i){
    metric.spline3 <- map_dfr(cross.val$splits,
                              function(x){
                                mod <- lm(mpg ~ ns(weight,df=i),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)
                                mae <- mlr3measures::mae(truth = truth,response = pred)
                                metric <- c(rmse,mae)
                                names(metric) <- c("rmse","mae")
                                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(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                          best.spline3 %>% slice_min(mae))) #berdasarkan mae

  df=basis[,"df"]
  knot1 = attr(ns(dataset$weight, df=df[1]),"knots")
  knot2 = attr(ns(dataset$weight, df=df[2]),"knots")
  length = c(paste0("length knot1 = ",length(knot1)),paste0("length knot2 = ",length(knot2)))
  return(list(basis=basis,knot1=round(knot1,2),knot2=round(knot2,2),length=length))
}

basisAuto = nspline(Auto);basisAuto

## $basis
##   df     rmse      mae
## 1  7 4.126169 3.066635
## 2  7 4.126169 3.066635
## 
## $knot1
## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429% 
##   2074.86   2285.43   2635.00   2986.57   3446.14   4096.29 
## 
## $knot2
## 14.28571% 28.57143% 42.85714% 57.14286% 71.42857% 85.71429% 
##   2074.86   2285.43   2635.00   2986.57   3446.14   4096.29 
## 
## $length
## [1] "length knot1 = 6" "length knot2 = 6"

Dari hasil di atas, berdasarkan rmse dan mae dipilih 6 knot (df=7).

Cross Validation

k1 = basisAuto$knot1
k2 = basisAuto$knot2
nsplineval = function(dataset,k1,k2){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10,strata = "mpg")
  
  metric1 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(weight, knots = k1),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm1 <- colMeans(metric1)# menghitung rata-rata untuk 10 folds
  
  
  metric2 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(weight, knots = k2),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm2 <- colMeans(metric2)# menghitung rata-rata untuk 10 folds
  
  
  nilai.cv.ns <- rbind("df1" = meanm1,"df2" = meanm2)
  
  nama.model.ncs <- c(paste0("knot=",length(k1)),paste0("knot=",length(k2)))
  
  model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)
  
  return(model.ncs)
}
model.ncs = nsplineval(Auto,k1,k2); model.ncs

##      Model     rmse    mae
## df1 knot=6 4.136846 3.0575
## df2 knot=6 4.136846 3.0575

Terlihat bahwa model ini memiliki rmse 4.136846 dan mae 3.0575. Adapun grafiknya adalah sebagai berikut.

nc6 <- ggplot(Auto, aes(x = weight, y = mpg)) + 
  geom_point(alpha = 0.55, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = k1), col = "blue", se = F) +
  theme_bw()+
  ggtitle("Regresi Natural Cubic Spline 6 Knot")

nc6

Perbandingan Model

Secara Empiris

a <- as.matrix(bestModel$tabel[2,])
b <- as.matrix(tangga$best[1,])
c <- as.matrix(model.ncs[1,])
perbandingan.model <- rbind(a,b,c)
rownames(perbandingan.model) <- c("Polinomial Derajat 2","Fungsi Tangga 12 Interval",
                                  "Natural Cubic Spline 6 Knot")
as.data.frame(perbandingan.model[,-1])

##                                         rmse              mae
## Polinomial Derajat 2                4.142592         3.074852
## Fungsi Tangga 12 Interval   4.13303464932709 3.03492445640165
## Natural Cubic Spline 6 Knot         4.136846           3.0575

Berdasarkan nilai rmse dan mae, model terbaik adalah Fungsi Tangga 12 Interval. Berikut grafiknya.

pc12

Pemodelan displacement dan mpg

Visualisasi Data

ggplot(Auto,aes(x=displacement, y=mpg)) +
  geom_point(alpha=2.55, color="coral") + 
  theme_bw()

Sebagaimana scatter plot sebelumnya, scatter plot antara displacement dan mpg juga terlihat tidak linear.

Pemodelan hubungan antara displacement dan mpg di sini juga menggunakan Regresi Polinomial, Piecewise Constant (Fungsi Tangga), atau Natural Cubic Splines.

Regresi Polinomial

Cross Validation: untuk memilih model terbaik

polinomial <- function(dataset){
  bestmodel = c()
  derajat = 1:10
  for (i in derajat){
    set.seed(101)
    cross_val <- vfold_cv(dataset,v=10) #dari rsample
    metric_poly <- map_dfr(cross_val$splits, #dari tydiverse
                           function(x){  
                             mod <- lm(mpg ~ poly(displacement,i,raw = T), #memodelkan data pada data training
                                       data=dataset[x$in_id,]) 
                             pred <- predict(mod, # nilai prediksi untuk data testing
                                             newdata=dataset[-x$in_id,]) 
                             truth <- dataset[-x$in_id,]$mpg
                             rmse <- mlr3measures::rmse(truth = truth,response = pred)
                             mae <- mlr3measures::mae(truth = truth,response = pred)
                             metric <- c(rmse,mae) 
                             names(metric) <- c("rmse","mae")
                             return(metric)
                           }
    )
    average = colMeans(metric_poly)
    bestmodel = c(bestmodel,average)
  }
  
  M = matrix(bestmodel,byrow = T,ncol = 2)
  colnames(M) <- c("rmse","mae")
  bestModel = as.data.frame(cbind(derajat,M))
  brmse <- as.data.frame(bestModel %>% slice_min(rmse))
  bmae <- as.data.frame(bestModel %>% slice_min(mae))
  return(list(tabel=bestModel,best=rbind(brmse,bmae)))
}
bestModel = polinomial(Auto); bestModel

## $tabel
##    derajat     rmse      mae
## 1        1 4.629839 3.533917
## 2        2 4.348269 3.180433
## 3        3 4.354788 3.193282
## 4        4 4.365079 3.199783
## 5        5 4.370460 3.202688
## 6        6 4.349644 3.225788
## 7        7 4.294501 3.225329
## 8        8 4.247949 3.234253
## 9        9 4.193352 3.171016
## 10      10 4.162002 3.149607
## 
## $best
##   derajat     rmse      mae
## 1      10 4.162002 3.149607
## 2      10 4.162002 3.149607

Dari tabel di atas, terlihat bahwa regresi polinomial dengan derajat 10 memiliki rmse maupun mae terkecil. Agar lebih mudah membandingkan, dapat dilihat pada plot berikut.

bestModel$tabel$derajat <- as.factor(bestModel$tabel$derajat)
ggplot(data=bestModel$tabel[-1,], aes(x=derajat, y=rmse, group=1)) +
  geom_line(linetype = "dashed")+
  geom_point()+
  ggtitle("K-Fold Cross Validation") +
  xlab("Degree of Polynomial") + 
  ylab("Root Mean Square Error")

Secara empiris regresi polinomial dengan derajat 10 memiliki rmse maupun mae terkecil. Namun, di sini dipilih regresi polinomial dengan derajat 2 karena menurut kebiasaan dalam regresi polinomial hanya digunakan sampai dejarat 4.

p2 <- ggplot(Auto,aes(x=displacement, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  ggtitle("Regresi Polinomial Derajat 2")

p10 <- ggplot(Auto,aes(x=displacement, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,10,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  ggtitle("Regresi Polinomial Derajat 10")

plot_grid(p2,p10)

Terlihat bahwa polinomial dengan derajat 10 memiliki cenderung overfitting dan ujung-ujung kurva yang tidak stabil.

Fungsi Tangga (Piecewise Constant)

Secara empiris

cv.pc <- function(dataset){
  set.seed(101)
  cross_val <- vfold_cv(dataset,v=10)
  interval <- 2:12
  best_tangga <- map_dfr(interval, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
                             function(x){
                               training <- dataset[x$in_id,]
                               training$displacement <- cut(training$displacement,i)
          
                               mod <- lm(mpg ~ displacement,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,]
                               displacement_new <- cut(testing$displacement,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
                               pred <- predict(mod,
                                               newdata=list(displacement=displacement_new))
                               truth <- testing$mpg
                               
                               data_eval <- na.omit(data.frame(truth,pred))
                               
                               rmse <- mlr3measures::rmse(truth = data_eval$truth,
                                                          response = data_eval$pred)
                               mae <- mlr3measures::mae(truth = data_eval$truth,
                                                        response = data_eval$pred)
                               metric <- c(rmse,mae)
                               names(metric) <- c("rmse","mae")
                               return(metric)
                             }
    )
    colMeans(metric_tangga) # menghitung rata-rata untuk 10 folds
  }
  )
  best_tangga <- cbind(interval=interval,best_tangga) # menampilkan hasil all breaks
  brmse <- as.data.frame(best_tangga %>% slice_min(rmse))
  bmae <- as.data.frame(best_tangga %>% slice_min(mae))
  return(list(tabel=best_tangga,best=rbind(brmse,bmae)))
}
tangga = cv.pc(Auto); tangga

## $tabel
##    interval     rmse      mae
## 1         2 5.961505 4.664498
## 2         3 4.936135 3.694522
## 3         4 4.866941 3.654998
## 4         5 4.710877 3.526767
## 5         6 4.649193 3.437803
## 6         7 4.613718 3.393595
## 7         8 4.399489 3.261074
## 8         9 4.242120 3.118296
## 9        10 4.354184 3.161662
## 10       11 4.394856 3.260741
## 11       12 4.510986 3.322842
## 
## $best
##   interval    rmse      mae
## 1        9 4.24212 3.118296
## 2        9 4.24212 3.118296

Secara empiris berdasarkan rmse dan mae banyaknya interval terbaik adalah 9. Berikut visualisasi dari regresi fungsi tangga 9 interval.

pc9 <- ggplot(Auto,aes(x=displacement, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,9), 
               lty = 1, col = "blue",se = F)+
  theme_bw() +
  ggtitle("Regresi Fungsi Tangga 9 interval")

pc9

Natural Cubic Spline

Menentukan Knot

Berikut saya mencoba 2 knot dan 3 knot

Cross Validation

k1 = c(180,280)
k2 = c(180,280,330)
nsplineval = function(dataset,k1,k2){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10,strata = "mpg")
  
  metric1 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(displacement, knots = k1),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm1 <- colMeans(metric1)# menghitung rata-rata untuk 10 folds
  
  
  metric2 <- map_dfr(cross.val$splits,
                     function(x){
                       mod <- lm(mpg ~ ns(displacement, knots = k2),
                                 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)
                       mae <- mlr3measures::mae(truth = truth,
                                                response = pred)
                       metric <- c(rmse,mae)
                       names(metric) <- c("rmse","mae")
                       return(metric)
                     }
  )
  meanm2 <- colMeans(metric2)# menghitung rata-rata untuk 10 folds
  
  
  nilai.cv.ns <- rbind("df1" = meanm1,"df2" = meanm2)
  
  nama.model.ncs <- c(paste0("knot=",length(k1)),paste0("knot=",length(k2)))
  
  model.ncs<-data.frame(Model = nama.model.ncs, nilai.cv.ns)
  
  return(model.ncs)
}
model.ncs = nsplineval(Auto,k1,k2); model.ncs

##      Model     rmse      mae
## df1 knot=2 4.310693 3.182200
## df2 knot=3 4.311277 3.175623

Terlihat bahwa berdasarkan rmse model terbaik adalah 2 knot dan berdasarkan mae model terbaik adalah 3 knot. Adapun grafiknya adalah sebagai berikut.

nc2 <- ggplot(Auto, aes(x = displacement, y = mpg)) + 
  geom_point(alpha = 0.55, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = k1), col = "blue", se = F) +
  theme_bw()+
  ggtitle("Regresi Natural Cubic Spline 2 Knot")

nc3 <- ggplot(Auto, aes(x = displacement, y = mpg)) + 
  geom_point(alpha = 0.55, color="black") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = k1), col = "blue", se = F) +
  theme_bw()+
  ggtitle("Regresi Natural Cubic Spline 3 Knot")

plot_grid(nc2,nc3)

Di sini, saya memilih Natural Cubic Spline 3 knot

Perbandingan Model

Secara Empiris

a <- as.matrix(bestModel$tabel[2,])
b <- as.matrix(tangga$best[1,])
c <- as.matrix(model.ncs[2,])
perbandingan.model <- rbind(a,b,c)
rownames(perbandingan.model) <- c("Polinomial Derajat 2","Fungsi Tangga 9 Interval",
                                  "Natural Cubic Spline 3 Knot")
as.data.frame(perbandingan.model[,-1])

##                                         rmse              mae
## Polinomial Derajat 2                4.348269         3.180433
## Fungsi Tangga 9 Interval    4.24211993929606 3.11829599136929
## Natural Cubic Spline 3 Knot         4.311277         3.175623

Berdasarkan nilai rmse dan mae, model terbaik adalah Fungsi Tangga 9 Interval. Namun, saya memilih Natural Cubic Spline 3 Knot karena cukup baik merepresentasikan data dari segi grafik. Dapat dilihat gambar berikut.

plot_grid(p2,pc9,nc3)