Tugas Praktikum Individu 1
STA1381 - Pengantar Sains Data

Data yang akan digunakan dalam pemodelan nonlinear ini adalah data spesifikasi kualitas mesin dalam standar horsepower dan tingkat kehematan penggunaan bensin dalam standar miles per galon (mpg) pada 392 kendaraan. Data yang dapat diakses pada software R ini merupakan data yang diambil dan diolah oleh Carnegie Mellon University dan digunakan dalam 1983 American Statistical Association Exposition.

Pra-Processing Data

Memanggil Package di R

Untuk melakukan proses pemodelan nonlinear, diperlukan beberapa library untuk menunjang proses analisis data sebagai berikut.

library("MultiKink")
library(tidyverse)
library(ggplot2)
library(dplyr)
library(purrr)
library(rsample)
library(ISLR)
library(splines)

Pengambilan Data

Proses pengambilan data dapat dilakukan dengan mudah karena sudah tersedia dalam fungsi Auto di software R sebagai berikut,

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

Pemodelan Nonlinear

Metode pemodelan nonlinear yang hendak dicobakan pada gugus data terdiri atas beberapa model regresi dan pemulusan sebagai berikut,

Regresi Linear

model.linear <- lm(mpg~horsepower, data = Auto)
summary(model.linear)

## 
## Call:
## lm(formula = mpg ~ horsepower, data = Auto)
## 
## 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

Mengacu pada keluaran fungsi lm di atas, seluruh parameter regresi linear signifikan pada taraf 5% sehingga dapat diperoleh model regresi linear untuk gugus data terkait sebagai berikut,

\[ y = 39.936 - 0.158x \]

Berdasarkan model regresi linear di atas, dapat dimaknai bahwa perubahan nilai besar horsepower sebesar satu satuan akan menurunkan dugaan rataan rata-rata besar miles per gallon sebesar 0.158 satuan pada tiap kendaraan. Interpretasi ini dapat diamati secara eksploratif melalui grafik sebagai berikut,

ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~x, 
               lty = 1, col = "red",se = T)+
  theme_bw()

Regresi Polinomial Ordo 2

model.poli <- lm(mpg~horsepower+I(horsepower^2), data = Auto)
summary(model.poli)

## 
## Call:
## lm(formula = mpg ~ horsepower + I(horsepower^2), data = Auto)
## 
## 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 ***
## horsepower      -0.4661896  0.0311246  -14.98   <2e-16 ***
## I(horsepower^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

Mengacu pada keluaran fungsi lm di atas, seluruh parameter regresi polinomial dengan ordo 2 signifikan pada taraf 5% sehingga dapat diperoleh model regresi polinomial dengan ordo 2 untuk gugus data terkait sebagai berikut,

\[ y = 56.9 + 0.001x^2-0.466x \]

Berdasarkan model regresi polinomial ordo dua di atas, dapat dimaknai bahwa perubahan nilai besar horsepower sebesar satu satuan akan menurunkan dugaan rataan rata-rata besar miles per gallon sebesar 0.466 satuan dan pada saat yang bersamaan meningkatkan dugaan rataannya secara kuadratik sebesar 0.001 pada tiap kendaraan. Interpretasi ini dapat diamati secara eksploratif melalui grafik sebagai berikut,

ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~x+I(x^2), 
               lty = 1, col = "red",se = T)+
  theme_bw()

Regresi Polinomial Ordo 3

model.poli3 <- lm(mpg~poly(horsepower,3,raw=T), data = Auto)
summary(model.poli3)

## 
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Auto)
## 
## 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

Mengacu pada keluaran fungsi lm di atas, tidak seluruh parameter regresi polinomial dengan ordo 3 signifikan pada taraf 5% sehingga model secara keseluruhan menjadi tidak signifikan untuk menggambarkan pola hubungan antara spesifikasi horsepower dengan tingkat miles per gallon. Meskipun demikian, secara eksploratif, model regresi ini dapat diamati melalui plot sebagai berikut,

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

Regresi Fungsi Tangga

range(Auto$horsepower)

## [1]  46 230

Gugus data yang diamati memiliki range data dari 46 hingga 230 sehingga secara keseluruhan, data dapat dibagi menjadi delapan selang yang masing-masing selangnya terbagi dalam nilai yang bulat.

model.step <- lm(Auto$mpg~cut(horsepower,8),data=Auto)
summary(model.step)

## 
## Call:
## lm(formula = Auto$mpg ~ cut(horsepower, 8), data = Auto)
## 
## 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

Mengacu pada keluaran fungsi lm di atas, seluruh parameter regresi fungsi tangga signifikan pada taraf 5%. Secara eksploratif, pola hubungan data yang digambarkan oleh model regresi ini dapat diamati melalui grafik sebagai berikut,

ggplot(Auto,aes(x=horsepower, y=mpg)) +
       geom_point(alpha=0.55, color="black") +
       stat_smooth(method = "lm", 
               formula = y~cut(x,8), 
               lty = 1, col = "red",se = T)+
  theme_bw()

Berdasarkan fungsi tangga yang berhasil dibuat pada grafik di atas, dapat diamati bahwa terdapat hubungan negatif antara nilai horsepower dengan tingkat miles per gallon. Hal ini dapat dimaknai bahwa peningkatan kualitas kendaraan berdasarkan standar horsepower akan menurunkan dugaan rataan tingkat kehematan penggunaan bensin dalam standar miles per gallon kendaraan.

Regresi Spline

Akan diujikan pemodelan regresi spline dengan derajat bebas 6 dengan pembagian proporsi sebagai berikut,

knot_value_pc_df_6 <- attr(bs(Auto$horsepower, df=6), "knots")
knot_value_pc_df_6

##   25%   50%   75% 
##  75.0  93.5 126.0

