Data bangkitan untuk peubah X yang digunakan dalam tugas ini terdiri dari 1000 observasi yang berdistribusi normal dengan nilai mean=1 dan standar deviasi=1. Selain itu, dibangkitkan pula data errror yang berdistribusi normal dengan nilai mean=0 dan standar deviasi=1 (default setting). Kemudian, untuk nilai peubah Y didapatkan dari hasil persamaan $y = 5+2x+3x^{2}+error$

set.seed(123)
data.x <- rnorm(1000,1,1)
err <- rnorm(1000)
y <- 5+2*data.x+3*data.x^2+err
datatable(data.frame(y,data.x,err))

Data Triceps

Data yang digunakan untuk ilustrasi berasal dari studi antropometri terhadap 892 perempuan di bawah 50 tahun di tiga desa Gambia di Afrika Barat. Data terdiri dari 3 Kolom yaitu Age, Intriceps dan tricepts. Berikut adalah penjelasannya pada masing-masing kolom:

age : umur responden
Intriceps : logaritma dari ketebalan lipatan kulit triceps
triceps: ketebalan lipatan kulit triceps

Lipatan kulit trisep diperlukan untuk menghitung lingkar otot lengan atas. Ketebalannya memberikan informasi tentang cadangan lemak tubuh, sedangkan massa otot yang dihitung memberikan informasi tentang cadangan protein.

datatable(triceps)

Data Auto

Data Auto diambil dari package ISLR. Package ini berisikan dataset yang digunakan dalam buku “An Introduction to Statistical Learning with Applications in R”. Data Auto yang digunakan dalam tugas ini terdiri dari 392 observasi dan 3 peubah yakni

mpg : miles per gallon
Horsepower : engine horsepower
origin : origin of car (1. American, 2. European, 3. Japanese)

AutoData = Auto %>% select(mpg,horsepower,origin)
mpg <- AutoData$mpg # Peubah Y
horsepower <- AutoData$horsepower # Peubah X
origin <- AutoData$origin
datatable(AutoData)

Data Boston

Boston Dataset berasal dari library MASS. Boston dataset merupakan data atau informasi yang dikumpulkan oleh US Census Service mengenai perumahan di area Boston MASS. Dataset ini terdiri dari 506 observasi dan 14 peubah. Namun dalam tugas ini, hanya akan digunakan 2 peubah, yakni

medv : Median value of owner-occupied homes in $1000’s
lstat : % lower status of the population

data("Boston",package = "MASS")
medv <- Boston$medv
lstat <- Boston$lstat
datatable(data.frame(medv,lstat))

Eksplorasi Data

Plot

plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70), xlab="x",main="Plot Antara Y dan X")

Dapat dilihat dari plot eksplorasi antara peubah x dan peubah y di atas tidak membentuk pola linear melainkan membentuk pola non-linear, yakni polinomial. Maka dari itu, dapat kita asumsikan bahwa regresi polinomial adalah metode yang lebih cocok daripada regresi linear.

Regresi Linear

# Regresi Linear
lin.mod.bangkit <-lm( y~data.x)
summary(lin.mod.bangkit)

## 
## Call:
## lm(formula = y ~ data.x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.686 -2.574 -1.428  1.195 27.185 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.6056     0.1902   24.22   <2e-16 ***
## data.x        8.3790     0.1340   62.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.2 on 998 degrees of freedom
## Multiple R-squared:  0.7967, Adjusted R-squared:  0.7965 
## F-statistic:  3911 on 1 and 998 DF,  p-value: < 2.2e-16

Dari summary model di atas dapat dilihat bahwa ukuran kebaikan model dengan menggunakan R-Square dan Adjusted R-Square. R-Square yang bernilai 0.7967 mengindikasikan bahwa peubah X secara simultan memiliki pengaruh sebesar 79.67 persen terhadap peubah Y. Sedangkan sisanya sebesar 30.33 persen dipengaruhi oleh peubah lain yang tidak digunakan dalam permodelan ini. Kemudian, nilai Adjusted R-Square sebesar 0.7965 mengindikasikan bahwa kemampuan peubah X dalam memengaruhi peubah Y sebesar 79.65. Sedangkan sisanya sebesar 30.35 persen dijelaskan oleh peubah lain yang tidak digunakan dalam permodelan ini.

# Plot
plot(data.x,y,xlim=c( -2,5), ylim=c( -10,70),xlab="X",main="Plot Antara Y dan X",
      title(sub ="Regresi Linear",adj=1))
lines(data.x,lin.mod.bangkit$fitted.values,col="red",lwd=2)

Dapat dilihat dari plot di atas, pola data yang terbentuk kurang dapat mengikuti garis regresi linear (warna merah).

Regresi Polynomial

# Regresi Polynomial
poly.mod.bangkit <- lm( y~data.x+I(data.x^2))
summary(poly.mod.bangkit)

## 
## Call:
## lm(formula = y ~ data.x + I(data.x^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.0319 -0.6942  0.0049  0.7116  3.2855 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.95193    0.04568  108.41   <2e-16 ***
## data.x       2.10732    0.05861   35.95   <2e-16 ***
## I(data.x^2)  2.99081    0.02338  127.93   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.007 on 997 degrees of freedom
## Multiple R-squared:  0.9883, Adjusted R-squared:  0.9883 
## F-statistic: 4.221e+04 on 2 and 997 DF,  p-value: < 2.2e-16

Berdasarkan sumarry model regresi polinomial dengan ordo=2 di atas memiliki nilai kebaikan model yang sudah sangat baik. Nilai R-Square 0.9883 (hampir mendekati 1) memiliki arti bahwa peubah X secara simultan memiliki pengaruh terhadap peubah Y sebesar 98.83 persen. Sedangkan, nilai Adjusted R-Square 0.9883 mengindikasikan bahwa kemampuan peubah X dalam memengaruhi peubah Y dalam permodelan ini sebesar 98.83 persen. Nilai kebaikan model dari model regresi polinomial ini lebih baik daripada model regresi linear.

# Sorting
ix <- sort( data.x,index.return=T)$ix

# Plot
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70),xlab="X",main="Plot Antara Y dan X",
      title(sub ="Regresi Polinomial",adj=1))
lines(data.x[ix], poly.mod.bangkit$fitted.values[ix],col="red", lwd=2)

Selain dapat dilihat dari summary model, secara eksploratif model regresi polinomial juga lebih cocok untuk diterapkan pada data X dan Y karena polanya yang sesuai.

Regresi Fungsi Tangga

Sebelum membentuk model, terlebih dahulu dibentuk breaks dengan memperhatikan nilai range data dari peubah X

range(data.x)

## [1] -1.809775  4.241040

Secara manual, dibentuk 3 breaks

c1 <- as.factor(ifelse(data.x<=0,1,0))
c2 <- as.factor(ifelse(data.x<=2 & data.x>0,1,0))
c3 <- as.factor(ifelse(data.x>2,1,0))

# Regresi Fungsi Tangga
stairs.mod.bangkit <- lm(y~c1+c2+c3)
summary(stairs.mod.bangkit)

## 
## Call:
## lm(formula = y ~ c1 + c2 + c3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.395  -3.534  -0.530   2.527  36.876 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  30.4499     0.4087   74.50   <2e-16 ***
## c11         -25.2395     0.5710  -44.21   <2e-16 ***
## c21         -19.4184     0.4536  -42.81   <2e-16 ***
## c31               NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.121 on 997 degrees of freedom
## Multiple R-squared:  0.698,  Adjusted R-squared:  0.6974 
## F-statistic:  1152 on 2 and 997 DF,  p-value: < 2.2e-16

