Tugas Mandiri P8 Sains Data

summary(cars)

##      speed           dist       
##  Min.   : 4.0   Min.   :  2.00  
##  1st Qu.:12.0   1st Qu.: 26.00  
##  Median :15.0   Median : 36.00  
##  Mean   :15.4   Mean   : 42.98  
##  3rd Qu.:19.0   3rd Qu.: 56.00  
##  Max.   :25.0   Max.   :120.00

Install Package yg diperlukan

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
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── 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)
library(ISLR)
library(cowplot)

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:lubridate':
## 
##     stamp

library(caret) # cross-validation

## 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)
library(DT)
library(kableExtra)

## 
## Attaching package: 'kableExtra'

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

library(gam) # generalized additive models

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loaded gam 1.22-5

Eksplorasi Data

Data Auto

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

mpg : miles per gallon

horsepower: Engine horsepower

origin: Origin of car (1. American, 2. European, 3. Japanese)

Auto = Auto %>%  select(mpg, horsepower, origin)
datatable(Auto)

ggplot(Auto,aes(x=horsepower,y=mpg))+geom_point(alpha=0.55,color="red")+theme()+
  labs(y="Mile per Gallon",x="Horse Power",title="Scatter Plot 'mpg' dan 'horsepower'",
       subtitle="World")

a. Data Asal Produsen Amerika

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

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

b. Data Asal Produsen Eropa

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

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

c. Data Asal Produsen Jepang

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

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

1a. Data Auto Regresi Polinomial

Secara Empiris

cross validation untuk menghasilkan pemodelan mpg vs horsepower optimal dengan regresi polinomial

set.seed(232)

cv.poly1 <- function(dataset){
  set.seed(232)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  RegLin<- colMeans(metric_poly1) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(RegLin) 
}

cv.poly2 <- function(dataset){
  set.seed(232)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(232)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(232)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      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      mae
## [1,] 4.865781 3.835360
## [2,] 4.321076 3.258824
## [3,] 4.320929 3.259139
## [4,] 4.328460 3.268511

Jika dilihat dari RMSE dan MAE terlihat bahwa polinomial pangkat 3 memberikan nilai RMSE dan MAE yang terkecil sehingga dipilih polinomial pangkat 3.

Secara Visual

poly1 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,1,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


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

poly3 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly4 <- ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,4,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

plot_grid(poly1, poly2, poly3, poly4, labels=c("RegLin","poly2", "poly3", "poly4"), label_size = 6)

1b. Data Auto Piecewise Constant (Fungsi Tangga)

Secara Empiris

set.seed(252)

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

breaks <- 3:12

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
                           )
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(breaks=breaks,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(rbind(basedonrmse,basedonmae))
}

cv.pc(Auto)

##   breaks     rmse      mae
## 1     12 4.438717 3.310112
## 2     12 4.438717 3.310112

Secara empiris berdasarkan RMSE banyaknya interval yang terbaik adalah 12

Secara Visual

pc12 <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,12), 
               lty = 1, col = "blue",se = F)+
  theme_bw()
pc12

1c. Data Auto Natural Cubic Spline

Menentukan Banyaknya Fungsi Basis

set.seed(123)
cross.val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 1:12
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)
                  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

basis

##   df     rmse     mae
## 1 10 4.260316 3.19201
## 2 10 4.260316 3.19201

Banyaknya fungsi basis terpilih ada 10. Menentukan Knot Berdasarkan Banyaknya Fungsi Basis Terpilih

df_10<-attr(ns(Auto$horsepower, df=10),"knots")
df_10

## [1]  67.0  72.0  80.0  88.0  93.5 100.0 110.0 140.0 157.7

Secara Empiris

set.seed(319)

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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
mean.metric.spline3.ns4 <- colMeans(metric.spline3.ns4)# menghitung rata-rata untuk 5 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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
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      mae
## NCS4 4 - NCS Knots 80, 120, 160, 200 4.353114 3.283937
## NCS9                 9 - NCS (df=10) 4.281758 3.205120

Secara empiris berdasarkan RMSE dan MAE, banyaknya knot terbaik adalah 9 knot, yaitu knot 67, 72, 80, 88, 93.5, 100, 110, 140, 157.7

Secara Visual

nc4 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(80, 120, 160, 200)), col = "blue", se = F) +
  theme_bw()


nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  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)+
  theme_bw()

plot_grid(nc4, nc9, labels = c("nc4", "nc9"), label_size = 10)

2. Perbandingan 3 Model Sebelumnya

(Regresi Polinomial, Piecewise Constant, dan Natural Cubic Spline (NCS))

Secara Empiris

