Moving Beyond Linearity

Soal

Data yang digunakan diperoleh dari package ISLR. Gunakan model-model non-linear pada data Auto dengan respon mpg dan pilih tiga kolom untuk dijadikan peubah prediktor.

Apakah ada bukti untuk hubungan non-linear terhadap peubah-peubah yang dipilih? Buatlah visualisasi untuk mendukung jawaban Anda!

Tentukan model non-linear terbaik untuk masing-masing pasangan peubah (secara visual dan/atau secara empiris)

Data yang digunakan pada pembahasan kali ini adalah data Auto dari package ISLR yang terdiri dari 392 observasi dengan 9 peubah. Namun, pada pembahasan soal kali ini, peubah yang digunakan ada 4, yaitu:

mpg (mile per gallon)
acceleration
displacement
weight

Pemodelan acceleration dan mpg

Menampilkan data yang digunakan

datatable(cbind(mpg, acceleration))

plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Scatter Plot acceleration vs mpg")

Grafik di atas merupakan scatter plot dari 397 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah kecepatan suatu kendaraan (acceleration), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).

Secara grafik, terlihat bahwa pola data menunjukkan pola yang tidak linear. Untuk menemukan solusi optimal dari data kita bisa menggunakan Model Regresi Polinomial, Fungsi Tangga (Piecewise Constant) dan Natural Cubic Splines.

Regresi Polinomial

Secara Empiris

Cross Validation Data Auto menggunakan Model Regresi Polinomial

Regresi Polinomial Derajat 2

set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly2 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(acceleration,2,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly2)

Regresi Polinomial Derajat 3

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly3 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(acceleration,3,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly3)

Regresi Polinomial Derajat 4

    set.seed(1501211024)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(1501211024)
    metric_poly4 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(acceleration,4,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly4)

Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.

poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)

cvpoly <- rbind (poly2,poly3,poly4)
cvpoly

##           rmse      mae
## poly2 6.988011 5.748240
## poly3 7.042807 5.813189
## poly4 7.004028 5.688537

Berdasarkan hasil rata-rata rsme dan mae regresi polinomial diatas, terlihat bahwa Regresi Polinomial Derajat 3 memiliki nilai rsme dan mae paling kecil. Artinya, Regresi Polinomial Derajat 3 lebih baik dibandingkan Regresi Polinomial Derajat 2 dan 4.

# menampilkan model
nilai.cvpoly <- rbind("Poly3" = poly3)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 3")
model_poly_a <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_a

##                              Model     rmse      mae
## Poly3 Regresi Polinomial Derajat 3 7.042807 5.813189

Secara Visualisasi

par(mfrow=c(2,2))
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Scatter Plot acceleration vs mpg")