model.spline = lm(mpg ~ bs(horsepower,knots =knot_value_pc_df_6),
                  data=Auto)
summary(model.spline)

## 
## Call:
## lm(formula = mpg ~ bs(horsepower, knots = knot_value_pc_df_6), 
##     data = Auto)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5656  -2.5966  -0.2456   2.2574  15.2709 
## 
## Coefficients:
##                                             Estimate Std. Error t value
## (Intercept)                                   32.591      1.931  16.878
## bs(horsepower, knots = knot_value_pc_df_6)1    5.801      3.293   1.761
## bs(horsepower, knots = knot_value_pc_df_6)2   -2.609      1.972  -1.323
## bs(horsepower, knots = knot_value_pc_df_6)3  -11.374      2.241  -5.076
## bs(horsepower, knots = knot_value_pc_df_6)4  -16.977      2.622  -6.475
## bs(horsepower, knots = knot_value_pc_df_6)5  -23.202      3.557  -6.523
## bs(horsepower, knots = knot_value_pc_df_6)6  -18.161      2.829  -6.420
##                                             Pr(>|t|)    
## (Intercept)                                  < 2e-16 ***
## bs(horsepower, knots = knot_value_pc_df_6)1    0.079 .  
## bs(horsepower, knots = knot_value_pc_df_6)2    0.187    
## bs(horsepower, knots = knot_value_pc_df_6)3 6.02e-07 ***
## bs(horsepower, knots = knot_value_pc_df_6)4 2.90e-10 ***
## bs(horsepower, knots = knot_value_pc_df_6)5 2.18e-10 ***
## bs(horsepower, knots = knot_value_pc_df_6)6 4.03e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.294 on 385 degrees of freedom
## Multiple R-squared:  0.702,  Adjusted R-squared:  0.6974 
## F-statistic: 151.2 on 6 and 385 DF,  p-value: < 2.2e-16

Mengacu pada keluaran fungsi lm di atas, tidak seluruh parameter regresi spline signifikan pada taraf 5% sehingga model secara keseluruhan menjadi tidak signifikan untuk menggambarkan pola hubungan antara spesifikasi horsepower dengan tingkat miles per gallon. Meskipun demikian, secara eksploratif, model regresi ini dapat diamati melalui plot sebagai berikut,

ggplot(Auto,aes(x=horsepower, y=mpg)) +
                 geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~bs(x, knots = knot_value_pc_df_6), 
               lty = 1,col="red",se = T)+
  theme_bw()

Pemulusan Spline

model.sms <- with(data=Auto, 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

pred_data <- broom::augment(model.sms)

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

model_sms <- with(data = Auto,smooth.spline(horsepower,mpg,df=7))
model_sms

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

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

Melalui metode pemulusan spline yang teramati pada grafik di atas, dapat diamati bahwa terdapat hubungan yang saling bertolak belakang antara spesifikasi mesin kendaraan dalam standar horsepower dan tingkat kehematan penggunaan bensin dalam standar mpg. Meskipun demikian, dapat diamati pula bahwa grafik semakin melandai sejak nilai horsepower sebesar 175.

Pemulusan LOESS

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

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

Berdasarkan keluaran fungsi loess di atas, digunakan parameter span sebesar 0.75 yang secara eksploratif dapat diamati melalui grafik sebagai berikut,

ggplot(Auto, aes(horsepower,mpg)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.75,
              col="red",
              lty=1,
              se=T)+
  theme_bw()

Melalui metode pemulusan LOESS yang teramati pada grafik di atas, dapat diamati bahwa terdapat hubungan negatif antara spesifikasi mesin kendaraan dalam standar horsepower dan tingkat kehematan penggunaan bensin dalam standar mpg. Meskipun demikian, dapat diamati pula bahwa grafik semakin melandai sejak nilai horsepower sebesar 175 di kisaran angka 15 mpg. Hal ini menunjukkan bahwa dari kendaraan-kendaraan yang teramati, batas bawah tingkat kehematan penggunaan bensinnya adalah 15 mil per galon.

Perbandingan Kebaikan Model

Berdasarkan model-model regresi yang telah dicobakan dan memiliki parameter pemodelan yang seluruhnya signifikan, perlu ditentukan model terbaik yang dapat menggambarkan hubungan antara kedua peubah menggunakan nilai AIC dengan tabulasi sebagai berikut,

aic.lin <- AIC(model.linear)
aic.poli2 <- AIC(model.poli)
aic.tangga <- AIC(model.step)

metode <- c("Regresi Linear", "Regresi Polinomial Ordo 2", 
            "Regresi Fungsi Tangga")
aic <- c(aic.lin, aic.poli2, aic.tangga)
comp <- cbind(Metode=metode, AIC=aic)

kableExtra::kable(comp, caption = "Perbandingan Kebaikan Model Regresi", align = rep("c",2))

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

Perbandingan Kebaikan Model Regresi
Metode	AIC
Regresi Linear	2363.3236578374
Regresi Polinomial Ordo 2	2274.35352235969
Regresi Fungsi Tangga	2289.76066853693

Kesimpulan

Berdasarkan perbandingan kebaikan model regresi, diketahui bahwa model regresi polinomial berordo 2 merupakan model terbaik karena memiliki nilai AIC yang terkecil di antara model-model lain. Model regresi ini memiliki persamaan dan plot sebagai berikut,

\[ y = 56.9 + 0.001x^2-0.466x \]

dengan \(x\) menunjukkan nilai horsepower dan \(y\) menunjukkan tingkat miles per gallon.

ggplot(Auto,aes(x=horsepower, y=mpg)) + 
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~x+I(x^2), 
               lty = 1, col = "red",se = T)+
  theme_bw()

Tugas Praktikum Individu 1 STA1381 - Pengantar Sains Data