Dapat dilihat dari summary model regresi fungsi tangga di atas, nilai R-Square dan Adjusted R-Square yang lebih rendah dari pada dua model sebelumnya (regresi linear dan regresi polinomial)

#Plot
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70),xlab="X",main="Plot Antara Y dan X",
      title(sub ="Regresi Fungsi Tangga",adj=1))
lines(data.x[ix], stairs.mod.bangkit$fitted.values[ix],col="red",lwd=2)

Komparasi Model

Plot

Secara eksploratif dibandingkan dengan menggunakan plot dengan garis regresi

plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70))
lines(data.x,lin.mod.bangkit$fitted.values,col="red",lwd=2)
lines(data.x[ix], poly.mod.bangkit$fitted.values[ix],col="blue", lwd=2)
lines(data.x[ix], stairs.mod.bangkit$fitted.values[ix],col="green",lwd=2)
legend("bottomright",legend=c("Regresi Linear","Regresi Polynomial","Regresi Fungsi Tangga"),
       col=c("red","blue","green"),lty=1,cex=0.7,lwd=2)

Dapat dilihat dari plot di atas, pola garis regresi yang paling sesuai adalah pola garis regresi polinomial (warna biru) dibandingkan dengan regresi yang lain. Hal ini mengindikasikan bahwa secara eksploratif, metode regresi polinomial adalah model regresi yang paling cocok diterapkan pada data di atas.

Nilai MSE

Kemudian dilihat pula perbandingan nilai MSE (Mean Square Error). MSE adalah salah satu kriteria yang baik dalam menentukan estimator yang ada.

# Fungsi MSE
MSE = function(pred,actual){
  mean((pred-actual)^2)
}

# Perbandingan
MSEBangkit <- rbind(MSE(predict(lin.mod.bangkit),y),
              MSE(predict(poly.mod.bangkit),y),
              MSE(predict(stairs.mod.bangkit),y))
ModelBangkit <- c("Regresi Linear","Regresi Polynomial","Regresi Fungsi Tangga")
datatable(data.frame("Model"=ModelBangkit,"MSE"=MSEBangkit))

Dari tabel di atas, terlihat bahwa model regresi polinomial memiliki nilai MSE yang paling rendah dengan nilai 1.01065. Hal ini mengindikasikan bahwa model regresi tersebut merupakan model terbaik.

Nilai AIC

AIC <- data.frame("Model Regresi"=c("Linear","Polynomial","Fungsi Tangga"),
                  "AIC"=c(AIC(lin.mod.bangkit),AIC(poly.mod.bangkit),AIC(stairs.mod.bangkit)))
datatable(AIC)

Selain dengan nilai MSE, dilihat pula perbandingan nilai AIC dari model. Semakin rendah nilai AIC, maka semakin baik model tersebut. Dari tabel di atas, model polinomial memiliki nilai AIC yang paling rendah yakni dengan nilai 2856.47001. Hal ini mengindikasikan bahwa model regresi polinomial merupakan model yang terbaik.

Piecewise Cubic Polynomial

dt.all <- cbind(y,data.x)

##knots=1
dt1 <- dt.all[data.x <=1,]
dim(dt1)

## [1] 495   2

Terlihat dimensi data yang kita gunakan, terdapat 495 data yang memiliki kriteria kurang dari sama dengan 1.

dt2 <- dt.all[data.x >1,]
dim(dt2)

## [1] 505   2

Sedangkan, terdapat 505 data yang memiliki kriteria lebih dari 1.

Kemudian, dibentuk fungsi untuk data dt2 dan dt2 yang akan diterapkan pada plot antara Y dan X dengan menggunakan spline berbasis cubic polynomial.

#Plot
plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70), pch=16, col="orange",main="Plot Antara Y dan X",
     xlab="x",title(sub="Piecewise Cubic Polynomial",adj=1))

cub.mod1 <- lm(dt1[,1]~dt1[,2]+I(dt1[,2]^2)+I(dt1[,2]^3))
ix <- sort(dt1[,2], index.return=T)$ix
lines(dt1[ix,2],cub.mod1$fitted.values[ix],col="blue", lwd=2)

cub.mod2 <- lm(dt2[,1]~dt2[,2]+I(dt2[,2]^2)+I(dt2[,2]^3))
ix <- sort(dt2[,2], index.return=T)$ix
lines(dt2[ix,2],cub.mod2$fitted.values[ix],col="blue", lwd=2)

Terlihat dari plot antara Y dan X di atas telah diterapkan lines yang berbasis cubic polynomial. Apabila dicermati, terdapat cut off pada nilai 1. Hal ini karena sebelumnya kita melakukan splitting data menjadi dua bagian dengan kriteria kurang dari sama dengan 1 dan lebih dari 1. Untuk menyambungkan antara fungsi lines dt1 dan dt2 digunakan Truncated Power Basis.

Truncated Power Basis

# Plot
plot(data.x,y,xlim=c( -2,5), ylim=c( -10,70), pch=16,col="orange",main="Plot Antara Y dan X",
     xlab="x",title(sub="Truncated Power Basis",adj=1))
abline(v=1,col="red", lty=2)

hx <- ifelse( data.x>1,(data.x-1)^3,0)

cubspline.mod <- lm( y~data.x+I(data.x^2)+I(data.x^3)+hx)
ix <- sort( data.x,index.return=T)$ix
lines( data.x[ix], cubspline.mod$fitted.values[ix],col="blue", lwd=2)

Terlihat pada plot, antara lines dt1 dan dt2 sudah tersambung. Penyambungan ini kita lakukan dengan membentuk hx.

K-Fold CV

plot( data.x,y,xlim=c(-2,5), ylim=c( -10,70), pch=16,col="orange",main="Plot Antara Y dan X",
     xlab="x",title(sub="K-Fold CV",adj=1))
abline(v=0,col="red", lty=2)
abline(v=2,col="red", lty=2)

hx1 <- ifelse( data.x>0,(data.x-0)^3,0)
hx2 <- ifelse( data.x>2,(data.x-2)^3,0)

cubspline.mod2 <- lm( y~data.x+I(data.x^2)+I(data.x^3)+hx1+hx2)
ix <- sort( data.x,index.return=T)$ix
lines( data.x[ix],cubspline.mod2$fitted.values[ix],col="blue", lwd=2)

Selain itu, kita juga bisa membuat dua knot (pemisah) seperti pada plot di atas. Pada plot di atas, dibentuk tiga fungsi dengan dua knot yang memiliki titik potong di 0 dan 2.

Perbandingan Knot

# 1 Knot dengan 5 fold CV
data.all <- cbind( y,data.x,hx)
set.seed(456)
ind <- sample(1:5,length( data.x),replace=T)
res <- c()
for( i in 1:5){
  dt.train <- data.all[ ind!=i,]
  x.test <- data.all[ ind==i, -1]
  y.test <- data.all[ ind==i,1]
  mod1 <- lm( dt.train[,1]~dt.train[,2]+I( dt.train[,2]^2)+I( dt.train[,2]^3)+dt.train[,3])
  x.test.olah <- cbind(1,x.test[,1], x.test[,1]^2,x.test[,1]^3,x.test[,2])
  beta <- coefficients(mod1)
  prediksi <- x.test.olah%*%beta
  res <- c( res,mean(( y.test-prediksi)^2))
}

