Regresi Nonlinear pada Data Informasi Kendaraan

Library

Terdapat 9 packages yang digunakan, yaitu:
ISLR, MultiKink, tidyverse, ggplot2, dplyr, purrr, rsample, splines, dan cowplot

Data

Dataset yang digunakan adalah Auto data set dari library ISLR, yaitu data yang memuat informasi tentang 392 kendaraan yang terdiri atas 9 peubah. Pada analisis regresi nonlinear ini, digunakan dua peubah:

mpg (miles per gallon): jarak yang ditempuh kendaraan per galon bahan bakar (peubah respon)
horsepower (engine horsepower): kekuatan mesin (peubah prediktor)

Pada analisis ini, akan ditentukan model nonlinear terbaik untuk pasangan peubah mpg dan horsepower.

library(ISLR)
autodata = Auto %>% select(mpg, horsepower)
tibble(autodata)

## # A tibble: 392 × 2
##      mpg horsepower
##    <dbl>      <dbl>
##  1    18        130
##  2    15        165
##  3    18        150
##  4    16        150
##  5    17        140
##  6    15        198
##  7    14        220
##  8    14        215
##  9    14        225
## 10    15        190
## # … with 382 more rows

y <- autodata$mpg
x <- autodata$horsepower

Eksplorasi Data

ggplot(autodata, aes(x = horsepower, y = mpg)) +
  geom_point(alpha = 2, color = "#966989") + 
  ggtitle("Scatterplot Hubungan mpg dan horsepower") +
  theme_grey() +
  theme(plot.title = element_text(hjust = 0.5))

Berdasarkan scatterplot di atas, terlihat bahwa terdapat hubungan negatif antara mpg dan horsepower. Artinya, semakin besar kekuatan mesin (horsepower), jarak yang ditempuh kendaraan per galon bahan bakarnya akan semakin kecil. Dalam arti lain, semakin kuat mesin akan mengakibatkan bahan bakar kendaraan semakin boros. Pola data menyebar melengkung sehingga dapat dikatakan hubungan antara kedua peubah bersifat tidak linear.

Regresi Linier

lin.mod1 <- lm(mpg~horsepower, data = autodata)
summary(lin.mod1)

## 
## Call:
## lm(formula = mpg ~ horsepower, data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5710  -3.2592  -0.3435   2.7630  16.9240 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 39.935861   0.717499   55.66   <2e-16 ***
## horsepower  -0.157845   0.006446  -24.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.906 on 390 degrees of freedom
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.6049 
## F-statistic: 599.7 on 1 and 390 DF,  p-value: < 2.2e-16

Berdasarkan output di atas, seluruh parameter signifikan pada taraf nyata 5%. Akan tetapi, terlihat bahwa hubungan linear tidak tepat digunakan karena akan ada kejadian di mana kendaraan dengan kekuatan mesin yang sangat besar atau lebih besar dari data akan menempuh jarak nol atau negatif per galon bahan bakarnya. Hal ini sejalan dengan plot di bawah dan tidak masuk akal untuk terjadi.

ggplot(autodata, aes(x = horsepower, y = mpg)) +
                 geom_point(alpha=0.5, color="#966989") +
  stat_smooth(method = "lm", 
               formula = y~x, 
               lty = 1, col = "maroon", se = F)+
  theme_grey() +
  theme(plot.title = element_text(hjust = 0.5))

Karena terdapat hubungan nonlinier antara peubah mpg dan horsepower, untuk memodelkan kedua peubah ini dapat menggunakan regresi nonlinear, seperti regresi polinomial, fungsi tangga (piecewise), atau metode splines.

Regresi Polinomial

Regresi polinomial merupakan regresi linier berganda yang dibentuk dengan menjumlahkan pengaruh variabel prediktor (X) yang dipangkatkan secara meningkat sampai orde ke-k. Struktur analisis regresi polinomial sama dengan model regresi linier berganda.

Ordo = 2

pol.mod2 = lm(mpg ~ poly(horsepower, 2, raw = T), data = autodata)
summary(pol.mod2)

## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = T), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7135  -2.5943  -0.0859   2.2868  15.8961 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   56.9000997  1.8004268   31.60   <2e-16 ***
## poly(horsepower, 2, raw = T)1 -0.4661896  0.0311246  -14.98   <2e-16 ***
## poly(horsepower, 2, raw = T)2  0.0012305  0.0001221   10.08   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared:  0.6876, Adjusted R-squared:  0.686 
## F-statistic:   428 on 2 and 389 DF,  p-value: < 2.2e-16