perbandingan.model <- rbind(cv.poly3(Auto), cv.pc(Auto)[1,-1], model.ncs[2,-1])
perbandingan.model$metode <- c("Polinomial Ordo 3","Piecewise Constant 8","NCS 9 Knot")
perbandingan.model

##          rmse      mae               metode
## 1    4.320929 3.259139    Polinomial Ordo 3
## 2    4.438717 3.310112 Piecewise Constant 8
## NCS9 4.281758 3.205120           NCS 9 Knot

Model terbaik adalah NCS dengan 9 knot atau Natural Cubic Spline dengan 9 Knot

Secara Visual

poly3 <- ggplot(Auto,aes(x=horsepower, y=mpg), add=TRUE) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


pc8 <- ggplot(Auto,aes(x=horsepower, y=mpg), add=TRUE) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,8), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg), add=TRUE) + 
  geom_point(alpha = 0.55, color="red") + 
  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)+
  theme_bw()

plot_grid(poly2, pc8, nc9, labels = c("poly3","pc8", "nc9"), label_size = 6)

3. Data Amerika

a. Regresi Polinomial Data Amerika

Secara Empiris

cvpoly <- rbind (cv.poly1(AutoAmerika),cv.poly2(AutoAmerika),cv.poly3(AutoAmerika),cv.poly4(AutoAmerika))
cvpoly

##          rmse      mae
## [1,] 4.229279 3.350898
## [2,] 3.800834 2.894532
## [3,] 3.812732 2.915620
## [4,] 3.842147 2.917785

Secara Visual

poly1 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,1,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


poly2 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly3 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly4 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,4,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

plot_grid(poly1, poly2, poly3, poly4, labels=c("RegLin","poly2", "poly3", "poly4"), label_size = 6)

b. Piecewise Constant Data Amerika (Fungsi Tangga)

Secara Empiris

cv.pc(AutoAmerika)

##   breaks     rmse      mae
## 1      9 3.716899 2.856646
## 2     12 3.766258 2.850190

Secara empiris berdasarkan RMSE dan MAE banyaknya interval yang terbaik adalah 11

Secara Visual

pc9 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,9), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

pc11 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,11), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


plot_grid(pc9, pc11,  labels = c( "pc9", "pc11"), label_size = 6)

c. Natural Cubic Spline Data Amerika

Menentukan Banyaknya Basis

set.seed(123)
cross.val <- vfold_cv(AutoAmerika,v=10,strata = "mpg")
df <- 1:13
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)
                  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

basis

##   df     rmse      mae
## 1  7 3.769022 2.814088
## 2  7 3.769022 2.814088

diperoleh df atau banyaknya fungsi basis terpilih sebesar 8. Menentukan Knot Berdasarkan Fungsi Basis Terpilih

df_8<-attr(ns(AutoAmerika$horsepower, df=8),"knots")
df_8

## [1]  80  88  95 105 130 150 170

Secara Empiris

set.seed(405)

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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
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(80,88,95,105,130,150,170)),
         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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
mean.metric.spline3.ns7 <- colMeans(metric.spline3.ns7)# menghitung rata-rata untuk 5folds



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

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

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

model.ncs

##                                Model     rmse      mae
## NCS4 4 - NCS Knots 80, 120, 160, 200 3.788757 2.898205
## NCS7                  7 - NCS (df=8) 3.780835 2.819855

Secara empiris berdasarkan RMSE dan MAE, banyaknya knot terbaik adalah 7 knot.

Secara Visual

nc4 <- ggplot(AutoAmerika, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(80, 120, 160, 200)), col = "blue", se = F) +
  theme_bw()


nc7 <- ggplot(AutoAmerika, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(80,88,95,105,130,150,170)), col = "blue", se = F)+
  theme_bw()

plot_grid(nc4, nc7, labels = c("nc4", "nc7"), label_size = 8)

3a. Perbandingan Model untuk Data Amerika

Secara Empiris

perbandingan.model <- rbind(cv.poly2(AutoAmerika), cv.pc(AutoAmerika)[1,-1], model.ncs[2,-1])
perbandingan.model$metode <- c("Polinomial Ordo 2","Piecewise Constant 11","NCS 7 Knot")
perbandingan.model

##          rmse      mae                metode
## 1    3.800834 2.894532     Polinomial Ordo 2
## 2    3.716899 2.856646 Piecewise Constant 11
## NCS7 3.780835 2.819855            NCS 7 Knot

Model terbaik untuk data Auto Amerika adalah Piecewise Constant dengan 11 interval.

Secara Visual

poly2 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


pc11 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,11), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

nc7 <- ggplot(AutoAmerika, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(80,88,95,105,130,150,170)), col = "blue", se = F)+
  theme_bw()