# 2 Knot 5 fold CV
data.all2 <- cbind(y,data.x,hx1,hx2)
set.seed(456)
ind2 <- sample(1:5,length( data.x),replace=T)
res2 <- c()
for( i in 1:5){
  dt.train2 <- data.all2[ind2!=i,]
  x.test2 <- data.all2[ind2==i,-1]
  y.test2 <- data.all2[ind2==i,1]
  mod2 <- lm(dt.train2[,1]~dt.train2[,2]+I(dt.train2[,2]^2)+I(dt.train2[,2]^3)+dt.train2[,3]+dt.train2[,4])
  x.test.olah2 <- cbind(1,x.test2[,1],x.test2[,1]^2,x.test2[,1]^3,x.test2[,2],x.test2[,3])
  beta2 <- coefficients(mod2)
  prediksi2 <- x.test.olah2%*%beta2
  res2 <- c(res2,mean((y.test2-prediksi2)^2))
}

datatable(data.frame("Knot"=c(1,2),"Mean"=c(mean(res),mean(res2))))

Terlihat bahwa knot 1 memiliki nilai rata-rata residual lebih kecil daripada knot 2. Hal ini mengindikasikan bahwa knot 1 lebih baik dari pada knot 2.

Smoothing Spline

ss1 <- smooth.spline(data.x,y,all.knots = T)
ss1

## Call:
## smooth.spline(x = data.x, y = y, all.knots = T)
## 
## Smoothing Parameter  spar= 1.499938  lambda= 0.0002745903 (26 iterations)
## Equivalent Degrees of Freedom (Df): 14.84845
## Penalized Criterion (RSS): 994.1601
## GCV: 176.167

Terlihat dari output di atas, nilai spar sebesar 1.4999 dan lambda sebesar 0.00027 dengan iterasi sebanyak 26 kali. Didapatkan juga nilai derajat bebas sebesar 14.85.

plot(data.x,y,xlim=c(-2,5), ylim=c(-10,70), pch=16,col="orange", main="Plot Antara Y dan X",
     xlab="x",title(sub="Smoothing Spline",adj=1))
lines(ss1,col="blue", lwd=2)

Eksplorasi Data

Korelasi

chart.Correlation(triceps)
title(main = "Plot Korelasi Data Triceps")

Terlihat dari plot korelasi di atas, peubah age dan triceps memiliki nilai korelasi 0.58. Peubah yang memiliki nilai korelasi paling tinggi atau memiliki hubungan linear yang kuat adalah peubah intriceps dan triceps dengan nilai korelasi sebesar 0.96.

Plot

# Plot Triceps dan Age
plot(triceps$age,triceps$triceps, main="Plot Antara Triceps dan Age", xlab="Age",ylab = "Triceps")

Terlihat dari plot hubungan antara peubah Age (x) dan peubah Triceps (y) di atas memiliki pola hubungan yang non-linear.

# Plot Triceps dan lntriceps
plot(triceps$lntriceps,triceps$triceps, main="Plot Antara Triceps dan Lntriceps", xlab="Lntriceps",ylab = "Triceps")

Terlihat dari plot hubungan antara peubah Lntriceps (x) dan peubah Triceps (y) di atas memiliki pola hubungan yang cenderung linear dari kiri bawah ke kanan atas.

Regresi Linear

# Regresi Linear
lin.mod.triceps <- lm(triceps~age,data=triceps)
summary(lin.mod.triceps)

## 
## Call:
## lm(formula = triceps ~ age, data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.9512  -2.3965  -0.5154   1.5822  25.1233 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.19717    0.21244   29.17   <2e-16 ***
## age          0.21584    0.01014   21.28   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.007 on 890 degrees of freedom
## Multiple R-squared:  0.3372, Adjusted R-squared:  0.3365 
## F-statistic: 452.8 on 1 and 890 DF,  p-value: < 2.2e-16

Dari summary model regresi linear di atas, terlihat bahwa Intercept dan peubah Age memiliki pengaruh yang signifikan (***) terhadap model. Kita juga dapat melihat ukuran kebaikan model berupa R-Square dan Adjusted R-Suare yag hampir sama. Nilai Square sebesar 0.3372 mengindikasikan bahwa peubah X secara simultan memiliki pengaruh sebesar 33.72 persen terhadap peubah Y. Sedangkan sisanya (1-0.3372) 66.28 persen dipengaruhi oleh peubah lain yang tidak digunakan dalam model. Sedangkan nilai Adjusted R-Square sebesar 0.3365 memiliki arti bahwa kemampuan peubah X memengaruhi peubah Y sebesar 33.65 persen. Sedangkan sisanya, 66.35 persen dijelaskan oleh peubah lain dalam permodelan. Nilai kebaikan model ini cenderung rendah, sehingga dapat dikatakan bahwa model regresi linear tidak cukup baik.

plot(triceps$age,triceps$triceps, main="Plot Antara Triceps dan Age", xlab="Age",ylab = "Triceps",
     title(sub="Regresi Linear",adj=1))
lines(triceps$age,lin.mod.triceps$fitted.values,col="red",lwd=2)

Dapat dilihat secara eksploratif, bahwa garis regresi linear (warna merah) kurang dapat menunjukkan pola dari data.

Regresi Polynomial (Ordo 2)

# Regresi Polynomial Ordo 2
poly2.mod.triceps <- lm(triceps~poly(age,2,raw=T),data=triceps)
summary(poly2.mod.triceps)

## 
## Call:
## lm(formula = triceps ~ poly(age, 2, raw = T), data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.5677  -2.4401  -0.4587   1.6368  24.9961 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.0229191  0.3063806  19.658  < 2e-16 ***
## poly(age, 2, raw = T)1  0.2434733  0.0364403   6.681 4.17e-11 ***
## poly(age, 2, raw = T)2 -0.0006257  0.0007926  -0.789     0.43    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.008 on 889 degrees of freedom
## Multiple R-squared:  0.3377, Adjusted R-squared:  0.3362 
## F-statistic: 226.6 on 2 and 889 DF,  p-value: < 2.2e-16

Berdasarkan summary model regresi polinomial ordo=2 di atas, memiliki nilai R-Squared dan Adjusted R-Squared yang tidak terlalu jauh berbeda, yakni berturun-turut 0.3377 dan 0.3362. Nilai ini menunjukkan bahwa model yang dibentuk tidak cukup baik karena nilai kebaikan modelnya yang cukup rendah.

#Plot
ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Polynomial (Ordo 2)")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,2,raw=T), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Terlihat dari plot di atas, secara eksploratif lines yang terbentuk dari model regresi polinomial ordo=2 kurang dapat menunjukkan pola dari data.

Regresi Polynomial (Ordo 3)

# Regresi Polynomial Ordo 3
poly3.mod.triceps <- lm(triceps~poly(age,3,raw=T),data=triceps)
summary(poly3.mod.triceps)

## 
## Call:
## lm(formula = triceps ~ poly(age, 3, raw = T), data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.5832  -1.9284  -0.5415   1.3283  24.4440 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             8.004e+00  3.831e-01  20.893  < 2e-16 ***
## poly(age, 3, raw = T)1 -3.157e-01  7.721e-02  -4.089 4.73e-05 ***
## poly(age, 3, raw = T)2  3.101e-02  3.964e-03   7.824 1.45e-14 ***
## poly(age, 3, raw = T)3 -4.566e-04  5.612e-05  -8.135 1.38e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.868 on 888 degrees of freedom
## Multiple R-squared:  0.3836, Adjusted R-squared:  0.3815 
## F-statistic: 184.2 on 3 and 888 DF,  p-value: < 2.2e-16