Berdasarkan summary model di atas, parameter pada model polinomial dengan ordo = 2 seluruhnya signifikan pada alpha = 5%.

Ordo = 3

pol.mod3 = lm(mpg ~ poly(horsepower, 3, raw = T), data = autodata)
summary(pol.mod3)

## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7039  -2.4491  -0.1519   2.2035  15.8159 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    6.068e+01  4.563e+00  13.298  < 2e-16 ***
## poly(horsepower, 3, raw = T)1 -5.689e-01  1.179e-01  -4.824 2.03e-06 ***
## poly(horsepower, 3, raw = T)2  2.079e-03  9.479e-04   2.193   0.0289 *  
## poly(horsepower, 3, raw = T)3 -2.147e-06  2.378e-06  -0.903   0.3673    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.375 on 388 degrees of freedom
## Multiple R-squared:  0.6882, Adjusted R-squared:  0.6858 
## F-statistic: 285.5 on 3 and 388 DF,  p-value: < 2.2e-16

Berdasarkan summary model di atas, terdapat satu parameter yang tidak signifikan pada model polinomial dengan ordo = 3 untuk alpha = 5%.

poly2 <- ggplot(autodata, aes(x=horsepower, y=mpg)) + 
  geom_point(alpha = 0.5, color = "#966989") +
  stat_smooth(method = "lm", 
               formula = y~poly(x, 2, raw=T), 
               lty = 1, col = "navy", se = F, lwd = 0.7) +
  theme_grey() +
  ggtitle("Regresi Polinomial Derajat 2") +
  theme(plot.title = element_text(hjust = 0.5))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

poly3 <- ggplot(autodata, aes(x=horsepower, y=mpg)) + 
  geom_point(alpha = 0.5, color = "#966989") +
  stat_smooth(method = "lm", 
               formula = y~poly(x, 3, raw=T), 
               lty = 1, col = "navy", se = F, lwd = 0.7) +
  theme_grey() +
  ggtitle("Regresi Polinomial Derajat 3") +
  theme(plot.title = element_text(hjust = 0.5))

plot_grid(poly2, poly3)

Pada plot di atas, garis regresi polinomial dengan derajat 2 dan derajat 3 tidak berbeda signifikan, keduanya cenderung melengkung mengikuti pola sebaran dari data.

Fungsi Tangga

Suatu fungsi dikatakan piecewise constant apabila fungsi tersebut adalah konstanta lokal di daerah-daerah terhubung yang dipisahkan oleh batas-batas dimensi bawah yang mungkin jumlahnya tak terhingga.

set.seed(123)
cross_val <- vfold_cv(Auto, v=10, strata = "mpg")
breaks <- 3:10
best_tangga <- map_dfr(breaks, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
    function(x){
        training <- Auto[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 <- Auto[-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)
      }
    )
  metric_tangga
  # menghitung rata-rata untuk 10 folds
  mean_metric_tangga <- colMeans(metric_tangga)
  mean_metric_tangga
  }
)

best_tangga <- cbind(breaks=breaks,best_tangga)
# menampilkan hasil all breaks
best_tangga

##   breaks     rmse      mae
## 1      3 5.775792 4.521829
## 2      4 4.985270 3.774549
## 3      5 4.712623 3.571614
## 4      6 4.644383 3.548375
## 5      7 4.554116 3.379980
## 6      8 4.437597 3.405433
## 7      9 4.549668 3.471999
## 8     10 4.589172 3.455531

# Berdasarkan RMSE
best_tangga %>% slice_min(rmse)

##   breaks     rmse      mae
## 1      8 4.437597 3.405433

Jumlah breaks yang terbaik berdasarkan RMSE adalah 8 breaks.

# Berdasarkan MAE
best_tangga %>% slice_min(mae)

##   breaks     rmse     mae
## 1      7 4.554116 3.37998

Jumlah breaks yang terbaik berdasarkan MAE adalah 7 breaks.

Breaks = 7

step.mod7 <- lm(mpg~cut(horsepower, 7), data = autodata)
summary(step.mod7)

