1. Soal

Dengan menggunakan data Auto dari package ISLR (pastikan sudah diinstall), gunakan peubah respon mpg dan pilih tiga kolom lainnya untuk dijadikan peubah prediktor.

Apakah ada bukti bahwa peubah-peubah yang dipilih tersebut memiliki pola hubungan non linear. Buktikan dengan visualisasi untuk mendukung jawaban tersebut.
Lakukan pendekatan Smoothing Spline da Local Regression untuk mengenal dengan baik pola hubungan ketiga pasangan peubah yang telah dipilih sebelumnya.

2. Load Packages

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --

## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.1     v forcats 0.5.1

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

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

## In order to avoid name clashes, do not attach 'mlr3measures'. Instead, only load the namespace with `requireNamespace("mlrmeasures")` and access the measures directly via `::`, e.g. `mlr3measures::auc()`.

library(ISLR)
library(grid)

3. Pembahasan

A. Pola Hubungan Pasangan Peubah

Auto <- Auto
data("Auto",package = "ISLR")
datatable(Auto)

Tabel di atas adalah data frame dengan 392 observasi yang terdiri dari 9 peubah. Pada pembahasan ini, dipilih 3 kolom yang akan menjadi peubah prediktor yaitu horsepower, weight, dan displacement. Pola hubungan antara peubah respon mpg dan masing-peubah prediktor akan disajikan secara visual pada penjelasan di bawah ini.

A.1. `mpg` vs `horsepower`