Dari summary di atas terlihat bahwa parameter-parameter signifikan terhadap model. Selain itu, terlihat pula nilai R-Squared dan Adjusted R-Squared yang tidak terlalu tinggi yakni berurut-turut hanya 0.3836 dan 0.3815.

# Plot
ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Polynomial (Ordo 3)")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~poly(x,3,raw=T), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Terlihat dari plot di atas, lines yang terbentuk dari model regresi polinomial ordo=3 lebih dapat mengikuti/menunjukkan pola dari data dibandingkan dengan lines model regresi polinomial ordo=2. Secara eksploratif, model regresi polinomial ordo=3 lebih baik daripada model regresi polinomial ordo=2.

Regresi Fungsi Tangga (Breaks=5)

stairs5.mod.triceps <- lm(triceps~cut(age,5),data=triceps)
summary(stairs5.mod.triceps)

## 
## Call:
## lm(formula = triceps ~ cut(age, 5), data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.5474  -2.0318  -0.4465   1.3682  23.3759 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              7.2318     0.1994  36.260  < 2e-16 ***
## cut(age, 5)(10.6,20.9]   1.6294     0.3244   5.023 6.16e-07 ***
## cut(age, 5)(20.9,31.2]   5.9923     0.4222  14.192  < 2e-16 ***
## cut(age, 5)(31.2,41.5]   7.5156     0.4506  16.678  < 2e-16 ***
## cut(age, 5)(41.5,51.8]   7.4561     0.5543  13.452  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.939 on 887 degrees of freedom
## Multiple R-squared:  0.3617, Adjusted R-squared:  0.3588 
## F-statistic: 125.7 on 4 and 887 DF,  p-value: < 2.2e-16

ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Fungsi Tangga (Breaks=5))")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,5), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Regresi Fungsi Tangga (Breaks=7)

# Regresi Fungsi Tangga Breaks=7
stairs7.mod.triceps <- lm(triceps~cut(age,7),data=triceps)
summary(stairs7.mod.triceps)

## 
## Call:
## lm(formula = triceps ~ cut(age, 7), data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.8063  -1.7592  -0.4366   1.2894  23.1461 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              7.5592     0.2219  34.060  < 2e-16 ***
## cut(age, 7)(7.62,15]    -0.6486     0.3326  -1.950   0.0515 .  
## cut(age, 7)(15,22.3]     3.4534     0.4239   8.146 1.27e-15 ***
## cut(age, 7)(22.3,29.7]   5.8947     0.4604  12.804  < 2e-16 ***
## cut(age, 7)(29.7,37]     7.8471     0.5249  14.949  < 2e-16 ***
## cut(age, 7)(37,44.4]     6.9191     0.5391  12.835  < 2e-16 ***
## cut(age, 7)(44.4,51.8]   6.3013     0.6560   9.606  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.805 on 885 degrees of freedom
## Multiple R-squared:  0.4055, Adjusted R-squared:  0.4014 
## F-statistic: 100.6 on 6 and 885 DF,  p-value: < 2.2e-16

Terlihat dari summary di atas, model regresi fungsi tangga dengan breaks=7 memiliki nilai kebaikan model yang paling tinggi dibandingkan model yang lain walaupun nilainya masih cukup rendah, yakni sekitar 0.4055 untuk R-Squared dan 0.4014 untuk Adjusted R-Squared. Terlihat pula signifikansi dari tiap-tiap potongan mulai dari potongan pertama hingga ketujuh.

ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Fungsi Tangga (Breaks=7))")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~cut(x,7), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Berdasarkan plot, terlihat bahwa lines yang terbentuk lebih dapat menunjukkan pola dari data yang digunakan apabila dibandingkan dengan garis model regresi yang lain.

Komparasi Model

Nilai MSE

# Fungsi MSE
MSE = function(pred,actual){
  mean((pred-actual)^2)
}

# Perbandingan
MSETriceps <- rbind(MSE(predict(lin.mod.triceps),triceps$triceps),
              MSE(predict(poly2.mod.triceps),triceps$triceps),
              MSE(predict(poly3.mod.triceps),triceps$triceps),
              MSE(predict(stairs5.mod.triceps),triceps$triceps),
              MSE(predict(stairs7.mod.triceps),triceps$triceps))
ModelTriceps <- c("Regresi Linear","Regresi Polynomial (Ordo=2)", "Regresi Polynomial (Ordo=3)",
                "Regresi Fungsi Tangga (Breaks=5)", "Regresi Fungsi Tangga (Breaks=7)")
datatable(data.frame("Model"=ModelTriceps,"MSE"=MSETriceps))

Terlihat dari tabel di atas, model regresi fungsi tangga dengan breaks=7 memiliki nilai MSE yang paling kecil.

Nilai AIC

AICTriceps <- data.frame("Model Regresi"=c("Linear","Polynomial (Ordo=2)","Polynomial (Ordo=3)",
                                           "Fungsi Tangga (Breaks=5)","Fungsi Tangga (Breaks=7)"),
                  "AIC"=c(AIC(lin.mod.triceps),AIC(poly2.mod.triceps),AIC(poly3.mod.triceps),
                          AIC(stairs5.mod.triceps),AIC(stairs7.mod.triceps)))
datatable(AICTriceps)

Sama halnya dengan perbandingan nilai MSE, pada tabel di atas model regresi fungsi tangga dengan breaks=7 memiliki nilai AIC yang paling kecil sehingga dapat dikatakan bahwa model tersebut merupakan model yang terbaik.

Cross Validation Evaluation

Regresi Linear