## 
## Call:
## lm(formula = mpg ~ cut(horsepower, 7), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.5280  -2.5530  -0.3758   2.2881  16.7482 
## 
## Coefficients:
##                               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    32.5280     0.5069   64.17   <2e-16 ***
## cut(horsepower, 7)(72.3,98.6]  -6.8162     0.6358  -10.72   <2e-16 ***
## cut(horsepower, 7)(98.6,125]  -12.1523     0.7591  -16.01   <2e-16 ***
## cut(horsepower, 7)(125,151]   -16.5763     0.7957  -20.83   <2e-16 ***
## cut(horsepower, 7)(151,177]   -18.5020     1.0831  -17.08   <2e-16 ***
## cut(horsepower, 7)(177,204]   -19.4447     1.4187  -13.71   <2e-16 ***
## cut(horsepower, 7)(204,230]   -19.6280     1.5375  -12.77   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.59 on 385 degrees of freedom
## Multiple R-squared:  0.6594, Adjusted R-squared:  0.6541 
## F-statistic: 124.2 on 6 and 385 DF,  p-value: < 2.2e-16

Model fungsi tangga dengan breaks = 7 menghasilkan seluruh parameter signifikan pada taraf nyata 5%.

Breaks = 8

step.mod8 <- lm(mpg~cut(horsepower, 8), data = autodata)
summary(step.mod8)

## 
## Call:
## lm(formula = mpg ~ cut(horsepower, 8), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.9471  -2.6757  -0.1533   2.4015  14.5529 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  33.9085     0.5771   58.76   <2e-16 ***
## cut(horsepower, 8)(69,92]    -6.9614     0.6910  -10.07   <2e-16 ***
## cut(horsepower, 8)(92,115]  -12.7551     0.7425  -17.18   <2e-16 ***
## cut(horsepower, 8)(115,138] -15.6656     1.1264  -13.91   <2e-16 ***
## cut(horsepower, 8)(138,161] -18.7799     0.8568  -21.92   <2e-16 ***
## cut(horsepower, 8)(161,184] -19.9735     1.1470  -17.41   <2e-16 ***
## cut(horsepower, 8)(184,207] -21.1228     1.7720  -11.92   <2e-16 ***
## cut(horsepower, 8)(207,230] -21.0085     1.5159  -13.86   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.433 on 384 degrees of freedom
## Multiple R-squared:  0.6832, Adjusted R-squared:  0.6774 
## F-statistic: 118.3 on 7 and 384 DF,  p-value: < 2.2e-16

Model fungsi tangga dengan breaks = 8 menghasilkan seluruh parameter signifikan pada taraf nyata 5%.

int.8 <- ggplot(autodata, aes(x = horsepower, y = mpg)) +
                 geom_point(alpha=0.5, color="#966989") +
  stat_smooth(method = "lm", 
               formula = y~cut(x, 8), 
               lty = 1, col = "#C21460", se = F)+
  theme_grey() +
  ggtitle("Fungsi Tangga dengan Breaks = 8") +
  theme(plot.title = element_text(hjust = 0.5))

int.7 <- ggplot(autodata,aes(x = horsepower, y = mpg)) +
                 geom_point(alpha=0.5, color="#966989") +
  stat_smooth(method = "lm", 
               formula = y~cut(x, 7), 
               lty = 1, col = "#C21460", se = F)+
  theme_grey() +
  ggtitle("Fungsi Tangga dengan Breaks = 7") +
  theme(plot.title = element_text(hjust = 0.5))

plot_grid(int.7, int.8)

Berdasarkan plot di atas, pola garis pada fungsi tangga dengan breaks 7 ataupun 8 sangat mirip atau hampir tidak terlihat perbedaan selain jumlah interval yang digunakan.

Perbandingan Model

Model1 <- c("Linear", "Polinomial (ordo = 2)", "Polinomial (ordo = 3)", "Piecewise (breaks = 7)", "Piecewise (breaks = 8")
AIC1 <- c(AIC(lin.mod1), AIC(pol.mod2), AIC(pol.mod3), AIC(step.mod7), AIC(step.mod8))
sig1 <- c("2 dari 2", "3 dari 3", "3 dari 4", "7 dari 7", "8 dari 8")
tabel1 <- data.frame(Model1, AIC1, sig1)
kableExtra::kable(tabel1, align = c(rep('c', times = 2)), 
                  col.names = c("Nama Model", "AIC", "Parameter Signifikan"))

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 5 > 1' in coercion to 'logical(1)'