plot_grid(poly2, pc11, nc7, labels = c("poly2","pc11", "nc7"), label_size = 6)

4. Data Eropa

a. Regresi Polinomial Data Eropa

Secara Empiris

set.seed(421)
cv.poly2 <- function(dataset){
  set.seed(421)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(421)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(421)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      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))
cvpoly

##          rmse      mae
## [1,] 4.882869 3.904732
## [2,] 4.917659 3.921019
## [3,] 5.003685 3.982591

Secara Visual

poly2 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly3 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly4 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,4,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

plot_grid(poly2, poly3, poly4, labels=c("poly2", "poly3", "poly4"), label_size = 6)

b. Piecewise Constant Data Eropa (Fungsi Tangga)

Secara Empiris

set.seed(424)

cv.pc2 <- function(dataset){
  set.seed(424)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)

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
                           )
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(breaks=breaks,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(rbind(basedonrmse,basedonmae))
}

cv.pc2(AutoEropa)

##   breaks     rmse      mae
## 1      5 5.022814 3.836602
## 2      5 5.022814 3.836602

Secara Visual

pc5 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,5), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

pc10 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,10), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


plot_grid(pc5, pc10,  labels = c( "pc5", "pc10"), label_size = 6)

c. Natural Cubic Spline Data Eropa

Menentukan Banyaknya Fungsi Basis

set.seed(430)
cross.val <- vfold_cv(AutoEropa,v=10,strata = NULL)

df <- 2:3
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)
                  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

basis

##   df     rmse      mae
## 1  2 5.049133 3.900407
## 2  2 5.049133 3.900407

Menentukan Knot Berdasarkan Banyaknya Fungsi Basis

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

## [1] 76.5

Secara Empiris

set.seed(433)

cross.val <- vfold_cv(AutoEropa,v=10,strata = NULL)

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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
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 80, 150")

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

model.ncs

##                      Model     rmse      mae
## NCS1    1 - NCS Knots 76.5 4.812342 3.877812
## NCS2 2 - NCS Knots 80, 150 4.833626 3.894311

Secara empiris berdasarkan RMSE dan MAE, banyaknya knot terbaik adalah 1 knot.

Secara Visual

nc1 <- ggplot(AutoEropa, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(76.5)), col = "blue", se = F) +
  theme_bw()


nc2 <- ggplot(AutoEropa, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(80,105)), col = "blue", se = F)+
  theme_bw()

plot_grid(nc1, nc2, labels = c("nc1", "nc2"), label_size = 8)

4a. Perbandingan Model untuk Data Eropa

Secara Empiris

perbandingan.model <- rbind(cv.poly2(AutoEropa), cv.pc2(AutoEropa)[1,-1], model.ncs[1,-1])
perbandingan.model$metode <- c("Polinomial Ordo 2","Piecewise Constant 5","NCS 1 Knot 105")
perbandingan.model

##          rmse      mae               metode
## 1    4.882869 3.904732    Polinomial Ordo 2
## 2    5.022814 3.836602 Piecewise Constant 5
## NCS1 4.812342 3.877812       NCS 1 Knot 105

Model terbaik untuk data Auto Eropa adalah Natural Cubic Spline (NCS) 1 Knot 105.

Secara Visual

poly2 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

pc5 <- ggplot(AutoEropa,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,5), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


nc1 <- ggplot(AutoEropa, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(76.5)), col = "blue", se = F) +
  theme_bw()


plot_grid(poly2, pc5,nc1, labels = c("poly2","pc5","nc1"), label_size = 6)

5. Data Jepang

a. Regresi Polinomial Data Jepang

Secara Empiris

set.seed(1021)

cv.poly2 <- function(dataset){
  set.seed(123)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly2<- colMeans(metric_poly2) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly2) 
}


cv.poly3 <- function(dataset){
  set.seed(1021)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      return(metric)
      }
  )
  poly3<- colMeans(metric_poly3) #menghitung rmse dan mae rata-rata untuk 10 folds
  return(poly3)
    }

cv.poly4 <- function(dataset){
  set.seed(1021)
  cross_val <- vfold_cv(dataset,v=10,strata = NULL)
  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)
      mae <- mlr3measures::mae(truth = truth,response = pred)
      metric <- c(rmse,mae) 
      names(metric) <- c("rmse","mae")
      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))
cvpoly

##          rmse      mae
## [1,] 4.524739 3.546178
## [2,] 3.818496 2.994819
## [3,] 4.349383 3.279423

Secara empiris jika dilihat dari RMSE dan MAE terlihat bahwa polinomial pangkat 3 memberikan nilai RMSE dan MAE yang terkecil sehingga dipilih polinomial pangkat 3.