set.seed(123)
cross_val <- vfold_cv(triceps,v=5,strata = "triceps")
metric_linear <- map_dfr(cross_val$splits,
    function(x){
    mod <- lm(triceps ~ age,data=triceps[x$in_id,])
    pred <- predict(mod,newdata=triceps[-x$in_id,])
    truth <- triceps[-x$in_id,]$triceps
    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_linear)

# menghitung rata-rata untuk 10 folds
mean_metric_linear <- colMeans(metric_linear)
mean_metric_linear

##     rmse      mae 
## 3.995795 2.835123

Didapatkan nilai rata-rata RMSE dan MAE dari regresi linear 5 folds yakni berturut-turut sebesar 3.996 dan 2.835.

Regresi Polynomial (Ordo=2)

set.seed(123)
cross_val <- vfold_cv(triceps,v=5,strata = "triceps")
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(triceps ~ poly(age,2,raw = T),data=triceps[x$in_id,])
  pred <- predict(mod,newdata=triceps[-x$in_id,])
  truth <- triceps[-x$in_id,]$triceps
  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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly2 <- colMeans(metric_poly2)
mean_metric_poly2

##     rmse      mae 
## 4.001146 2.855938

Dengan menggunakan Cross Validation Evaluation didapatkan rata-rata nilai dari RMSE dan MAE dari regresi polinomial ordo=2 5 folds berturut-turut sebesar 4.001 dan 2.856.

Regresi Polynomial (Ordo=3)

set.seed(123)
cross_val <- vfold_cv(triceps,v=5,strata = "triceps")

metric_poly3 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(triceps ~ poly(age,3,raw = T),data=triceps[x$in_id,])
  pred <- predict(mod,newdata=triceps[-x$in_id,])
  truth <- triceps[-x$in_id,]$triceps
  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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly3 <- colMeans(metric_poly3)
mean_metric_poly3

##     rmse      mae 
## 3.867164 2.635003

Dengan menggunakan Cross Validation Evaluation didapatkan rata-rata nilai dari RMSE dan MAE regresi polinomial ordo=3 5 folds berturut-turut sebesar 3.867 dan 2.635.

Regresi Fungsi Tangga

set.seed(123)
cross_val <- vfold_cv(triceps,v=5,strata = "triceps")
breaks <- 3:10
best_tangga <- map_dfr(breaks, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
    function(x){
        training <- triceps[x$in_id,]
        training$age <- cut(training$age,i)
        mod <- lm(triceps ~ age,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 <- triceps[-x$in_id,]
        age_new <- cut(testing$age,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
        pred <- predict(mod,newdata=list(age=age_new))
        truth <- testing$triceps
        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)
      }
    )
  datatable(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
datatable(best_tangga)

Terlihat dari tabel di atas, nilai breaks yang memiliki nilai MAE terkecil adalah breaks=9.

Komparasi Model

nilai_metricTriceps <- rbind(mean_metric_linear,
                      mean_metric_poly2,
                      mean_metric_poly3,
                      best_tangga %>% select(-1) %>% slice_min(mae))

nama_modelTriceps <- c("Linear","Polynomial Ordo=2","Polynomial Ordo=3","Fungsi Tangga (Breaks=9)")
datatable(data.frame("Model"=nama_modelTriceps,nilai_metricTriceps))

Terlihat dari tabel di atas, regresi fungsi tangga dengan breaks=9 memiliki nilai RMSE dan MAE terendah di antara model-model yang lain. Hal ini mengindikasikan bahwa model tersebut merupakan model yang terbaik.

Regresi Spline

# Menentukan banyaknya fungsi basis dan knots
dim(bs(triceps$age, df=6))

## [1] 892   6

terlihat dari dimensi peubah Age.

#nilai knots yang ditentukan oleh komputer
attr(bs(triceps$age, df=6),"knots")

##     25%     50%     75% 
##  5.5475 12.2100 24.7275

Didapatkan nilai knots=3.

B-Spline

knot_value_manual_3 = c(10, 20,40)
spline1.mod = lm(triceps ~ bs(age, knots =knot_value_manual_3 ),data=triceps)
summary(spline1.mod)

## 
## Call:
## lm(formula = triceps ~ bs(age, knots = knot_value_manual_3), 
##     data = triceps)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.385  -1.682  -0.393   1.165  22.855 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            8.22533    0.61873  13.294  < 2e-16 ***
## bs(age, knots = knot_value_manual_3)1 -0.06551    1.14972  -0.057    0.955    
## bs(age, knots = knot_value_manual_3)2 -4.30051    0.72301  -5.948 3.90e-09 ***
## bs(age, knots = knot_value_manual_3)3  7.80435    1.17793   6.625 6.00e-11 ***
## bs(age, knots = knot_value_manual_3)4  6.14890    1.27439   4.825 1.65e-06 ***
## bs(age, knots = knot_value_manual_3)5  5.56640    1.42225   3.914 9.78e-05 ***
## bs(age, knots = knot_value_manual_3)6  7.90178    1.54514   5.114 3.87e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.76 on 885 degrees of freedom
## Multiple R-squared:  0.4195, Adjusted R-squared:  0.4155 
## F-statistic: 106.6 on 6 and 885 DF,  p-value: < 2.2e-16

Terlihat summary dari model spline di atas.

# Plot
ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi B-Spline")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~bs(x, knots = knot_value_manual_3), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

# Knot dari komputer
knot_value_pc_df_6 = attr(bs(triceps$age, df=6),"knots")
mod_spline_1 = lm(triceps ~ bs(age, knots =knot_value_pc_df_6 ),data=triceps)
summary(mod_spline_1)

## 
## Call:
## lm(formula = triceps ~ bs(age, knots = knot_value_pc_df_6), data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.0056  -1.7556  -0.2944   1.2008  23.0695 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                            6.5925     0.8770   7.517 1.37e-13 ***
## bs(age, knots = knot_value_pc_df_6)1   3.7961     1.5015   2.528  0.01164 *  
## bs(age, knots = knot_value_pc_df_6)2  -2.0749     0.8884  -2.335  0.01974 *  
## bs(age, knots = knot_value_pc_df_6)3   1.5139     1.1645   1.300  0.19391    
## bs(age, knots = knot_value_pc_df_6)4  11.6394     1.3144   8.855  < 2e-16 ***
## bs(age, knots = knot_value_pc_df_6)5   5.9680     1.5602   3.825  0.00014 ***
## bs(age, knots = knot_value_pc_df_6)6   8.9127     1.4053   6.342 3.60e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.757 on 885 degrees of freedom
## Multiple R-squared:  0.4206, Adjusted R-squared:  0.4167 
## F-statistic: 107.1 on 6 and 885 DF,  p-value: < 2.2e-16

# Plot
ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi B-Spline (Knot dari komputer)")+
  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,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Natural Spline

knot_value_manual_3 = c(10, 20,40)
mod_spline3ns = lm(triceps ~ ns(age, knots = knot_value_manual_3),data=triceps)
summary(mod_spline3ns)

## 
## Call:
## lm(formula = triceps ~ ns(age, knots = knot_value_manual_3), 
##     data = triceps)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.7220  -1.7640  -0.3985   1.1908  22.9684 
## 
## Coefficients:
##                                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                             8.8748     0.3742  23.714  < 2e-16 ***
## ns(age, knots = knot_value_manual_3)1   7.0119     0.6728  10.422  < 2e-16 ***
## ns(age, knots = knot_value_manual_3)2   6.0762     0.8625   7.045 3.72e-12 ***
## ns(age, knots = knot_value_manual_3)3   2.0780     1.0632   1.954    0.051 .  
## ns(age, knots = knot_value_manual_3)4   8.8616     0.9930   8.924  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.762 on 887 degrees of freedom
## Multiple R-squared:  0.4176, Adjusted R-squared:  0.415 
## F-statistic:   159 on 4 and 887 DF,  p-value: < 2.2e-16

# Plot
ggplot(triceps,aes(x=age, y=triceps)) +
        xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Natural Spline")+
  geom_point(alpha=0.55, color="black") +
  stat_smooth(method = "lm", 
               formula = y~ns(x, knots = knot_value_manual_3), 
               lty = 1, col = "red",se = T,lwd=2)+
        theme_update(plot.title = element_text(hjust = 0.5))

Smoothing Splin

smoothspline.mod <- with(data = triceps,smooth.spline(age,triceps))
smoothspline.mod

## Call:
## smooth.spline(x = age, y = triceps)
## 
## Smoothing Parameter  spar= 0.5162704  lambda= 1.232259e-06 (13 iterations)
## Equivalent Degrees of Freedom (Df): 50.00639
## Penalized Criterion (RSS): 8591.581
## GCV: 13.77835

# Plot
pred_data <- broom::augment(smoothspline.mod)

ggplot(pred_data,aes(x=x,y=y))+
  geom_point(alpha=0.55, color="black")+
  geom_line(aes(y=.fitted),col="red",
            lwd=2)+
  xlab("Age")+
  ylab("Triceps")+
   ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Regresi Natural Spline")+
  theme_update(plot.title = element_text(hjust = 0.5))

# Melihat Pengaruh Parameter Lambda
smoothspline.mod.lambda <- data.frame(lambda=seq(0,5,by=0.5)) %>% 
  group_by(lambda) %>% 
  do(broom::augment(with(data = triceps,smooth.spline(age,triceps,lambda = .$lambda))))

p <- ggplot(smoothspline.mod.lambda,
       aes(x=x,y=y))+
  geom_line(aes(y=.fitted),
            col="red",
            lwd=1
            )+
  facet_wrap(~lambda)
p

# df=7
smoothspline7.mod <- with(data = triceps,smooth.spline(age,triceps,df=7))
smoothspline7.mod

## Call:
## smooth.spline(x = age, y = triceps, df = 7)
## 
## Smoothing Parameter  spar= 1.049874  lambda= 0.008827582 (11 iterations)
## Equivalent Degrees of Freedom (Df): 7.000747
## Penalized Criterion (RSS): 10112.16
## GCV: 14.20355

# Plot
ggplot(triceps, aes(age,triceps)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.75,
              col="red",
              lwd=2,
              se=F)+
  xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "Smoothing Spline (df=7")

LOESS

loess.mod <- loess(triceps~age,
                     data = triceps)
summary(loess.mod)

## Call:
## loess(formula = triceps ~ age, data = triceps)
## 
## Number of Observations: 892 
## Equivalent Number of Parameters: 4.6 
## Residual Standard Error: 3.777 
## Trace of smoother matrix: 5.01  (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
model_loess_span <- data.frame(span=seq(0.1,5,by=0.5)) %>% 
  group_by(span) %>% 
  do(broom::augment(loess(triceps~age,
                     data = triceps,span=.$span)))

ggplot(model_loess_span,
       aes(x=age,y=triceps))+
  geom_line(aes(y=.fitted),
            col="red",
            lty=1
            )+
  facet_wrap(~span)

# Plot
ggplot(triceps, aes(age,triceps)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
               formula=y~x,
              span = 0.75,
              col="red",
              lwd=2,
              se=F) +
  xlab("Age")+
        ylab("Triceps")+
        ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "LOESS (Span=0.75)")

Tuning Span

# Tuning Parameter Span
set.seed(123)
cross_val <- vfold_cv(triceps,v=5,strata = "triceps")

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(triceps ~ age,span = i,
         data=triceps[x$in_id,])
pred <- predict(mod,
               newdata=triceps[-x$in_id,])
truth <- triceps[-x$in_id,]$triceps

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
datatable(best_loess)

# Plot
ggplot(triceps, aes(age,triceps)) +
  geom_point(alpha=0.5,color="black") +
  stat_smooth(method='loess',
              formula=y~x,
              span = 0.412244897959184  ,
              col="red",
              lwd=2,
              se=F) +
  xlab("Age")+
  ylab("Triceps")+
   ggtitle("Plot Antara Triceps dan Age")+
        labs(caption = "LOESS (Tuning Span)")+
  theme_update(plot.title = element_text(hjust = 0.5))

Eksplorasi Data

Korelasi

# Perfomance Analytics
chart.Correlation(AutoData)
title(main="Plot Korelasi AutoData")

Berdasarkan plot korelasi di atas, peubah Y (MPG) dan peubah X (Horsepower) memiliki korelasi negatif yang cukup besar yakni sekitar -0.78. Sedangkan antara peubah Y dengan peubah Origin memiliki korelasi positif sekitar 0.57.

Plot

# Plot MPG dan Horsepower
plot(horsepower,mpg, main = "Plot Antara MPG dan Horsepower",ylab="MPG",xlab="Horsepower")

# Plot MPG dan Origin
plot(origin,mpg, main = "Plot Antara MPG dan Origin",ylab="MPG",xlab="Origin")

Terlihat dari plot di atas, pada data tidak terlihat membentuk garis lurus melainkan membentuk seperti lengkungan dari kiri atas ke kanan bawah.

Regresi Linear

# Regresi Linear
lin.mod.auto <- lm(mpg~horsepower)
summary(lin.mod.auto)

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

# Plot
plot(horsepower,mpg, main = "Plot Antara MPG dan Horsepower",ylab="MPG",xlab="Horsepower",
     title(sub ="Regresi Linear",adj=1))
lines(horsepower,lin.mod.auto$fitted.values,col="red",lwd=2)

Regresi Polynomial

# Regresi Polynomial
poly.mod.auto <- lm(mpg~horsepower+I(horsepower^2))
summary(poly.mod.auto)

## 
## Call:
## lm(formula = mpg ~ horsepower + I(horsepower^2))
## 
## 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

# Sorting
ix<-sort(horsepower,index.return=T)$ix

# Plot
plot(horsepower,mpg, main = "Plot Antara MPG dan Horsepower",ylab="MPG",xlab="Horsepower",
     title(sub ="Regresi Polynomial",adj=1))
lines(horsepower[ix],poly.mod.auto$fitted.values[ix],col="red",lwd=2)

Regresi Fungsi Tangga

# Formatting
range(horsepower)

## [1]  46 230

c1 <- as.factor(ifelse(horsepower<=100,1,0))
c2 <- as.factor(ifelse(horsepower<=200,1,0))
c3 <- as.factor(ifelse(horsepower>200,1,0))
datatable(data.frame(mpg,c1,c2,c3))

# Regresi Fungsi Tangga
stairs.mod.auto <- lm(mpg~c1+c2+c3)
summary(stairs.mod.auto)

## 
## Call:
## lm(formula = mpg ~ c1 + c2 + c3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.5917  -3.7167  -0.5917   3.4083  19.0083 
## 
## Coefficients: (1 not defined because of singularities)
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  12.9000     1.8132   7.114 5.44e-12 ***
## c11          10.5589     0.6089  17.342  < 2e-16 ***
## c21           4.1329     1.8769   2.202   0.0283 *  
## c31               NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.734 on 389 degrees of freedom
## Multiple R-squared:  0.4631, Adjusted R-squared:  0.4603 
## F-statistic: 167.7 on 2 and 389 DF,  p-value: < 2.2e-16

# Plot
plot(horsepower,mpg, main = "Plot Antara MPG dan Horsepower",ylab="MPG",xlab="Horsepower",
     title(sub ="Regresi Fungsi Tangga",adj=1))
lines(horsepower[ix],stairs.mod.auto$fitted.values[ix],col="red",lwd=2)

Komparasi Model

Plot

plot(horsepower,mpg, main = "Plot Antara MPG dan Horsepower",ylab="MPG",xlab="Horsepower",
     title(sub ="Perbandingan Model Regresi",adj=1))
lines(horsepower,lin.mod.auto$fitted.values,col="red",lwd=2)
lines(horsepower[ix],poly.mod.auto$fitted.values[ix],col="blue",lwd=2)
lines(horsepower[ix],stairs.mod.auto$fitted.values[ix],col="green",lwd=2)
legend("topright",legend=c("Regresi Linear","Regresi Polynomial","Regresi Fungsi Tangga"),
       col=c("red","blue","green"),lty=1,cex=0.7,lwd=2)

Nilai MSE

# Fungsi MSE
MSE = function(pred,actual){
  mean((pred-actual)^2)
}

# Perbandingan
MSEAuto <- rbind(MSE(predict(lin.mod.auto),mpg),
              MSE(predict(poly.mod.auto),mpg),
              MSE(predict(stairs.mod.auto),mpg))
ModelAuto <- c("Regresi Linear","Regresi Polynomial","Regresi Fungsi Tangga")
datatable(data.frame("Model"=ModelAuto,"MSE"=MSEAuto))

Nilai AIC

AICAuto <- data.frame("Model Regresi"=c("Linear","Polynomial","Fungsi Tangga"),
                  "AIC"=c(AIC(lin.mod.auto),AIC(poly.mod.auto),AIC(stairs.mod.auto)))
datatable(AICAuto)

AIC terendah dimiliki oleh model regresi polynomial dengan nilai sebesar 2274.354. Sehingga dapat disimpulkan jika model regresi polynomial merupakan model terbaik.

Cross Validation Evaluation

Regresi Linear

set.seed(123)
cross_val <- vfold_cv(AutoData,v=5,strata = "mpg")
metric_linear <- map_dfr(cross_val$splits,
    function(x){
    mod <- lm(mpg ~ horsepower,data=triceps[x$in_id,])
    pred <- predict(mod,newdata=AutoData[-x$in_id,])
    truth <- AutoData[-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_linear)

# menghitung rata-rata untuk 10 folds
mean_metric_linear <- colMeans(metric_linear)
mean_metric_linear

##     rmse      mae 
## 4.883215 3.828614

Dari output cross validation di atas terlihat bahwa model regresi linear memiliki rata-rata nilai RMSE dan MAE 5 folds sebesar 4.883215 dan 3.828614.

Regresi Polynomial (Ordo=2)

set.seed(123)
cross_val <- vfold_cv(AutoData,v=5,strata = "mpg")
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(mpg ~ poly(horsepower,2,raw = T),data=AutoData[x$in_id,])
  pred <- predict(mod,newdata=AutoData[-x$in_id,])
  truth <- AutoData[-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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly2 <- colMeans(metric_poly2)
mean_metric_poly2

##     rmse      mae 
## 4.370631 3.284748

Dari output cross validation di atas terlihat bahwa model regresi polinomial ordo=2 memiliki rata-rata nilai RMSE dan MAE 5 folds sebesar 4.370631 dan 3.284748.

Regresi Polynomial (Ordo=3)

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

metric_poly3 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(mpg ~ poly(horsepower,3,raw = T),data=AutoData[x$in_id,])
  pred <- predict(mod,newdata=AutoData[-x$in_id,])
  truth <- AutoData[-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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly3 <- colMeans(metric_poly3)
mean_metric_poly3

##     rmse      mae 
## 4.372623 3.280150

Dari output cross validation di atas terlihat bahwa model regresi polinomial ordo=3 memiliki rata-rata nilai RMSE dan MAE 5 folds sebesar 4.372623 dan 3.280150.

Regresi Fungsi Tangga

set.seed(123)
cross_val <- vfold_cv(AutoData,v=5,strata = "mpg")
breaks <- 3:10
best_tangga <- map_dfr(breaks, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
    function(x){
        training <- AutoData[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 <- AutoData[-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)
      }
    )
  datatable(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
datatable(best_tangga)

Dari output cross validation di atas, terlihat bahwa breaks=7 memiliki nilai MAE yang paling rendah diantara breaks yang lain.

Komparasi Model

nilai_metricAuto <- rbind(mean_metric_linear,
                      mean_metric_poly2,
                      mean_metric_poly3,
                      best_tangga %>% select(-1) %>% slice_min(mae))

nama_modelAuto <- c("Linear","Polynomial Ordo=2","Polynomial Ordo=3","Fungsi Tangga (Breaks=7)")
datatable(data.frame("Model"=nama_modelAuto,nilai_metricAuto))

Berdasarkan tabel di atas, terlihat bahwa model regresi polinomial dengan ordo=3 memiliki nilai MAE yang paling kecil dibandingkan dengan model-model yang lain. Hal ini mengindikasikan bahwa berdasarakan Cross Validation Evalution, model tersebut merupakan model yang terbaik.

Eksplorasi Data

Korelasi

DataBoston <- cbind(medv,lstat)
DataBoston <- data.frame(DataBoston)
chart.Correlation(DataBoston)
title(main = "Plot Korelasi Data Boston")

Plot

plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv")

Regresi Linear

lin.mod.boston <- lm(medv~lstat)
summary(lin.mod.boston)

## 
## Call:
## lm(formula = medv ~ lstat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.168  -3.990  -1.318   2.034  24.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ***
## lstat       -0.95005    0.03873  -24.53   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.216 on 504 degrees of freedom
## Multiple R-squared:  0.5441, Adjusted R-squared:  0.5432 
## F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

# Plot
plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv",
     title(sub ="Regresi Linear",adj=1))
lines(lstat,lin.mod.boston$fitted.values,lwd=2,col="red")

Regresi Polynomial (Ordo=2)

poly2.mod.boston <- lm(medv~poly(lstat,2,raw=T))
summary(poly2.mod.boston)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 2, raw = T))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2834  -3.8313  -0.5295   2.3095  25.4148 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              42.862007   0.872084   49.15   <2e-16 ***
## poly(lstat, 2, raw = T)1 -2.332821   0.123803  -18.84   <2e-16 ***
## poly(lstat, 2, raw = T)2  0.043547   0.003745   11.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.524 on 503 degrees of freedom
## Multiple R-squared:  0.6407, Adjusted R-squared:  0.6393 
## F-statistic: 448.5 on 2 and 503 DF,  p-value: < 2.2e-16

# Plot
plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv",
     title(sub ="Regresi Polynomial (Ordo=2)",adj=1))
lines(lstat,poly2.mod.boston$fitted.values,lwd=2,col="red")

Regresi Polynomial (Ordo=3)

poly3.mod.boston <- lm(medv~poly(lstat,3,raw=T))
summary(poly3.mod.boston)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 3, raw = T))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5441  -3.7122  -0.5145   2.4846  26.4153 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              48.6496253  1.4347240  33.909  < 2e-16 ***
## poly(lstat, 3, raw = T)1 -3.8655928  0.3287861 -11.757  < 2e-16 ***
## poly(lstat, 3, raw = T)2  0.1487385  0.0212987   6.983 9.18e-12 ***
## poly(lstat, 3, raw = T)3 -0.0020039  0.0003997  -5.013 7.43e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.396 on 502 degrees of freedom
## Multiple R-squared:  0.6578, Adjusted R-squared:  0.6558 
## F-statistic: 321.7 on 3 and 502 DF,  p-value: < 2.2e-16

# Plot
plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv",
     title(sub ="Regresi Polynomial (Ordo=3)",adj=1))