Nama Model	AIC	Parameter Signifikan
Linear	2363.324	2 dari 2
Polinomial (ordo = 2)	2274.354	3 dari 3
Polinomial (ordo = 3)	2275.531	3 dari 4
Piecewise (breaks = 7)	2316.152	7 dari 7
Piecewise (breaks = 8	2289.761	8 dari 8

Berdasarkan tabel di atas, nilai AIC minimum adalah 2274.354, yang dimiliki oleh model yang menggunakan regresi polinomial dengan ordo = 2. Selain itu, seluruh parameter signifikan pada taraf nyata 5%. Dengan demikian, model regresi polinomial dengan ordo = 2 merupakan model yang paling baik di antara lima model yang telah dicobakan. Plot dari model tersebut adalah sebagai berikut:

poly2

Natural Cubic Spline

nspline = function(dataset){
  set.seed(101)
  cross.val <- vfold_cv(dataset,v=10,strata = "mpg")
  df <- 1:20
  best.spline3 <- map_dfr(df, function(i){
    metric.spline3 <- map_dfr(cross.val$splits,
                              function(x){
                                mod <- lm(mpg ~ ns(horsepower,df=i),data=dataset[x$in_id,])
                                pred <- predict(mod,newdata=dataset[-x$in_id,])
                                truth <- dataset[-x$in_id,]$mpg
                                rmse <- mlr3measures::rmse(truth = truth,response = pred)
                                mae <- mlr3measures::mae(truth = truth,response = pred)
                                metric <- c(rmse,mae)
                                names(metric) <- c("rmse","mae")
                                return(metric)
                              }
    )
    metric.spline3
    mean.metric.spline3 <- colMeans(metric.spline3) #rata-rata untuk 10 folds
    mean.metric.spline3
  }
  )
  best.spline3 <- cbind(df=df,best.spline3)
  basis<-data.frame(rbind(best.spline3 %>% slice_min(rmse), #berdasarkan rmse
                          best.spline3 %>% slice_min(mae))) #berdasarkan mae

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

basisAuto = nspline(autodata)
basisAuto

## $basis
##   df     rmse      mae
## 1 10 4.249102 3.184156
## 2 12 4.275312 3.182377
## 
## $knot1
##   10%   20%   30%   40%   50%   60%   70%   80%   90% 
##  67.0  72.0  80.0  88.0  93.5 100.0 110.0 140.0 157.7 
## 
## $knot2
## 8.333333% 16.66667%       25% 33.33333% 41.66667%       50% 58.33333% 66.66667% 
##     65.00     70.00     75.00     84.00     88.00     93.50    100.00    110.00 
##       75% 83.33333% 91.66667% 
##    126.00    150.00    165.83 
## 
## $length
## [1] "length knot1 = 9"  "length knot2 = 11"

Berdasarkan nilai RMSE, dipilih 9 knot (df = 10). Berdasarkan nilai MAE dipilih 11 knot (df = 12).

k1 = basisAuto$knot1
k2 = basisAuto$knot2

NCS 9 Knot

model.ncs9 <- lm(mpg~ns(horsepower, knots = k1), data = autodata)
summary(model.ncs9)

## 
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = k1), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5003  -2.4956  -0.2086   2.2577  14.6400 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    32.462      1.680  19.324  < 2e-16 ***
## ns(horsepower, knots = k1)1    -4.600      1.789  -2.572 0.010495 *  
## ns(horsepower, knots = k1)2    -3.007      2.555  -1.177 0.239996    
## ns(horsepower, knots = k1)3    -8.663      2.043  -4.239 2.82e-05 ***
## ns(horsepower, knots = k1)4    -7.170      2.188  -3.278 0.001142 ** 
## ns(horsepower, knots = k1)5   -14.722      2.125  -6.927 1.84e-11 ***
## ns(horsepower, knots = k1)6    -8.446      2.467  -3.424 0.000684 ***
## ns(horsepower, knots = k1)7   -16.947      2.080  -8.148 5.38e-15 ***
## ns(horsepower, knots = k1)8   -21.715      1.891 -11.483  < 2e-16 ***
## ns(horsepower, knots = k1)9   -14.384      4.116  -3.494 0.000531 ***
## ns(horsepower, knots = k1)10  -22.326      1.946 -11.476  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.24 on 381 degrees of freedom
## Multiple R-squared:  0.7124, Adjusted R-squared:  0.7049 
## F-statistic: 94.39 on 10 and 381 DF,  p-value: < 2.2e-16

Untuk model NCS dengan 9 knot, terdapat satu parameter yang tidak signifikan dari keseluruhan 11 parameter pada taraf nyata 5%.

NCS 11 Knot

model.ncs11 <- lm(mpg~ns(horsepower, knots = k2), data = autodata)
summary(model.ncs11)

## 
## Call:
## lm(formula = mpg ~ ns(horsepower, knots = k2), data = autodata)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.9346  -2.4576  -0.0912   2.4092  15.2377 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    32.890      1.770  18.585  < 2e-16 ***
## ns(horsepower, knots = k2)1    -1.614      1.840  -0.878 0.380727    
## ns(horsepower, knots = k2)2    -5.920      2.580  -2.295 0.022298 *  
## ns(horsepower, knots = k2)3    -4.144      2.442  -1.697 0.090576 .  
## ns(horsepower, knots = k2)4    -9.549      2.183  -4.375 1.57e-05 ***
## ns(horsepower, knots = k2)5    -6.790      2.325  -2.920 0.003709 ** 
## ns(horsepower, knots = k2)6   -16.160      2.386  -6.772 4.86e-11 ***
## ns(horsepower, knots = k2)7    -9.387      2.382  -3.941 9.68e-05 ***
## ns(horsepower, knots = k2)8   -15.476      2.679  -5.777 1.59e-08 ***
## ns(horsepower, knots = k2)9   -17.894      2.134  -8.383 1.02e-15 ***
## ns(horsepower, knots = k2)10  -21.849      2.034 -10.740  < 2e-16 ***
## ns(horsepower, knots = k2)11  -16.455      4.489  -3.666 0.000282 ***
## ns(horsepower, knots = k2)12  -21.675      2.060 -10.520  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.251 on 379 degrees of freedom
## Multiple R-squared:  0.7125, Adjusted R-squared:  0.7034 
## F-statistic: 78.27 on 12 and 379 DF,  p-value: < 2.2e-16

Untuk model NCS dengan 11 knot, terdapat dua parameter yang tidak signifikan dari keseluruhan 13 parameter pada taraf nyata 5%.

ncs9 <- ggplot(autodata, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.5, color = "#966989") + 
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = k1), 
              col = "#B2D732", se = F, lwd = 0.8) +
  theme_grey()+
  ggtitle("Regresi Natural Cubic Spline 9 Knot") +
  theme(plot.title = element_text(hjust = 0.5))