Secara Visual

poly2 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

poly4 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,4,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

plot_grid(poly2, poly3, poly4, labels=c("poly2", "poly3", "poly4"), label_size = 6)

b. Piecewise Constant Data Jepang (Fungsi Tangga)

Secara Empiris

set.seed(1035)

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

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
                           )
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(breaks=breaks,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(rbind(basedonrmse,basedonmae))
}

cv.pc2(AutoJepang)

##   breaks     rmse      mae
## 1      3 4.205446 3.344684
## 2      2 4.319524 3.323871

Secara empiris berdasarkan RMSE banyaknya interval yang terbaik adalah 3

Secara Visual

pc2 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,2), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

pc3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,3), 
               lty = 1, col = "blue",se = F)+
  theme_bw()


plot_grid(pc2, pc3,  labels = c( "pc2", "pc3"), label_size = 6)

c. Natural Cubic Spline Data Jepang (NCS)

Menentukan Banyaknya Fungsi Basis

set.seed(123)
cross.val <- vfold_cv(AutoJepang,v=10,strata = NULL)

df <- 2:3
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)
                  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

basis

##   df     rmse      mae
## 1  3 4.162289 3.187387
## 2  3 4.162289 3.187387

Menentukan Knot Berdasarkan Fungsi Basis Terpilih

df_3<-attr(ns(AutoJepang$horsepower, df=3),"knots")
df_3

## [1] 68 90

Secara Empiris

set.seed(1051)

cross.val <- vfold_cv(AutoJepang,v=10,strata = NULL)

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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
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
                           )
mae <- mlr3measures::mae(truth = truth,
                           response = pred
                           )
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
mean.metric.spline3.ns3 <- colMeans(metric.spline3.ns3)# menghitung rata-rata untuk 5 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      mae
## NCS2       2 - NCS Knots 68,90 4.101217 3.239065
## NCS3 3 - NCS Knots 60, 88, 105 4.104026 3.201125

Secara empiris berdasarkan RMSE dan MAE, banyaknya knot terbaik adalah 2 knot.

Secara Visual

nc2 <- ggplot(AutoJepang, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(68,90)), col = "blue", se = F) +
  theme_bw()


nc3 <- ggplot(AutoJepang, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(60, 88, 105)), col = "blue", se = F)+
  theme_bw()

plot_grid(nc2, nc3, labels = c("nc2", "nc3"), label_size = 8)

5a. Perbandingan Model untuk Data Jepang

Secara Empiris

perbandingan.model <- rbind(cv.poly3(AutoJepang), cv.pc2(AutoJepang)[1,-1], model.ncs[1,-1])
perbandingan.model$metode <- c("Polinomial Ordo 3","Piecewise Constant 3","NCS 2 Knot 68,90")
perbandingan.model

##          rmse      mae               metode
## 1    3.818496 2.994819    Polinomial Ordo 3
## 2    4.205446 3.344684 Piecewise Constant 3
## NCS2 4.101217 3.239065     NCS 2 Knot 68,90

Model terbaik untuk data Auto Jepang adalah Model Regresi polinomial Ordo 3.

Secara Visual

poly3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "blue",se = F)+
  theme_bw()

pc3 <- ggplot(AutoJepang,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,3), 
               lty = 1, col = "blue",se = F)+
    theme_bw()

nc2 <- ggplot(AutoJepang, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(68,90)), col = "blue", se = F) +
  theme_bw()

plot_grid(poly3,pc3,nc2, labels=c("poly3","pc3","nc2"), label_size = 6)

6. Plot Model Terbaik

nc9 <- ggplot(Auto, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  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)+
  theme_bw()+
  labs(title = "NCS 9 Knots (df=10)", 
       subtitle = "Auto Data Set - World")

pc11 <- ggplot(AutoAmerika,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="red") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,11), 
               lty = 1, col = "blue",se = F)+
  theme_bw()+
  labs(title = "Piecewise Constant 11", 
       subtitle = "Auto Data Set - Amerika")

nc1 <- ggplot(AutoEropa, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(76.5)), col = "blue", se = F) +
  theme_bw()+
  labs(title = "Natural Cubic Splines Knots 76.5 (df=2)", 
       subtitle = "Auto Data Set - Eropa")

nc2 <- ggplot(AutoJepang, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.55, color="red") + 
  stat_smooth(method = "lm", formula = y~ns(x, knots = c(68,90)), col = "blue", se = F) +
  theme_bw()+
  labs(title = "Natural Cubic Splines Knots 86 90 (df=3)", 
       subtitle = "Auto Data Set - Jepang")


plot_grid(nc9, pc11, nc1, nc2)