lines(lstat,poly3.mod.boston$fitted.values,lwd=2,col="red")

Regresi Fungsi Tangga (Breaks=3)

stairs3.mod.boston <- lm(medv~cut(lstat,3))
summary(stairs3.mod.boston)

## 
## Call:
## lm(formula = medv ~ cut(lstat, 3))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.985  -4.863  -2.072   3.228  23.115 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               26.8847     0.4117   65.30   <2e-16 ***
## cut(lstat, 3)(13.8,25.9] -10.7131     0.7050  -15.20   <2e-16 ***
## cut(lstat, 3)(25.9,38]   -15.0492     1.3715  -10.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.284 on 503 degrees of freedom
## Multiple R-squared:  0.3753, Adjusted R-squared:  0.3728 
## F-statistic: 151.1 on 2 and 503 DF,  p-value: < 2.2e-16

# Plot
plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv",
     title(sub ="Regresi Fungsi Tangga (Breaks=3)",adj=1))
lines(lstat,stairs3.mod.boston$fitted.values,lwd=2,col="red")

Regresi Fungsi Tangga (Breaks=5)

stairs5.mod.boston <- lm(medv~cut(lstat,5))
summary(stairs5.mod.boston)

## 
## Call:
## lm(formula = medv ~ cut(lstat, 5))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.6055  -3.8549  -0.7199   2.4383  29.0656 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               30.5055     0.4730  64.495   <2e-16 ***
## cut(lstat, 5)(8.98,16.2]  -9.5710     0.6689 -14.308   <2e-16 ***
## cut(lstat, 5)(16.2,23.5] -15.1438     0.8119 -18.651   <2e-16 ***
## cut(lstat, 5)(23.5,30.7] -18.6190     1.1533 -16.143   <2e-16 ***
## cut(lstat, 5)(30.7,38]   -18.9166     2.1846  -8.659   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.398 on 501 degrees of freedom
## Multiple R-squared:  0.5198, Adjusted R-squared:  0.516 
## F-statistic: 135.6 on 4 and 501 DF,  p-value: < 2.2e-16