ncs11 <- ggplot(autodata, aes(x = horsepower, y = mpg)) + 
  geom_point(alpha = 0.5, color = "#966989") + 
  stat_smooth(method = "lm", 
              formula = y~ns(x, knots = k2), 
              col = "#B2D732", se = F, lwd = 0.8) +
  theme_grey()+
  ggtitle("Regresi Natural Cubic Spline 11 Knot") +
  theme(plot.title = element_text(hjust = 0.5))

plot_grid(ncs9, ncs11)

Pada plot di atas, terlihat bahwa garis pada model dengan 9 knot menyebar dengan pola yang mirip dengan garis pada model dengan 11 knot.

Perbandingan Model

Model2 <- c("NCS Knot = 9", "NCS Knot = 11")
AIC2 <- c(AIC(model.ncs9), AIC(model.ncs11))
sig2 <- c("10 dari 11", "11 dari 13")
tabel2 <- data.frame(Model2, AIC2, sig2)
kableExtra::kable(tabel2, align = c(rep('c', times = 2)), 
                  col.names = c("Nama Model", "AIC", "Parameter Signifikan"))

Nama Model	AIC	Parameter Signifikan
NCS Knot = 9	2257.847	10 dari 11
NCS Knot = 11	2261.762	11 dari 13

Berdasarkan nilai AIC-nya, dapat dilihat bahwa model Natural Cubic Spline dengan 9 knot memiliki AIC minimum. Selain itu, jumlah parameter yang signifikan adalah 10 dari 11 parameter (terdapat satu parameter yang tidak signifikan pada taraf nyata 5%). Oleh karena itu, model dengan 9 knot lebih baik digunakan pada data yang dimiliki.

Smoothing Spline

Regresi smoothing spline menduga fungsi regresi dengan meminimumkan kriteria kuadrat terkecil tersanksi (penalized least squares). Pendugaan terhadap fungsi f bergantung kepada parameter pemulus λ yang nilainya secara umum diperoleh dari data.

model_sms <- with(data = autodata, 
                  smooth.spline(horsepower, mpg))
model_sms