polinom_d2 <- lm(mpg ~ poly(acceleration,2))
prediksi <- predict(polinom_d2, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 2")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d3 <- lm(mpg ~ poly(acceleration,3))
prediksi <- predict(polinom_d3, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d4 <- lm(mpg ~ poly(acceleration,4))
prediksi <- predict(polinom_d4, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 4")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

Fungsi Tangga (Piecewise Constant)

Secara Empiris

Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

breaks <- 3:12
set.seed(179)
best_tangga <- map_dfr(breaks, function(i){
  metric_tangga <- map_dfr(cross_val$splits,
                           function(x){
                             training <- Auto[x$in_id,]
                             training$acceleration <- cut(training$acceleration,i)
                             
                             mod <- lm(mpg ~ acceleration,
                                       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 <- Auto[-x$in_id,]
                             
                             acceleration_new <- cut(testing$acceleration,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
                             
                             pred <- predict(mod,
                                             newdata=list(acceleration=acceleration_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)
                           }
  )
  
  metric_tangga
  
  # menghitung rata-rata untuk 8 folds
  mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)

# menampilkan semua jumlah interval
best_tangga

## # A tibble: 10 x 2
##     rmse   mae
##    <dbl> <dbl>
##  1  6.96  5.63
##  2  7.22  5.90
##  3  7.26  6.01
##  4  7.01  5.69
##  5  7.23  5.89
##  6  7.16  5.79
##  7  7.05  5.67
##  8  7.25  5.90
##  9  7.17  5.81
## 10  7.16  5.72

# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   breaks     rmse      mae
## 1      3 6.957824 5.626536
## 2      3 6.957824 5.626536

Berdasarkan rmse dan mae jumlah interval terbaik adalah 9, maka pada kasus ini breaks 9 atau Fungsi Tangga dengan 8 knots adalah interval terbaik. ### Secara Visualiasi

tangga9 <- lm(mpg ~ cut(acceleration,8))
prediksi <- predict(tangga9, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(acceleration, index.return=T)$ix
lines(acceleration[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 8 knots")
model_tangga_a <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_a

##                        Model     rmse      mae
## tangga Fungsi Tangga 8 knots 6.957824 5.626536

Regresi Natural Cubic Spline

Secara Empiris

Menentukan basis Model Regresi Natural Cubic Spline.

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

df <- 6:12
set.seed(179)
best_spline3 <- map_dfr(df, function(i){
  metric_spline3 <- map_dfr(cross_val$splits,
                            function(x){
                              mod <- lm(mpg ~ ns(acceleration,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
  
  ### menghitung rata-rata untuk 8 folds
  mean_metric_spline3 <- colMeans(metric_spline3)
  
  mean_metric_spline3
}
)

best_spline <- cbind(df=df,best_spline3)

# menampilkan hasil semua basis
best_spline

##   df     rmse      mae
## 1  6 7.028307 5.671402
## 2  7 7.048462 5.701358
## 3  8 7.035331 5.682875
## 4  9 7.060425 5.698845
## 5 10 7.091596 5.741730
## 6 11 7.102720 5.733232
## 7 12 7.092593 5.739242

# menampilkan basis berdasarkan rsme dan mae yang terpilih 
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   df     rmse      mae
## 1  6 7.028307 5.671402
## 2  6 7.028307 5.671402

Berdasarkan rsme basis yang terpilih adalah df=11, dan berdasarkan mae basis yang terpilih adalah df=12. Namun pada kasus ini basis yang dipilih berdasarkan rsme terkecil yaitu df=11. Selanjutnya, menampilkan vektor knots dengan df=11 dan 12.

attr(ns(Auto$acceleration, df=11), "knots")

## 9.090909% 18.18182% 27.27273% 36.36364% 45.45455% 54.54545% 63.63636% 72.72727% 
##  12.00000  13.20000  14.00000  14.50000  15.00000  15.62727  16.40000  17.00000 
## 81.81818% 90.90909% 
##  18.00000  19.20000

# menampilkan model
nilai.cv.ns <- rbind("NCS10" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("10 - NCS (df=11)")
model_ncs_a <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_a

##                  Model     rmse      mae
## NCS10 10 - NCS (df=11) 7.028307 5.671402

Secara Visualisasi

regnspline <- lm(mpg ~ ns(acceleration, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(acceleration), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(acceleration), interval =
## "predict"): prediction from a rank-deficient fit may be misleading

plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "NCS (df=11) acceleration vs mpg")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Perbandingan Model

Secara Empiris

rbind(model_poly_a, model_tangga_a, model_ncs_a)

##                               Model     rmse      mae
## Poly3  Regresi Polinomial Derajat 3 7.042807 5.813189
## tangga        Fungsi Tangga 8 knots 6.957824 5.626536
## NCS10              10 - NCS (df=11) 7.028307 5.671402

Secara empiris jika melihat nilai rsme dan mae, Model Natural Cubic Spline memiliki nilai rsme terkecil, Fungsi Tangga 8 knots memiliki nilai mae terkecil. Jika kita memilih Model terbaik berdasarkan nilai rsme terkecil maka Model Natural Cubic Spline adalah Model terbaik untuk data Auto untuk peubah acceleration vs mpg.

Secara Visualisasi

par(mfrow=c(2,2))

plot(displacement, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Scatter Plot acceleration vs mpg")

polinom_d3 <- lm(mpg ~ poly(acceleration,3))
prediksi <- predict(polinom_d3, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

tangga9 <- lm(mpg ~ cut(acceleration,8))
prediksi <- predict(tangga9, data.frame(acceleration), interval = "predict")
plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(acceleration, index.return=T)$ix
lines(acceleration[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

regnspline <- lm(mpg ~ ns(acceleration, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(acceleration), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(acceleration), interval =
## "predict"): prediction from a rank-deficient fit may be misleading

plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "NCS (df=11) acceleration vs mpg")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Sedangkan secara visualisasi, Model Fungsi Tangga 8 knots memiliki grafik yang lebih baik dibandingkan Model lainnya.

Pemodelan displacement dan mpg

Menampilkan Data yang Digunakan

datatable(cbind(mpg, displacement))

plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Scatter Plot displacement vs mpg")

Grafik di atas merupakan scatter plot dari 392 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah jarak tempuh suatu kendaraan (displacement), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).

Secara grafik, terlihat bahwa pola data menunjukkan ke pola nonlinear. Untuk menemukan solusi optimal dari data kita bisa menggunakan Model Regresi Polinomial, Fungsi Tangga (Piecewise Constant) dan Natural Cubic Splines.

Regresi Polinomial

###nSecara Empiris Cross Validation Data Auto menggunakan Model Regresi Polinomial - Regresi Polinomial Derajat 2

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly2 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(displacement,2,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly2)

Regresi Polinomial Derajat 3

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly3 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(displacement,3,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly3)

Regresi Polinomial Derajat 4

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly4 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(displacement,4,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly4)

Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.

poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)

cvpoly <- rbind (poly2,poly3,poly4)
cvpoly

##           rmse      mae
## poly2 4.316667 3.163100
## poly3 4.327419 3.177663
## poly4 4.344343 3.183435

# menampilkan model
nilai.cvpoly <- rbind("Poly3" = poly3)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 3")
model_poly_d <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_d

##                              Model     rmse      mae
## Poly3 Regresi Polinomial Derajat 3 4.327419 3.177663

Secara Visualisasi

par(mfrow=c(2,2))
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Scatter Plot displacement vs mpg")

polinom_d2 <- lm(mpg ~ poly(displacement,2))
prediksi <- predict(polinom_d2, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 2")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d3 <- lm(mpg ~ poly(displacement,3))
prediksi <- predict(polinom_d3, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d4 <- lm(mpg ~ poly(displacement,4))
prediksi <- predict(polinom_d4, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 4")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

Fungsi Tangga (Piecewise Constant)

Secara Empiris

Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

breaks <- 3:12
set.seed(179)
best_tangga <- map_dfr(breaks, function(i){
  metric_tangga <- map_dfr(cross_val$splits,
                           function(x){
                             training <- Auto[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 <- Auto[-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)
                           }
  )
  
  metric_tangga
  
  # menghitung rata-rata untuk 8 folds
  mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)

# menampilkan semua jumlah interval
best_tangga

## # A tibble: 10 x 2
##     rmse   mae
##    <dbl> <dbl>
##  1  4.91  3.70
##  2  4.83  3.63
##  3  4.66  3.49
##  4  4.59  3.41
##  5  4.56  3.37
##  6  4.40  3.25
##  7  4.25  3.12
##  8  4.28  3.14
##  9  4.43  3.29
## 10  4.48  3.32

# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   breaks     rmse      mae
## 1      9 4.253154 3.122767
## 2      9 4.253154 3.122767

Berdasarkan rmse dan mae jumlah interval terbaik adalah 9, maka pada kasus ini breaks 9 atau Fungsi Tangga dengan 8 knots adalah interval terbaik.

Secara Visualiasi

tangga9 <- lm(mpg ~ cut(displacement,8))
prediksi <- predict(tangga9, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(displacement, index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 8 knots")
model_tangga_d <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_d

##                        Model     rmse      mae
## tangga Fungsi Tangga 8 knots 4.253154 3.122767

Regresi Natural Cubic Spline

Secara Empiris

Menentukan basis Model Regresi Natural Cubic Spline.

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

df <- 6:12
set.seed(179)
best_spline3 <- map_dfr(df, function(i){
  metric_spline3 <- map_dfr(cross_val$splits,
                            function(x){
                              mod <- lm(mpg ~ ns(displacement,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
  
  ### menghitung rata-rata untuk 8 folds
  mean_metric_spline3 <- colMeans(metric_spline3)
  
  mean_metric_spline3
}
)

best_spline <- cbind(df=df,best_spline3)

# menampilkan hasil semua basis
best_spline

##   df     rmse      mae
## 1  6 4.289092 3.164622
## 2  7 4.259549 3.158275
## 3  8 4.253233 3.177852
## 4  9 4.236064 3.171636
## 5 10 4.211634 3.145508
## 6 11 4.170381 3.098630
## 7 12 4.161753 3.087814

# menampilkan basis berdasarkan rsme dan mae yang terpilih 
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   df     rmse      mae
## 1 12 4.161753 3.087814
## 2 12 4.161753 3.087814

attr(ns(Auto$displacement, df=11), "knots")

## 9.090909% 18.18182% 27.27273% 36.36364% 45.45455% 54.54545% 63.63636% 72.72727% 
##        90        97       107       120       140       168       231       258 
## 81.81818% 90.90909% 
##       318       351

# menampilkan model
nilai.cv.ns <- rbind("NCS10" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("10 - NCS (df=11)")
model_ncs_d <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_d

##                  Model     rmse      mae
## NCS10 10 - NCS (df=11) 4.161753 3.087814

Secara Visualisasi

regnspline <- lm(mpg ~ ns(displacement, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "NCS (df=11) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Perbandingan Model

Secara Empiris

rbind(model_poly_a, model_tangga_a, model_ncs_d)

##                               Model     rmse      mae
## Poly3  Regresi Polinomial Derajat 3 7.042807 5.813189
## tangga        Fungsi Tangga 8 knots 6.957824 5.626536
## NCS10              10 - NCS (df=11) 4.161753 3.087814

Secara Visualisasi

par(mfrow=c(2,2))

plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Scatter Plot displacement vs mpg")

polinom_d3 <- lm(mpg ~ poly(displacement,3))
prediksi <- predict(polinom_d3, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

tangga9 <- lm(mpg ~ cut(displacement,8))
prediksi <- predict(tangga9, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(displacement, index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

regnspline <- lm(mpg ~ ns(displacement, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")
plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "NCS (df=11) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Sedangkan secara visualisasi, Model Fungsi Tangga 8 knots memiliki grafik yang lebih baik dibandingkan Model lainnya.

Pemodelan weight dan mpg

Menampilkan Data yang Digunakan

attach(Auto)

## The following objects are masked from Auto (pos = 3):
## 
##     acceleration, cylinders, displacement, horsepower, mpg, name,
##     origin, weight, year

## The following object is masked from package:ggplot2:
## 
##     mpg

datatable(cbind(mpg, weight))

plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Scatter Plot weight vs mpg")

Grafik di atas merupakan scatter plot dari 392 kendaraan yang diukur berdasarkan dua peubah. Peubah pada sumbu x adalah jarak tempuh suatu kendaraan (weight), sedangkan peubah pada sumbu y adalah jarak yang ditempuh kendaraan untuk setiap galon bahan bakarnya (mile per gallon).

Regresi Polinomial

Secara Empiris

Cross Validation Data Auto menggunakan Model Regresi Polinomial

Regresi Polinomial Derajat 2

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly2 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(weight,2,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly2)

Regresi Polinomial Derajat 3

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly3 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(weight,3,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly3)

Regresi Polinomial Derajat 4

    set.seed(179)
    cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

    set.seed(179)
    metric_poly4 <- map_dfr(cross_val$splits,
                        function(x){
                          #data training
                          mod <- lm(mpg ~ poly(weight,4,raw = T),
                                    data=Auto[x$in_id,])
                          #data testing
                          pred <- predict(mod, newdata=Auto[-x$in_id,])
                          #nilai y pada data testing
                          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)
                        }
    )
    datatable(metric_poly4)

Pada Regresi Polinomial Derajat 2 sampai 4 menghasilkan rsme dan mae. Selanjutnya, dihitung rata-rata rsme dan mae setiap derajatnya untuk mengetahui rsme dan mae yang mana yang paling optimal.

poly2 <- colMeans(metric_poly2)
poly3 <- colMeans(metric_poly3)
poly4 <- colMeans(metric_poly4)

cvpoly <- rbind (poly2,poly3,poly4)
cvpoly

##           rmse      mae
## poly2 4.160664 3.064187
## poly3 4.172153 3.077568
## poly4 4.181472 3.084418

# menampilkan model
nilai.cvpoly <- rbind("Poly3" = poly3)
nama.model.cvpoly <- c("Regresi Polinomial Derajat 3")
model_poly_w <-data.frame(Model = nama.model.cvpoly, nilai.cvpoly)
model_poly_w

##                              Model     rmse      mae
## Poly3 Regresi Polinomial Derajat 3 4.172153 3.077568

Secara Visualisasi

par(mfrow=c(2,2))
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Scatter Plot weight vs mpg")

polinom_d2 <- lm(mpg ~ poly(weight,2))
prediksi <- predict(polinom_d2, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 2")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d3 <- lm(mpg ~ poly(weight,3))
prediksi <- predict(polinom_d3, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)


polinom_d4 <- lm(mpg ~ poly(weight,4))
prediksi <- predict(polinom_d4, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 4")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

Fungsi Tangga (Piecewise Constant)

Secara Empiris

Cross Validation Data Auto menggunakan Model Fungsi Tangga (Piecewise Constant)

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

breaks <- 3:12
set.seed(179)
best_tangga <- map_dfr(breaks, function(i){
  metric_tangga <- map_dfr(cross_val$splits,
                           function(x){
                             training <- Auto[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 <- Auto[-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)
                           }
  )
  
  metric_tangga
  
  # menghitung rata-rata untuk 8 folds
  mean_metric_tangga <- colMeans(metric_tangga)
})
best_tangga1 <- cbind(breaks=breaks,best_tangga)

# menampilkan jumlah interval terbaik berdasarkan rsme dan mae
bestrmse <- as.data.frame(best_tangga1 %>% slice_min(rmse))
bestmae <- as.data.frame(best_tangga1 %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   breaks     rmse      mae
## 1     12 4.239457 3.124184
## 2     12 4.239457 3.124184

Berdasarkan rmse dan mae jumlah interval terbaik adalah 9, maka pada kasus ini breaks 9 atau Fungsi Tangga dengan 8 knots adalah interval terbaik.

Secara Visualiasi

tangga9 <- lm(mpg ~ cut(weight,8))
prediksi <- predict(tangga9, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(weight, index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

# menampilkan model
nilai.tangga <- rbind("tangga" = as.data.frame(best_tangga %>% slice_min(rmse)))
nama.model.tangga <- c("Fungsi Tangga 8 knots")
model_tangga_w <-data.frame(Model = nama.model.tangga, nilai.tangga)
model_tangga_w

##                        Model     rmse      mae
## tangga Fungsi Tangga 8 knots 4.239457 3.124184

Regresi Natural Cubic Spline

Secara Empiris

Menentukan basis Model Regresi Natural Cubic Spline.

set.seed(179)
cross_val <- vfold_cv(Auto,v=8,strata = "mpg")

df <- 6:12
set.seed(179)
best_spline3 <- map_dfr(df, function(i){
  metric_spline3 <- map_dfr(cross_val$splits,
                            function(x){
                              mod <- lm(mpg ~ ns(weight,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
  
  ### menghitung rata-rata untuk 8 folds
  mean_metric_spline3 <- colMeans(metric_spline3)
  
  mean_metric_spline3
}
)

best_spline <- cbind(df=df,best_spline3)

# menampilkan hasil semua basis
best_spline

##   df     rmse      mae
## 1  6 4.169388 3.097331
## 2  7 4.147863 3.058274
## 3  8 4.160682 3.072711
## 4  9 4.189059 3.083213
## 5 10 4.181294 3.077955
## 6 11 4.176538 3.089075
## 7 12 4.175139 3.096803

# menampilkan basis berdasarkan rsme dan mae yang terpilih 
bestrmse <- as.data.frame(best_spline %>% slice_min(rmse))
bestmae <- as.data.frame(best_spline %>% slice_min(mae))
rbind(bestrmse,bestmae)

##   df     rmse      mae
## 1  7 4.147863 3.058274
## 2  7 4.147863 3.058274

attr(ns(Auto$weight, df=11), "knots")

## 9.090909% 18.18182% 27.27273% 36.36364% 45.45455% 54.54545% 63.63636% 72.72727% 
##  1985.000  2130.000  2264.636  2491.818  2676.364  2936.273  3228.364  3521.818 
## 81.81818% 90.90909% 
##  3896.545  4317.909

# menampilkan model
nilai.cv.ns <- rbind("NCS10" = as.data.frame(best_spline3 %>% slice_min(rmse)))
nama.model.ncs <- c("10 - NCS (df=11)")
model_ncs_w <-data.frame(Model = nama.model.ncs, nilai.cv.ns)
model_ncs_w

##                  Model     rmse      mae
## NCS10 10 - NCS (df=11) 4.147863 3.058274

Secara Visualisasi

regnspline <- lm(mpg ~ ns(weight, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(weight), interval = "predict"):
## prediction from a rank-deficient fit may be misleading

plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "NCS (df=11) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Perbandingan Model

Secara Empiris

rbind(model_poly_w, model_tangga_w, model_ncs_w)

##                               Model     rmse      mae
## Poly3  Regresi Polinomial Derajat 3 4.172153 3.077568
## tangga        Fungsi Tangga 8 knots 4.239457 3.124184
## NCS10              10 - NCS (df=11) 4.147863 3.058274

Secara Visualisasi

par(mfrow=c(2,2))

plot(displacement, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Scatter Plot weight vs mpg")

polinom_d3 <- lm(mpg ~ poly(weight,3))
prediksi <- predict(polinom_d3, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Regresi Polinomial Ordo 3")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix],
      main = "Polynomial Regression", col="blue", type="l",lwd=2)

tangga9 <- lm(mpg ~ cut(weight,8))
prediksi <- predict(tangga9, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(weight, index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l",lwd=2)

regnspline <- lm(mpg ~ ns(weight, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(weight), interval = "predict"):
## prediction from a rank-deficient fit may be misleading

plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "NCS (df=11) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)

Sedangkan secara visualisasi, Model Fungsi Tangga 8 knots memiliki grafik yang lebih baik dibandingkan Model lainnya.

Perbandingan semua model dan semua data

Secara Empiris

acceleration1 <- rbind(acceleration1 = rbind(model_poly_a, model_tangga_a, model_ncs_a) %>% slice_min(rmse) )

displacement1 <- rbind(displacement1 = rbind(model_poly_d, model_tangga_d, model_ncs_d) %>% slice_min(rmse))

weight1 <- rbind(weight1 = rbind(model_poly_w, model_tangga_w, model_ncs_w) %>% slice_min(rmse))

rbind(acceleration1,displacement1,weight1)

##                               Model     rmse      mae
## acceleration1 Fungsi Tangga 8 knots 6.957824 5.626536
## displacement1      10 - NCS (df=11) 4.161753 3.087814
## weight1            10 - NCS (df=11) 4.147863 3.058274

Bedasarkan perbandingan model dengan melihat rsme, Model Natural Cubic Spline adalah model terbaik untuk pada data Auto dengan peubah acceleration vs mpg, displacement vs mpg, dan weight vs mpg.

Secara Visualisasi

par(mfrow=c(2,2))

regnspline <- lm(mpg ~ ns(acceleration, knots=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7)))
prediksi <- predict(regnspline, data.frame(acceleration), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(acceleration), interval =
## "predict"): prediction from a rank-deficient fit may be misleading

plot(acceleration, mpg, pch=20, col="coral",
     xlab = "acceleration", ylab = "mile per gallon",
     main = "NCS (df=10) acceleration vs mpg")
ix <- sort(acceleration,index.return=T)$ix
lines(acceleration[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(67, 72, 80, 88, 93.5, 100, 110, 140, 157.7),col="grey50",lty=2)


regnspline <- lm(mpg ~ ns(displacement, knots=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286)))
prediksi <- predict(regnspline, data.frame(displacement), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(displacement), interval =
## "predict"): prediction from a rank-deficient fit may be misleading

plot(displacement, mpg, pch=20, col="coral",
     xlab = "displacement", ylab = "mile per gallon",
     main = "NCS (df=7) displacement vs mpg")
ix <- sort(displacement,index.return=T)$ix
lines(displacement[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(2074.857, 2285.429, 2635, 2986.571, 3446.143, 4096.286),col="grey50",lty=2)


tangga9 <- lm(mpg ~ cut(weight,8))
prediksi <- predict(tangga9, data.frame(weight), interval = "predict")
plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "Fungsi Tangga dengan 8 knots")
ix <- sort(weight, index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l",lwd=2)


regnspline <- lm(mpg ~ ns(weight, knots=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351)))
prediksi <- predict(regnspline, data.frame(weight), interval = "predict")

## Warning in predict.lm(regnspline, data.frame(weight), interval = "predict"):
## prediction from a rank-deficient fit may be misleading

plot(weight, mpg, pch=20, col="coral",
     xlab = "weight", ylab = "mile per gallon",
     main = "NCS (df=11) weight vs mpg")
ix <- sort(weight,index.return=T)$ix
lines(weight[ix], prediksi[ix], col="blue",
      type="l", lwd=2)
abline(v=c(90, 97, 107, 120, 140, 168, 231, 258, 318, 351), col="grey50",lty=2)