# Plot
plot(lstat,medv,main="Plot Antara Medv dan Lstat",xlab="Lstat",ylab="Medv",
     title(sub ="Regresi Fungsi Tangga (Breaks=3)",adj=1))
lines(lstat,stairs5.mod.boston$fitted.values,lwd=2,col="red")

Komparasi Model

Nilai MSE

# Fungsi MSE
MSE = function(pred,actual){
  mean((pred-actual)^2)
}

# Perbandingan
MSEBoston <- rbind(MSE(predict(lin.mod.boston),medv),
              MSE(predict(poly2.mod.boston),medv),
              MSE(predict(poly3.mod.boston),medv),
              MSE(predict(stairs3.mod.boston),medv),
              MSE(predict(stairs5.mod.boston),medv))
ModelBoston <- c("Regresi Linear","Regresi Polynomial (Ordo=2)","Regresi Polynomial (Ordo=3)",
                 "Regresi Fungsi Tangga (Breaks=3)","Regresi Fungsi Tangga (Breaks=5)")
datatable(data.frame("Model"=ModelBoston,"MSE"=MSEBoston))

Nilai AIC

AICBoston <- data.frame("Model Regresi"=c("Linear","Polynomial (Ordo=2)","Polynomial (Ordo=3)",
                                          "Fungsi Tangga (Breaks=3)","Fungsi Tangga (Breaks=5)"),
                  "AIC"=c(AIC(lin.mod.boston),AIC(poly2.mod.boston),AIC(poly3.mod.boston),
                          AIC(stairs3.mod.boston),AIC(stairs5.mod.boston)))