## Call:
## smooth.spline(x = horsepower, y = mpg)
## 
## Smoothing Parameter  spar= 0.2648644  lambda= 5.101567e-07 (12 iterations)
## Equivalent Degrees of Freedom (Df): 42.14987
## Penalized Criterion (RSS): 1117.657
## GCV: 18.64132

Berdasarkan output di atas, parameter lambda didapatkan menggunakan 12 kali iterasi. Selain lambda, didapatkan parameter spar yang digunakan untuk mengatur tingkat smoothness dari kurva, yaitu spar = 0.26486.

Fungsi smooth.spline secara otomatis memilih parameter lambda terbaik dengan menggunakan Cross-Validation (CV) dan Generalized Cross-Validation (GCV). Perbedaan utama antara keduanya adalah pada ukuran kebaikan model yang digunakan. CV menggunakan MSE sebagai ukuran kebaikan model, sedangkan GCV menggunakan Weighted-MSE. Banyaknya folds (lipatan) pada CV dan GCV adalah sejumlah observasinya.

pred_data <- broom::augment(model_sms)

ggplot(pred_data,aes(x=x, y=y))+
  geom_point(alpha = 0.5, color = "#966989") +
  geom_line(aes(y=.fitted),col = "#5f8cb9",
            lwd = 1) +
  xlab("horsepower") +
  ylab("mpg") +
  theme_grey()

Pengaruh parameter lambda terhadap tingkat smoothness kurva bisa dilihat dari ilustrasi berikut:

model_sms_lambda <- data.frame(lambda = seq(0, 5, by = 0.5)) %>% 
  group_by(lambda) %>% 
  do(broom::augment(with(data = autodata,
                         smooth.spline(horsepower, mpg, lambda = .$lambda))))

p <- ggplot(model_sms_lambda,
       aes(x = x, y = y)) +
  geom_line(aes(y = .fitted),
            col = "#5f8cb9",
            lty = 1) +
  facet_wrap(~lambda)
p

Jika kita tentukan df = 7, maka hasil kurva model smooth.spline akan lebih merepresentasikan data.

model_sms1 <- with(data = autodata, 
                   smooth.spline(horsepower, mpg, df = 7))
model_sms1

## Call:
## smooth.spline(x = horsepower, y = mpg, df = 7)
## 
## Smoothing Parameter  spar= 0.7927002  lambda= 0.003320238 (12 iterations)
## Equivalent Degrees of Freedom (Df): 6.999082
## Penalized Criterion (RSS): 2424.819
## GCV: 18.84972

pred_data <- broom::augment(model_sms1)

ggplot(pred_data, aes(x = x, y = y)) +
  geom_point(alpha = 0.5, color = "#966989") +
  geom_line(aes(y = .fitted), col = "#5f8cb9", lwd = 1) +
  xlab("horsepower") +
  ylab("mpg") +
  theme_grey()

LOESS

Metode LOESS merupakan metode pemulus yang secara luas sering digunakan dengan sifat-sifat kekekaran (robustness) yang baik. LOESS merupakan running-line smoother terboboti dan setiap garis lokal diduga mengunakan suatu metode robust menggantikan kuadrat terkecil.

model_loess <- loess(mpg~horsepower,
                     data = autodata)
summary(model_loess)

## Call:
## loess(formula = mpg ~ horsepower, data = autodata)
## 
## Number of Observations: 392 
## Equivalent Number of Parameters: 4.97 
## Residual Standard Error: 4.32 
## Trace of smoother matrix: 5.43  (exact)
## 
## Control settings:
##   span     :  0.75 
##   degree   :  2 
##   family   :  gaussian
##   surface  :  interpolate      cell = 0.2
##   normalize:  TRUE
##  parametric:  FALSE
## drop.square:  FALSE

Pengaruh parameter span terhadap tingkat smoothness kurva bisa dilihat dari ilustrasi berikut:

model_loess_span <- data.frame(span = seq(0.1, 5, by = 0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~horsepower,
                     data = autodata, span = .$span)))

p2 <- ggplot(model_loess_span,
       aes(x = horsepower, y = mpg)) +
  geom_line(aes(y = .fitted),
            col = "#5f8cb9",
            lty = 1) +
  facet_wrap(~span)
p2

Pendekatan ini juga dapat dilakukan dengan fungsi stat_smooth() pada package ggplot2