ggplot(Auto,aes(x=horsepower, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("Horse Power VS mile per gallon") +
  xlab("Horse Power")+
  ylab("mile per gallon")+
  theme_grey()

Dari scatter plot di atas, secara umum terlihat bahwa semakin besar kekuatan mesin (horsepower) suatu mobil, semakin pendek jarak yang ditempuh kendaraan itu per galon bahan bakarnya. Artinya semakin kuat mesinnya, maka semakin boros. Dapat dilihat pula bahwa hubungan antara kekuatan mesin kendaraan (pada sumbu X) dan jarak yang dapat ditempuh oleh kendaraan per galon bahan bakarnya (pada sumbu Y) bersifat tidak linear.

A.2. `mpg` vs `weight`

ggplot(Auto,aes(x=weight, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("Weight VS mile per gallon") +
  xlab("Weight")+
  ylab("mile per gallon")+
  theme_grey()

Dari scatter plot di atas, secara umum terlihat bahwa semakin besar berat (weight) suatu mobil, semakin banyak bahan bakar yang dikonsumsi, namun jarak yang ditempuh kendaraan itu per galon bahan bakarnya semakin pendek. Artinya semakin berat mesinnya, maka semakin boros. Dapat dilihat pula bahwa hubungan antara berat mobil (pada sumbu X) dan jarak yang dapat ditempuh oleh kendaraan per galon bahan bakarnya (pada sumbu Y) bersifat tidak linear.

A.3. `mpg` vs `displacement`

ggplot(Auto,aes(x=displacement, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("displacement VS mile per gallon") +
  xlab("displacement")+
  ylab("mile per gallon")+
  theme_grey()

Dari scatter plot di atas, hubungan antara displacement (pada sumbu X) dan jarak yang dapat ditempuh oleh kendaraan per galon bahan bakarnya (pada sumbu Y) bersifat tidak linear.

Jika hubungan masing-masing ketiga peubah prediktor tersebut dengan peubah respon mpg dipaksakan dengan regresi linear sederhana, maka akan terlihat seperti berikut.

plot_hp <- ggplot(Auto,aes(x=horsepower, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("Horse Power VS mile per gallon") +
  xlab("Horse Power")+
  ylab("mile per gallon")+
  stat_smooth(method = "lm", 
               formula = y~x,lty = 1,
               col = "blue",se = F)+
  theme_grey() 

plot_weight <- ggplot(Auto,aes(x=weight, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("Weight VS mile per gallon") +
  xlab("Weight")+
  ylab("mile per gallon")+
  stat_smooth(method = "lm", 
               formula = y~x,lty = 1,
               col = "blue",se = F)+
  theme_grey() 

plot_displacement <- ggplot(Auto,aes(x=displacement, y=mpg)) +
  geom_point(alpha=0.55, color="black") + 
  ggtitle("Displacement VS mile per gallon") +
  xlab("Displacement")+
  ylab("mile per gallon")+
  stat_smooth(method = "lm", 
               formula = y~x,lty = 1,
               col = "blue",se = F)+
  theme_grey() 

plot_grid(plot_hp, plot_weight, plot_displacement)

Tampak bahwa hubungan linear tidak masuk akal karena seolah jika mengikuti garis linearnya maka kendaraan dengan kekuatan mesin yang lebih besar lagi akan menempuh jarak nol atau bahkan negatif per galon bahan bakarnya. Ini jelas tidak mungkin. Demikian halnya untuk hubungan antara peubah mpg vs weight dan mpg vs displacement.

Oleh karena itu, untuk memetakan data dengan baik agar dapat digunakan dalam analisis lanjutan, akan dilakukan pendekatan Smoothing Spline dan Local Regression pada data mpg vs horsepower, mpg vs weight dan mpg vs displacement agar dapat mengenali pola hubungan data yang tidak linear tersebut.

B. Metode Smoothing Spline

Metode Smoothing Spline pada R dilakukan dengan menggunakan fungsi smooth.spline(). Fungsi smooth.spline() tersebut 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. Kalau CV menggunakan MSE sebagai ukuran kebaikan model, sedangkan GCV menggunakan Weighted-MSE. Banyaknya folds (lipatan) pada CV dan GCV ini adalah sejumlah observasinya.

B.1. `mpg` vs `horsepower`

model_sms_hp <- with(data = Auto,smooth.spline(horsepower,mpg))
model_sms_hp

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

Output data di atas menunjukkan bahwa dengan pendekatan Smoothing Spline pada data mpg vs horsepower diperoleh secara otomatis dengan GCV nilai lambda yang sangat kecil yaitu mendekati 0 dengan 12 kali iterasi sebagai parameter lambda terbaik. Dari hasil penentuan model smoothing splines terbaik tersebut, berikut diperlihatkan secara visual hasil Smoothing Spline pada data mpg vs horsepower.

pred_data <- broom::augment(model_sms_hp)

ggplot(pred_data,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("horsepower")+
  ylab("mpg")+
  theme_bw()

Plot data di atas menunjukkan bahwa model smoothing spline yang diperoleh masih terlihat “jumpy” atau masih belum mulus, sehingga berikut akan dilakukan simulasi dengen mengubah parameter lambda agar diproleh model yang lebih mulus dan mampu menggambarkan pola hubungan data antara mpg vs horsepower dengan sangat baik.

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

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

Ilustrasi di atas menunjukkan bahwa semakin besar parameter lambda maka model regresi yang dihasilkan akan semakin mulus, bahkan cenderung linear.

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

model_sms_hp <- with(data = Auto,smooth.spline(horsepower,mpg,df=5))
model_sms_hp

## Call:
## smooth.spline(x = horsepower, y = mpg, df = 5)
## 
## Smoothing Parameter  spar= 0.8891706  lambda= 0.01650666 (14 iterations)
## Equivalent Degrees of Freedom (Df): 5.000611
## Penalized Criterion (RSS): 2564.889
## GCV: 19.02216

pred_data <- broom::augment(model_sms_hp)

ggplot(pred_data,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("horsepower")+
  ylab("mpg")+
  theme_bw()

B.2. `mpg` vs `weight`

model_sms_weight <- with(data = Auto,smooth.spline(weight,mpg))
model_sms_weight

## Call:
## smooth.spline(x = weight, y = mpg)
## 
## Smoothing Parameter  spar= 1.138481  lambda= 0.08299391 (16 iterations)
## Equivalent Degrees of Freedom (Df): 3.674305
## Penalized Criterion (RSS): 6000.624
## GCV: 17.63865

Output data di atas menunjukkan bahwa dengan pendekatan Smoothing Spline pada data mpg vs weight diperoleh secara otomatis dengan GCV nilai lambda terbaik yaitu 0,083 dengan 16 kali iterasi. Dari hasil penentuan model smoothing splines terbaik tersebut, berikut diperlihatkan secara visual hasil Smoothing Spline pada data mpg vs weight.

pred_data <- broom::augment(model_sms_weight)

ggplot(pred_data,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("weight")+
  ylab("mpg")+
  theme_bw()

Plot data di atas menunjukkan bahwa model smoothing spline yang diperoleh sudah cukup mulus untuk dapat mereprentasikan pola hubungan data antara mpg vs weight.

B.3. `mpg` vs `displacement`

model_sms_displacement <- with(data = Auto,smooth.spline(displacement,mpg))
model_sms_displacement

## Call:
## smooth.spline(x = displacement, y = mpg)
## 
## Smoothing Parameter  spar= 0.08978477  lambda= 6.674266e-09 (13 iterations)
## Equivalent Degrees of Freedom (Df): 52.73308
## Penalized Criterion (RSS): 471.1873
## GCV: 16.58384

Output data di atas menunjukkan bahwa dengan pendekatan Smoothing Spline pada data mpg vs displacement diperoleh secara otomatis dengan GCV nilai lambda terbaik yaitu 0 dengan 13 kali iterasi. Dari hasil penentuan model smoothing splines terbaik tersebut, berikut diperlihatkan secara visual hasil Smoothing Spline pada data mpg vs displacement.

pred_data <- broom::augment(model_sms_displacement)

ggplot(pred_data,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("displacement")+
  ylab("mpg")+
  theme_bw()

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

p <- ggplot(model_sms_lambda1,
       aes(x=x,y=y))+
  geom_line(aes(y=.fitted),
            col="blue",
            lty=1
            )+
  facet_wrap(~lambda)
p

Ilustrasi di atas menunjukkan bahwa semakin besar parameter lambda maka model regresi yang dihasilkan akan semakin mulus, bahkan cenderung linear.

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

model_sms_displacement <- with(data = Auto,smooth.spline(displacement,mpg,df=5))
model_sms_displacement

## Call:
## smooth.spline(x = displacement, y = mpg, df = 5)
## 
## Smoothing Parameter  spar= 0.9815859  lambda= 0.01851043 (12 iterations)
## Equivalent Degrees of Freedom (Df): 5.000745
## Penalized Criterion (RSS): 2936.755
## GCV: 19.19855

pred_data1 <- broom::augment(model_sms_displacement)

ggplot(pred_data1,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="blue",
            lty=1)+
  xlab("displacement")+
  ylab("mpg")+
  theme_bw()

C. Metode Local Regression

Masih menggunakan data Auto seperti pada ilustrasi sebelumnya, kali ini kita akan mencoba melakukan pendekatan local regression dengan fungsi loess().

C.1. `mpg` vs `horsepower`

Model_loess_hp <- loess(mpg~horsepower,
                     data = Auto)
summary(Model_loess_hp)

## Call:
## loess(formula = mpg ~ horsepower, data = Auto)
## 
## 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_hp <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~horsepower,
                     data = Auto,span=.$span)))

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

Tuning span dapat dilakukan dengan menggunakan cross-validation

set.seed(123)
cross_val <- vfold_cv(Auto,v=5,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

#berdasarkan rmse
best_loess %>% slice_min(rmse)

#berdasarkan mae
best_loess %>% slice_min(mae)

library(ggplot2)
ggplot(Auto, aes(horsepower,mpg)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.3769231,
              col="blue",
              lty=1,
              se=F)

C.2. `mpg` vs `weight`

Model_loess_weight <- loess(mpg~weight,
                     data = Auto)
summary(Model_loess_weight)

## Call:
## loess(formula = mpg ~ weight, data = Auto)
## 
## Number of Observations: 392 
## Equivalent Number of Parameters: 4.57 
## Residual Standard Error: 4.188 
## Trace of smoother matrix: 4.98  (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_weight <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~weight,
                     data = Auto,span=.$span)))

p2 <- ggplot(model_loess_span_weight,
       aes(x=weight,y=mpg))+
  geom_line(aes(y=.fitted),
            col="blue",
            lty=1
            )+
  facet_wrap(~span)
p2

Tuning span dapat dilakukan dengan menggunakan cross-validation

set.seed(123)
cross_val <- vfold_cv(Auto,v=5,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 ~ weight,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

#berdasarkan rmse
best_loess %>% slice_min(rmse)

#berdasarkan mae
best_loess %>% slice_min(mae)

library(ggplot2)
ggplot(Auto, aes(weight,mpg)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 1,
              col="blue",
              lty=1,
              se=F)

C.3. `mpg` vs `displacement`

Model_loess_displacement <- loess(mpg~displacement,
                     data = Auto)
summary(Model_loess_displacement)

## Call:
## loess(formula = mpg ~ displacement, data = Auto)
## 
## Number of Observations: 392 
## Equivalent Number of Parameters: 4.42 
## Residual Standard Error: 4.372 
## Trace of smoother matrix: 4.82  (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_displacement <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(mpg~displacement,
                     data = Auto,span=.$span)))

p2 <- ggplot(model_loess_span_displacement,
       aes(x=displacement,y=mpg))+
  geom_line(aes(y=.fitted),
            col="blue",
            lty=1
            )+
  facet_wrap(~span)
p2

Tuning span dapat dilakukan dengan menggunakan cross-validation

set.seed(123)
cross_val <- vfold_cv(Auto,v=5,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 ~ displacement,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

#berdasarkan rmse
best_loess %>% slice_min(rmse)

#berdasarkan mae
best_loess %>% slice_min(mae)

library(ggplot2)
ggplot(Auto, aes(displacement,mpg)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.3020408,
              col="blue",
              lty=1,
              se=F)

IPB University-Prodi Statistika dan Sains Data 2021↩︎

Smoothing Spline

Titin Yuniarty1

1. Soal

2. Load Packages

3. Pembahasan

A. Pola Hubungan Pasangan Peubah

A.1. mpg vs horsepower

A.2. mpg vs weight

A.3. mpg vs displacement

B. Metode Smoothing Spline

B.1. mpg vs horsepower

B.2. mpg vs weight

B.3. mpg vs displacement

C. Metode Local Regression

C.1. mpg vs horsepower

C.2. mpg vs weight

C.3. mpg vs displacement

Titin Yuniarty¹

A.1. `mpg` vs `horsepower`

A.2. `mpg` vs `weight`

A.3. `mpg` vs `displacement`

B.1. `mpg` vs `horsepower`

B.2. `mpg` vs `weight`

B.3. `mpg` vs `displacement`

C.1. `mpg` vs `horsepower`

C.2. `mpg` vs `weight`

C.3. `mpg` vs `displacement`