datatable(AICBoston)

Cross Validation Evaluation

Regresi Linear

set.seed(123)
cross_val <- vfold_cv(DataBoston,v=5,strata = "medv")
metric_linear <- map_dfr(cross_val$splits,
    function(x){
    mod <- lm(medv ~ lstat,data=DataBoston[x$in_id,])
    pred <- predict(mod,newdata=DataBoston[-x$in_id,])
    truth <- DataBoston[-x$in_id,]$medv
    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_linear)

# menghitung rata-rata untuk 10 folds
mean_metric_linear <- colMeans(metric_linear)
mean_metric_linear

##     rmse      mae 
## 6.204479 4.511336

Regresi Polynomial (Ordo=2)

set.seed(123)
cross_val <- vfold_cv(DataBoston,v=5,strata = "medv")
metric_poly2 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(medv ~ poly(lstat,2,raw = T),data=DataBoston[x$in_id,])
  pred <- predict(mod,newdata=DataBoston[-x$in_id,])
  truth <- DataBoston[-x$in_id,]$medv
  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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly2 <- colMeans(metric_poly2)
mean_metric_poly2

##     rmse      mae 
## 5.489874 4.070041

Regresi Polynomial (Ordo=3)

set.seed(123)
cross_val <- vfold_cv(DataBoston,v=5,strata = "medv")

metric_poly3 <- map_dfr(cross_val$splits,
function(x){
  mod <- lm(medv ~ poly(lstat,3,raw = T),data=DataBoston[x$in_id,])
  pred <- predict(mod,newdata=DataBoston[-x$in_id,])
  truth <- DataBoston[-x$in_id,]$medv
  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)

# menghitung rata-rata untuk 5 folds
mean_metric_poly3 <- colMeans(metric_poly3)
mean_metric_poly3

##     rmse      mae 
## 5.356209 4.015634

Regresi Fungsi Tangga

set.seed(123)
cross_val <- vfold_cv(DataBoston,v=5,strata = "medv")
breaks <- 3:10
best_tangga <- map_dfr(breaks, function(i){
    metric_tangga <- map_dfr(cross_val$splits,
    function(x){
        training <- DataBoston[x$in_id,]
        training$lstat <- cut(training$lstat,i)
        mod <- lm(medv ~ lstat,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 <- DataBoston[-x$in_id,]
        lstat_new <- cut(testing$lstat,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
        pred <- predict(mod,newdata=list(lstat=lstat_new))
        truth <- testing$medv
        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)
      }
    )
  datatable(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
datatable(best_tangga)

Komparasi Model

nilai_metricBoston <- rbind(mean_metric_linear,
                      mean_metric_poly2,
                      mean_metric_poly3,
                      best_tangga %>% select(-1) %>% slice_min(mae))

nama_modelBoston <- c("Linear","Polynomial Ordo=2","Polynomial Ordo=3","Fungsi Tangga (Breaks=10)")
datatable(data.frame("Model"=nama_modelBoston,nilai_metricBoston))

Laporan Responsi I (P08 dan P09)

Faisal Arkan

11/3/2022

Menu

Laporan Responsi I (P08 dan P09)

Faisal Arkan

11/3/2022

Menu

Pendahuluan

Library

Cakupan Materi

Dataset

Data Bangkitan

Data Triceps

Data Auto

Data Boston

Data Bangkitan

Eksplorasi Data

Plot

Regresi Linear

Regresi Polynomial

Regresi Fungsi Tangga

Komparasi Model

Plot

Nilai MSE

Nilai AIC

Piecewise Cubic Polynomial

Truncated Power Basis

K-Fold CV

Perbandingan Knot

Smoothing Spline

Data Triceps

Eksplorasi Data

Korelasi

Plot

Regresi Linear

Regresi Polynomial (Ordo 2)

Regresi Polynomial (Ordo 3)

Regresi Fungsi Tangga (Breaks=5)

Regresi Fungsi Tangga (Breaks=7)

Komparasi Model

Nilai MSE

Nilai AIC

Cross Validation Evaluation

Regresi Linear

Regresi Polynomial (Ordo=2)

Regresi Polynomial (Ordo=3)

Regresi Fungsi Tangga

Komparasi Model

Regresi Spline

B-Spline

Natural Spline

Smoothing Splin

LOESS

Tuning Span

Data Auto

Eksplorasi Data

Korelasi

Plot

Regresi Linear

Regresi Polynomial

Regresi Fungsi Tangga

Komparasi Model

Plot

Nilai MSE

Nilai AIC

Cross Validation Evaluation

Regresi Linear

Regresi Polynomial (Ordo=2)

Regresi Polynomial (Ordo=3)

Regresi Fungsi Tangga

Komparasi Model

Data Boston

Eksplorasi Data

Korelasi

Plot

Regresi Linear

Regresi Polynomial (Ordo=2)

Regresi Polynomial (Ordo=3)

Regresi Fungsi Tangga (Breaks=3)

Regresi Fungsi Tangga (Breaks=5)

Komparasi Model

Nilai MSE

Nilai AIC