library(ggplot2)
ggplot(autodata, aes(horsepower, mpg)) +
  geom_point(alpha = 0.5, color = "#966989") +
  stat_smooth(method = 'loess',
              formula = y~x,
              span = 0.75, col = "orange", lty = 1, se = F, lwd = 0.8)

Tuning span dapat dilakukan dengan menggunakan cross-validation:

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

span <- seq(0.1,1,length.out=50)

best_loess <- map_dfr(span, function(i){
metric_loess <- map_dfr(cross_val$splits,
    function(x){
mod <- loess(mpg ~ horsepower,span = i,
         data = Auto[x$in_id,])
pred <- predict(mod,
               newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg

data_eval <- na.omit(data.frame(pred=pred,
                                truth=truth))


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

head(metric_loess,20)

# menghitung rata-rata untuk 10 folds
mean_metric_loess <- colMeans(metric_loess)

mean_metric_loess
}
)

best_loess <- cbind(span=span,best_loess)
# menampilkan hasil all breaks
best_loess

##         span     rmse      mae
## 1  0.1000000 4.279961 3.212481
## 2  0.1183673 4.279930 3.199382
## 3  0.1367347 4.282181 3.206301
## 4  0.1551020 4.271935 3.206052
## 5  0.1734694 4.280129 3.216263
## 6  0.1918367 4.291612 3.230315
## 7  0.2102041 4.304326 3.241927
## 8  0.2285714 4.306384 3.242725
## 9  0.2469388 4.299509 3.226658
## 10 0.2653061 4.294184 3.215434
## 11 0.2836735 4.295722 3.222993
## 12 0.3020408 4.289119 3.217525
## 13 0.3204082 4.292602 3.222667
## 14 0.3387755 4.284627 3.213171
## 15 0.3571429 4.284646 3.210649
## 16 0.3755102 4.286424 3.214093
## 17 0.3938776 4.286272 3.214380
## 18 0.4122449 4.290260 3.218785
## 19 0.4306122 4.295797 3.225056
## 20 0.4489796 4.291814 3.224207
## 21 0.4673469 4.293737 3.225118
## 22 0.4857143 4.298144 3.228398
## 23 0.5040816 4.301433 3.232241
## 24 0.5224490 4.302167 3.232329
## 25 0.5408163 4.303161 3.232579
## 26 0.5591837 4.304828 3.234132
## 27 0.5775510 4.307357 3.236528
## 28 0.5959184 4.308655 3.237146
## 29 0.6142857 4.307372 3.236599
## 30 0.6326531 4.306337 3.235772
## 31 0.6510204 4.305561 3.234744
## 32 0.6693878 4.306828 3.235237
## 33 0.6877551 4.307049 3.236061
## 34 0.7061224 4.308953 3.237984
## 35 0.7244898 4.310339 3.239400
## 36 0.7428571 4.318304 3.244223
## 37 0.7612245 4.324715 3.247615
## 38 0.7795918 4.330090 3.250112
## 39 0.7979592 4.331711 3.249172
## 40 0.8163265 4.335098 3.249201
## 41 0.8346939 4.338779 3.250354
## 42 0.8530612 4.338608 3.249875
## 43 0.8714286 4.339232 3.249140
## 44 0.8897959 4.340280 3.249661
## 45 0.9081633 4.345057 3.253641
## 46 0.9265306 4.349224 3.257563
## 47 0.9448980 4.352061 3.259012
## 48 0.9632653 4.354034 3.260748
## 49 0.9816327 4.355639 3.260118
## 50 1.0000000 4.355334 3.260488

# Berdasarkan RMSE
best_loess %>% slice_min(rmse)

##       span     rmse      mae
## 1 0.155102 4.271935 3.206052

# Berdasarkan MAE
best_loess %>% slice_min(mae)

##        span    rmse      mae
## 1 0.1183673 4.27993 3.199382

Span terbaik berdasarkan nilai RMSE adalah span = 0.155102, sedangkan berdasarkan nilai MAE adalah span = 0.1183673. Akan tetapi, selisih nilai RMSE dan MAE keduanya sangat kecil sehingga dilakukan visualisasi terhadap salah satunya saja. Didapatkan hasil visualisasi seperti ini:

library(ggplot2)
ggplot(autodata, aes(horsepower, mpg)) +
  geom_point(alpha = 0.5, color = "#966989") +
  stat_smooth(method = 'loess',
              formula = y~x,
              span = 0.155102, col = "orange", lty=1, se=F, lwd = 0.8)