Soal
- Lakukan cross validasi untuk menghasilkan pemodelan mpg vs horsepower optimal menggunakan 3 metode berikut, kemudian bandingkan hasilnya
- Reg polinomial
- Piecewise constant
- Natural cubic splines
- Lakukan hal di atas untuk masing-masing sub-set data berdasarkan asal negara produsennya (Amerika, Eropa, Jepang). Plot model terbaik dalam satu
frame. - Berikan ulasan terhadap hasil anda.
A. Define Data
Berikut ini adalah beberapa table yang akan dipakai untuk penelitian ini diantaranya tabel mpg dan horsepower untuk semua negara, wilayah Amerika, Eropa, dan Jepang.
library(ISLR)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(agricolae) # paket untuk uji lanjut
library(ggpubr) # paket menyusun grafik## Loading required package: ggplot2attach(Auto)## The following object is masked from package:ggplot2:
##
## mpg#Memanggil data yang akan dismulasikan
View(Auto)
Data.Amerika<-Auto %>% filter(origin==1)
Data.Eropa<-Auto %>% filter(origin==2)
Data.Jepang<-Auto %>% filter(origin==3)
Data.komparasi<- Auto%>% mutate("HP/mpg" = horsepower/Auto$mpg)
Data.Anova<-Data.komparasi%>% select ("HP/mpg")
Perlakuan<-Data.komparasi%>% select ("origin")
Data.test.anova<- cbind(Perlakuan, Data.Anova)B. Visualize Data Linier
Membuat pola linier antara variabel horsepower dan mpg serta menghitung nilai rmse nya dengan validasi silang dengan v=5.
library(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.3 v purrr 0.3.4
## 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()
## x car::recode() masks dplyr::recode()
## x purrr::some() masks car::some()library(splines)
library(rsample)
#Membuat plot HP vs MpG
#Cobakan dengan model linier
mod_linear = lm(mpg~horsepower,data=Auto)
summary(mod_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#Plot model linier
plot.linierall<-ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Linier Plot")+
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = T)+
theme_bw()
plot.linierall
#Uji Empiris Regresi Linier
set.seed(123)
cross_valda <- vfold_cv(Auto,v=5,strata = "mpg")
metric_linear1ya <- map_dfr(cross_valda$splits,
function(x){
modnya <- lm(mpg ~ horsepower,
data=Auto[x$in_id,])
preds <- predict(modnya,
newdata=Auto[-x$in_id,])
truths <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truths,
response = preds
)
mae <- mlr3measures::mae(truth = truths,
response = preds
)
metricnya <- c(rmse,mae)
names(metricnya) <- c("rmse","mae")
return(metricnya)
}
)
metric_linear1ya# menghitung rata-rata untuk 5 folds
mean_metric_linear1ya <- colMeans(metric_linear1ya)
mean_metric_linear1ya ## rmse mae
## 4.898592 3.835773Berdasarkan output di atas, terlihat dari plot bahwa secara visual variabel dan mpg dan horsepower tidak memiliki hubungan yang linier. hal ini diperkuat dengan nilai rata-rata mrse nya sebesar 4,89 (setelah ini akan dibandingkan dengan model lainnya)
C. Menentukan Model
Tentukan model non-linear terbaik untuk pasangan peubah mpg dan horsepower (Secara visual dan secara empiris)
1. Regresi Polinomial
Pada pendekatan polinomial ini akan diuji dengan pendekatan polinomial berderajat 2,3, dan 4 untuk dilihat secara visual dan juga secara empiris dengan nilai mrse dan mae nya, namun pada case ini nilai rmse yang akan dijadikan acuan sebagai dasar pengambilan keputusan.
library(tidyverse)
library(splines)
library(rsample)
#Regresi polinomial derajat 2
mod_polinomial2 = lm(mpg ~ poly(horsepower,2,raw = T),
data=Auto)
summary(mod_polinomial2)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = T), 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 ***
## poly(horsepower, 2, raw = T)1 -0.4661896 0.0311246 -14.98 <2e-16 ***
## poly(horsepower, 2, raw = T)2 0.0012305 0.0001221 10.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.374 on 389 degrees of freedom
## Multiple R-squared: 0.6876, Adjusted R-squared: 0.686
## F-statistic: 428 on 2 and 389 DF, p-value: < 2.2e-16#Plot model polinom derajat 2
plot.poli2all<- ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = T)+
labs(title= "Polinom d=2")+
theme_bw()
plot.poli2all 
#Regresi polinomial derajat 3
mod_polinomial3 = lm(mpg ~ poly(horsepower,3,raw = T),
data=Auto)
summary(mod_polinomial3)##
## 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#Plot model polinom derajat 3
plot.poli3all<- 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 = "blue",se = T)+labs(title= "Polinom d=3")+
theme_bw()
plot.poli3all
#Regresi polinomial derajat 4
mod_polinomial4 = lm(mpg ~ poly(horsepower,3,raw = T),
data=Auto)
summary(mod_polinomial4)##
## 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#Plot model polinom derajat 4
plot.poli4all<-ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "blue",se = T)+labs(title= "Polinom d=4")+
theme_bw()
plot.poli4all 
#Perbandingan Model dengan RMSE dan MAE polinom derajat 2
set.seed(123)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_poly_d2 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(horsepower,2,raw = T),
data=Auto[x$in_id,])
pred <- predict(mod,
newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = pred
)
mae <- mlr3measures::mae(truth = truth,
response = pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_poly_d2# menghitung rata-rata untuk 10 folds
mean_metric_poly_d2 <- colMeans(metric_poly_d2)
mean_metric_poly_d2## rmse mae
## 4.351129 3.266460#Perbandingan Model dengan RMSE dan MAE polinom derajat 3
set.seed(123)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_poly_d3 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(horsepower,3,raw = T),
data=Auto[x$in_id,])
pred <- predict(mod,
newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = pred
)
mae <- mlr3measures::mae(truth = truth,
response = pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_poly_d3# menghitung rata-rata untuk 10 folds
mean_metric_poly_d3 <- colMeans(metric_poly_d3)
mean_metric_poly_d3## rmse mae
## 4.355430 3.266945#Perbandingan Model dengan RMSE dan MAE polinom derajat 4
set.seed(333)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
metric_poly_d4 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ poly(horsepower,4,raw = T),
data=Auto[x$in_id,])
pred <- predict(mod,
newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = pred
)
mae <- mlr3measures::mae(truth = truth,
response = pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_poly_d4# menghitung rata-rata untuk 10 folds
mean_metric_poly_d4 <- colMeans(metric_poly_d4)
mean_metric_poly_d4## rmse mae
## 4.368638 3.2901332. Regresi Tangga
Pada pendekatan fungsi tangga ini akan dicari pada knots berapa sehingga memiliki nilai mrse dan mae yang paling kecil dan juga akan ditunjukkan visualisasi datanya.
library(tidyverse)
library(splines)
library(rsample)
#Uji Dengan regresi fungsi tangga
mod_tangga1 = lm(mpg ~ cut(horsepower,5),data=Auto)
summary(mod_tangga1)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 5), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.3033 -3.0220 -0.5413 2.4074 16.6394
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.3033 0.4272 73.27 <2e-16 ***
## cut(horsepower, 5)(82.8,120] -8.2813 0.5642 -14.68 <2e-16 ***
## cut(horsepower, 5)(120,156] -15.2427 0.7210 -21.14 <2e-16 ***
## cut(horsepower, 5)(156,193] -17.5922 1.0036 -17.53 <2e-16 ***
## cut(horsepower, 5)(193,230] -18.5340 1.3767 -13.46 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.719 on 387 degrees of freedom
## Multiple R-squared: 0.6382, Adjusted R-squared: 0.6345
## F-statistic: 170.7 on 4 and 387 DF, p-value: < 2.2e-16#Sebaran plotnya
plot.tanggaall<-ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "blue",se = F)+labs(title= "Fungsi Tangga")+
theme_bw()
plot.tanggaall
mod1 <- lm(mpg ~ horsepower,
data=Auto)
#Jika diuji secara empirin dengan RMSE dan MAE
set.seed(223)
cross_val1 <- vfold_cv(Auto,v=10,strata = "mpg")
az<-cross_val1$splits[[1]]
#Misal ingin dicoba-coba pada interval berapa model terbaiknya, misal dicoba intervalnya 3 -10
breaks <- 3:10
best_tangga1 <- map_dfr(breaks, function(i){
metric_tangga1 <- map_dfr(cross_val1$splits,
function(x){
training1 <- Auto[x$in_id,]
training1$horsepower <- cut(training1$horsepower,i)
mod1 <- lm(mpg ~ horsepower,
data=training1)
labs_x <- levels(mod1$model[,2])
labs_x_breaks <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_x) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_x) ))
testing1 <- Auto[-x$in_id,]
horsepower_new <- cut(testing1$horsepower,c(labs_x_breaks[1,1],labs_x_breaks[,2]))
pred1 <- predict(mod1,
newdata=list(horsepower=horsepower_new))
truth <- testing1$mpg
data_eval1 <- na.omit(data.frame(truth,pred1))
rmse <- mlr3measures::rmse(truth = data_eval1$truth,
response = data_eval1$pred1
)
mae <- mlr3measures::mae(truth = data_eval1$truth,
response = data_eval1$pred1
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_tangga1
# menghitung rata-rata untuk 10 folds
mean_metric_tangga1 <- colMeans(metric_tangga1)
mean_metric_tangga1
})
best_tangga1 <- cbind(breaks=breaks,best_tangga1)
# menampilkan hasil all breaks
best_tangga1#berdasarkan rmse
best_tangga1 %>% slice_min(rmse)#berdasarkan mae
best_tangga1 %>% slice_min(mae)3. Regresi Spline
Pada pendekatan regresi spline ini akan dicoba dengan menggunakan base spline dan juga natural spline untuk kemudian dibandingkan mana yang memiliki nilai rmse yang paling kecil.
library(tidyverse)
library(splines)
library(rsample)
#b-spline
dim(bs(Auto$horsepower, df=10))## [1] 392 10mod_spline3 = lm(mpg ~ bs(horsepower, df=10),data=Auto)
summary(mod_spline3)##
## Call:
## lm(formula = mpg ~ bs(horsepower, df = 10), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.1697 -2.5584 -0.1747 2.3976 14.8299
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 33.9183 2.3274 14.574 < 2e-16 ***
## bs(horsepower, df = 10)1 -0.8287 5.0040 -0.166 0.86855
## bs(horsepower, df = 10)2 3.6912 2.5511 1.447 0.14875
## bs(horsepower, df = 10)3 -7.4827 2.8842 -2.594 0.00984 **
## bs(horsepower, df = 10)4 -6.5317 2.5197 -2.592 0.00990 **
## bs(horsepower, df = 10)5 -12.5247 2.6031 -4.811 2.16e-06 ***
## bs(horsepower, df = 10)6 -13.3021 2.6731 -4.976 9.83e-07 ***
## bs(horsepower, df = 10)7 -14.9004 2.8721 -5.188 3.46e-07 ***
## bs(horsepower, df = 10)8 -22.9147 3.5196 -6.511 2.37e-10 ***
## bs(horsepower, df = 10)9 -20.7134 4.2981 -4.819 2.08e-06 ***
## bs(horsepower, df = 10)10 -20.5287 3.3528 -6.123 2.29e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.288 on 381 degrees of freedom
## Multiple R-squared: 0.7059, Adjusted R-squared: 0.6982
## F-statistic: 91.47 on 10 and 381 DF, p-value: < 2.2e-16#natural spline
mod_spline3ns = lm(mpg ~ ns(horsepower, df=10),data=Auto)
summary(mod_spline3ns)##
## Call:
## lm(formula = mpg ~ ns(horsepower, df = 10), data = Auto)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.5003 -2.4956 -0.2086 2.2577 14.6400
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.462 1.680 19.324 < 2e-16 ***
## ns(horsepower, df = 10)1 -4.600 1.789 -2.572 0.010495 *
## ns(horsepower, df = 10)2 -3.007 2.555 -1.177 0.239996
## ns(horsepower, df = 10)3 -8.663 2.043 -4.239 2.82e-05 ***
## ns(horsepower, df = 10)4 -7.170 2.188 -3.278 0.001142 **
## ns(horsepower, df = 10)5 -14.722 2.125 -6.927 1.84e-11 ***
## ns(horsepower, df = 10)6 -8.446 2.467 -3.424 0.000684 ***
## ns(horsepower, df = 10)7 -16.947 2.080 -8.148 5.38e-15 ***
## ns(horsepower, df = 10)8 -21.715 1.891 -11.483 < 2e-16 ***
## ns(horsepower, df = 10)9 -14.384 4.116 -3.494 0.000531 ***
## ns(horsepower, df = 10)10 -22.326 1.946 -11.476 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.24 on 381 degrees of freedom
## Multiple R-squared: 0.7124, Adjusted R-squared: 0.7049
## F-statistic: 94.39 on 10 and 381 DF, p-value: < 2.2e-16#Buat plot bs
plot.splinebsall<-ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Basic Spline")+
stat_smooth(method = "lm",
formula = y~bs(x, df=10),
lty = 1,se=F)
plot.splinebsall 
#Buat plot ns
plot.splinensall<-ggplot(Auto,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Natural Spline")+
stat_smooth(method = "lm",
formula = y~ns(x, df=10),
lty = 1,se=F)
plot.splinensall
#Secara empiris Bs-Spline
set.seed(123)
cross_val <- vfold_cv(Auto,v=10,strata = "mpg")
df <- 6:12
best_spline3 <- map_dfr(df, function(i){
metric_spline3 <- map_dfr(cross_val$splits,
function(x){
mod <- lm(mpg ~ bs(horsepower,df=i),
data=Auto[x$in_id,])
pred <- predict(mod,
newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = pred
)
mae <- mlr3measures::mae(truth = truth,
response = pred
)
metric <- c(rmse,mae)
names(metric) <- c("rmse","mae")
return(metric)
}
)
metric_spline3
# menghitung rata-rata untuk 10 folds
mean_metric_spline3 <- colMeans(metric_spline3)
mean_metric_spline3
}
)## Warning in bs(horsepower, degree = 3L, knots = c(`25%` = 75, `50%` = 92.5, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`20%` = 72, `40%` = 88, : some
## 'x' values beyond boundary knots may cause ill-conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`16.66667%` = 70, `33.33333%` =
## 84, : some 'x' values beyond boundary knots may cause ill-conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`14.28571%` = 68, `28.57143%`
## = 78.8571428571428, : some 'x' values beyond boundary knots may cause ill-
## conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`12.5%` = 67, `25%` = 75, :
## some 'x' values beyond boundary knots may cause ill-conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`11.11111%` = 67, `22.22222%` =
## 75, : some 'x' values beyond boundary knots may cause ill-conditioned bases## Warning in bs(horsepower, degree = 3L, knots = c(`10%` = 66.3, `20%` = 72, :
## some 'x' values beyond boundary knots may cause ill-conditioned basesbest_spline3 <- cbind(df=df,best_spline3)
# menampilkan hasil all breaks
best_spline3#berdasarkan rmse
best_spline3 %>% slice_min(rmse)#berdasarkan mae
best_spline3 %>% slice_min(mae)#Secara empiris Ns-Spline
set.seed(123)
cross_valns <- vfold_cv(Auto,v=10,strata = "mpg")
#Mencoba dengan jumlah basis 6-12, jumlah basis 6 maka knots = 3 (dengan rumus basis = knots +3)
df <- 6:12
best_spline3ns <- map_dfr(df, function(i){
metric_spline3ns <- map_dfr(cross_valns$splits,
function(x){
modns <- lm(mpg ~ ns(horsepower,df=i),
data=Auto[x$in_id,])
predns <- predict(modns,
newdata=Auto[-x$in_id,])
truth <- Auto[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predns
)
mae <- mlr3measures::mae(truth = truth,
response = predns
)
metricns <- c(rmse,mae)
names(metricns) <- c("rmse","mae")
return(metricns)
}
)
metric_spline3ns
# menghitung rata-rata untuk 10 folds
mean_metric_spline3ns <- colMeans(metric_spline3ns)
mean_metric_spline3ns
}
)
best_spline3ns <- cbind(df=df,best_spline3ns)
# menampilkan hasil all breaks
best_spline3ns#berdasarkan rmse
best_spline3ns %>% slice_min(rmse)#berdasarkan mae
best_spline3ns %>% slice_min(mae)4. Komparasi Model
Setelah didapatkan hasil pendekatan dengan berbagai macam pendekatan, berikut ini adalah hasil komparasi terhadap nilai rmse dan mae nya untuk kemudian ditentukan mana pendekatan yang paling baik untuk mengetahui pola hubungan antara variabel mpg dan horsepower.
library(tidyverse)
library(splines)
library(rsample)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(reshape)##
## Attaching package: 'reshape'## The following objects are masked from 'package:tidyr':
##
## expand, smiths## The following object is masked from 'package:dplyr':
##
## renamelibrary(plyr)## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:reshape':
##
## rename, round_any## The following object is masked from 'package:purrr':
##
## compact## The following object is masked from 'package:ggpubr':
##
## mutate## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize#Perbandingan semua model
nilai_metric <- rbind(mean_metric_linear1ya,
mean_metric_poly_d2,
mean_metric_poly_d3,
mean_metric_poly_d4,
best_tangga1 %>% select(-1) %>% slice_min(mae),
best_spline3 %>% select(-1)%>% slice_min(mae)
)
nama_model <- c("Linear",
"Poly2",
"Poly3","Poly4","Tangga(breaks=8)",
"spline(df=10)"
)
data.frame(nama_model,nilai_metric)#Grafik Gabungan komparasi
grid.arrange(plot.linierall, plot.poli2all, plot.poli3all, plot.poli4all, plot.tanggaall, plot.splinensall, newpage = T, ncol = 2 )
Dari output tersebut dapat disimpulkan bahwa regresi natural cubic spline adalah pendekatan yang paling baik dalam menjelaskan hubungan antara variabel mpg dan horsepower, hal itu terlihat dari nilai rmse nya paling kecil.
D. Amerika Visualize Data Linier
Membuat pola linier antara variabel horsepower dan mpg hanya pada mobil buatan Amerika serta menghitung nilai rmse nya dengan validasi silang.
library(tidyverse)
library(splines)
library(rsample)
#Membuat plot HP vs MpG
#Cobakan dengan model linier
mod_linear1ar = lm(mpg~horsepower,data=Data.Amerika)
summary(mod_linear1ar)##
## Call:
## lm(formula = mpg ~ horsepower, data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7415 -3.1643 -0.6389 2.4849 13.8357
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.476496 0.857484 40.21 <2e-16 ***
## horsepower -0.121320 0.006831 -17.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.257 on 243 degrees of freedom
## Multiple R-squared: 0.5649, Adjusted R-squared: 0.5631
## F-statistic: 315.4 on 1 and 243 DF, p-value: < 2.2e-16#Plot model linier
plot.linierar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = T)+labs(title= "Plot Linier Amerika")+
theme_bw()
plot.linierar
#Uji Empiris Regresi Linier
set.seed(123)
cross_valar <- vfold_cv(Data.Amerika,v=5,strata = "mpg")
metric_linear1ar <- map_dfr(cross_valar$splits,
function(x){
modar <- lm(mpg ~ horsepower,
data=Data.Amerika[x$in_id,])
predar <- predict(modar,
newdata=Data.Amerika[-x$in_id,])
truth <- Data.Amerika[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predar
)
mae <- mlr3measures::mae(truth = truth,
response = predar
)
metricar <- c(rmse,mae)
names(metricar) <- c("rmse","mae")
return(metricar)
}
)
metric_linear1ar# menghitung rata-rata untuk 10 folds
mean_metric_linear1ar <- colMeans(metric_linear1ar)
mean_metric_linear1ar ## rmse mae
## 4.241013 3.326802Berdasarkan output plot di atas, terlihat bahwa hubungan antara variabel mpg dan horsepower cenderung tidak linier terlihat pada nilai X yang berada di ujung atas dan bawah.
E. Amerika Menentukan Model
Tentukan model non-linear terbaik untuk masing pasangan peubah (Secara visual dan/atau secara empiris) untuk mobil buatan amerika
1. Regresi Polinomial Amerika
Pada pendekatan polinomial ini akan diuji dengan pendekatan polinomial berderajat 2,3, dan 4 untuk dilihat secara visual dan juga secara empiris dengan nilai mrse dan mae nya.
library(tidyverse)
library(splines)
library(rsample)
#Regresi polinomial derajat 2
mod_polinomial2ar = lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Amerika)
summary(mod_polinomial2ar)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = T), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9497 -2.5745 -0.0895 2.1583 13.1866
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 51.9150173 2.3893530 21.73 < 2e-16 ***
## poly(horsepower, 2, raw = T)1 -0.4104369 0.0379990 -10.80 < 2e-16 ***
## poly(horsepower, 2, raw = T)2 0.0010776 0.0001398 7.71 3.26e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.822 on 242 degrees of freedom
## Multiple R-squared: 0.6507, Adjusted R-squared: 0.6478
## F-statistic: 225.4 on 2 and 242 DF, p-value: < 2.2e-16#Plot model polinom derajat 2
plot.polinomd2ar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = T)+
labs(title= "Polinom d=2 Amerika")+
theme_bw()
plot.polinomd2ar
#Regresi polinomial derajat 3
mod_polinomial3ar = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Amerika)
summary(mod_polinomial3ar)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.3092 -2.3199 -0.2081 2.0794 13.2928
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.224e+01 7.181e+00 8.667 6.56e-16 ***
## poly(horsepower, 3, raw = T)1 -6.634e-01 1.703e-01 -3.896 0.000127 ***
## poly(horsepower, 3, raw = T)2 2.995e-03 1.266e-03 2.366 0.018796 *
## poly(horsepower, 3, raw = T)3 -4.531e-06 2.974e-06 -1.524 0.128867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.812 on 241 degrees of freedom
## Multiple R-squared: 0.654, Adjusted R-squared: 0.6497
## F-statistic: 151.8 on 3 and 241 DF, p-value: < 2.2e-16#Plot model polinom derajat 3
plot.polinomd3ar<- ggplot(Data.Amerika,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 = "blue",se = T)+labs(title= "Polinom d=3 Amerika")+
theme_bw()
plot.polinomd3ar 
#Regresi polinomial derajat 4
mod_polinomial4ar = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Amerika)
summary(mod_polinomial4ar)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.3092 -2.3199 -0.2081 2.0794 13.2928
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.224e+01 7.181e+00 8.667 6.56e-16 ***
## poly(horsepower, 3, raw = T)1 -6.634e-01 1.703e-01 -3.896 0.000127 ***
## poly(horsepower, 3, raw = T)2 2.995e-03 1.266e-03 2.366 0.018796 *
## poly(horsepower, 3, raw = T)3 -4.531e-06 2.974e-06 -1.524 0.128867
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.812 on 241 degrees of freedom
## Multiple R-squared: 0.654, Adjusted R-squared: 0.6497
## F-statistic: 151.8 on 3 and 241 DF, p-value: < 2.2e-16#Plot model polinom derajat 4
plot.polinomd4ar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "blue",se = T)+labs(title= "Polinom d=4 Amerika")+
theme_bw()
plot.polinomd4ar 
#Perbandingan Model dengan RMSE dan MAE polinom derajat 2
set.seed(123)
cross_valar2 <- vfold_cv(Data.Amerika,v=10,strata = "mpg")
metric_poly_d2ar <- map_dfr(cross_valar2$splits,
function(x){
modar2 <- lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Amerika[x$in_id,])
predar2 <- predict(modar2,
newdata=Data.Amerika[-x$in_id,])
truth <- Data.Amerika[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predar2
)
mae <- mlr3measures::mae(truth = truth,
response = predar2
)
metricar2 <- c(rmse,mae)
names(metricar2) <- c("rmse","mae")
return(metricar2)
}
)
metric_poly_d2ar# menghitung rata-rata untuk 10 folds
mean_metric_poly_d2ar <- colMeans(metric_poly_d2ar)
mean_metric_poly_d2ar## rmse mae
## 3.799433 2.887688#Perbandingan Model dengan RMSE dan MAE polinom derajat 3
set.seed(123)
cross_valar3 <- vfold_cv(Data.Amerika,v=10,strata = "mpg")
metric_poly_d3ar <- map_dfr(cross_valar3$splits,
function(x){
modar3 <- lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Amerika[x$in_id,])
predar3 <- predict(modar3,
newdata=Data.Amerika[-x$in_id,])
truth <- Data.Amerika[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predar3
)
mae <- mlr3measures::mae(truth = truth,
response = predar3
)
metricar3 <- c(rmse,mae)
names(metricar3) <- c("rmse","mae")
return(metricar3)
}
)
metric_poly_d3ar# menghitung rata-rata untuk 10 folds
mean_metric_poly_d3ar <- colMeans(metric_poly_d3ar)
mean_metric_poly_d3ar## rmse mae
## 3.804542 2.901911#Perbandingan Model dengan RMSE dan MAE polinom derajat 4
set.seed(123)
cross_valar4 <- vfold_cv(Data.Amerika,v=10,strata = "mpg")
metric_poly_d4ar <- map_dfr(cross_valar4$splits,
function(x){
modar4 <- lm(mpg ~ poly(horsepower,4,raw = T),
data=Data.Amerika[x$in_id,])
predar4 <- predict(modar4,
newdata=Data.Amerika[-x$in_id,])
truth <- Data.Amerika[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predar4
)
mae <- mlr3measures::mae(truth = truth,
response = predar4
)
metricar4 <- c(rmse,mae)
names(metricar4) <- c("rmse","mae")
return(metricar4)
}
)
metric_poly_d4ar# menghitung rata-rata untuk 10 folds
mean_metric_poly_d4ar <- colMeans(metric_poly_d4ar)
mean_metric_poly_d4ar## rmse mae
## 3.823584 2.9008482. Regresi Tangga Amerika
Pada pendekatan fungsi tangga ini akan dicari pada knots berapa sehingga memiliki nilai mrse dan mae yang paling kecil dan juga akan ditunjukkan visualisasi datanya.
library(tidyverse)
library(splines)
library(rsample)
#Uji Dengan regresi fungsi tangga
mod_tangga1ar = lm(mpg ~ cut(horsepower,9),data=Data.Amerika)
summary(mod_tangga1ar)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 9), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.946 -1.946 -0.375 2.077 13.054
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.5933 0.9521 34.23 < 2e-16 ***
## cut(horsepower, 9)(71.8,91.6] -7.6477 1.0519 -7.27 5.24e-12 ***
## cut(horsepower, 9)(91.6,111] -12.6916 1.0682 -11.88 < 2e-16 ***
## cut(horsepower, 9)(111,131] -13.8267 1.3465 -10.27 < 2e-16 ***
## cut(horsepower, 9)(131,151] -17.1706 1.1025 -15.57 < 2e-16 ***
## cut(horsepower, 9)(151,171] -18.3933 1.2892 -14.27 < 2e-16 ***
## cut(horsepower, 9)(171,190] -18.9010 1.3973 -13.53 < 2e-16 ***
## cut(horsepower, 9)(190,210] -21.2600 1.7813 -11.94 < 2e-16 ***
## cut(horsepower, 9)(210,230] -19.2183 1.6144 -11.90 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.688 on 236 degrees of freedom
## Multiple R-squared: 0.6829, Adjusted R-squared: 0.6722
## F-statistic: 63.53 on 8 and 236 DF, p-value: < 2.2e-16#Sebaran plotnya
plot.tanggaar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,9),
lty = 1, col = "blue",se = F)+labs(title= "Fungsi Tangga Amerika")+
theme_bw()
plot.tanggaar
mod1ar <- lm(mpg ~ horsepower,
data=Data.Amerika)
#Jika diuji secara empirin dengan RMSE dan MAE
set.seed(223)
cross_val1ar <- vfold_cv(Data.Amerika,v=10,strata = "mpg")
azar<-cross_val1ar$splits[[1]]
#Misal ingin dicoba-coba pada interval berapa model terbaiknya, misal dicoba intervalnya 3 -10
breaks <- 3:10
best_tangga1ar <- map_dfr(breaks, function(i){
metric_tangga1ar <- map_dfr(cross_val1ar$splits,
function(x){
training1ar <- Data.Amerika[x$in_id,]
training1ar$horsepower <- cut(training1ar$horsepower,i)
mod1art <- lm(mpg ~ horsepower,
data=training1ar)
labs_xar <- levels(mod1art$model[,2])
labs_x_breaksar <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_xar) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_xar) ))
testing1ar <- Data.Amerika[-x$in_id,]
horsepower_newar <- cut(testing1ar$horsepower,c(labs_x_breaksar[1,1],labs_x_breaksar[,2]))
pred1ar <- predict(mod1art,
newdata=list(horsepower=horsepower_newar))
truth <- testing1ar$mpg
data_eval1ar <- na.omit(data.frame(truth,pred1ar))
rmse <- mlr3measures::rmse(truth = data_eval1ar$truth,
response = data_eval1ar$pred1ar
)
mae <- mlr3measures::mae(truth = data_eval1ar$truth,
response = data_eval1ar$pred1ar
)
metricar <- c(rmse,mae)
names(metricar) <- c("rmse","mae")
return(metricar)
}
)
metric_tangga1ar
# menghitung rata-rata untuk 10 folds
mean_metric_tangga1ar <- colMeans(metric_tangga1ar)
mean_metric_tangga1ar
})
best_tangga1ar <- cbind(breaks=breaks,best_tangga1ar)
# menampilkan hasil all breaks
best_tangga1ar#berdasarkan rmse
best_tangga1ar %>% slice_min(rmse)#berdasarkan mae
best_tangga1ar %>% slice_min(mae)3. Regresi Spline Amerika
Pada pendekatan regresi spline ini akan dicoba dengan menggunakan base spline dan juga natural spline untuk kemudian dibandingkan mana yang memiliki nilai rmse dan mae yang paling kecil.
library(tidyverse)
library(splines)
library(rsample)
#b-spline
dim(bs(Data.Amerika$horsepower,df=19))## [1] 245 19mod_spline3ar = lm(mpg ~ bs(horsepower,df=19),data=Data.Amerika)
summary(mod_spline3ar)##
## Call:
## lm(formula = mpg ~ bs(horsepower, df = 19), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.824 -2.021 -0.076 1.998 11.713
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 27.493 3.639 7.556 1.04e-12 ***
## bs(horsepower, df = 19)1 14.864 6.568 2.263 0.024569 *
## bs(horsepower, df = 19)2 -1.059 5.030 -0.210 0.833526
## bs(horsepower, df = 19)3 -2.646 4.352 -0.608 0.543811
## bs(horsepower, df = 19)4 1.715 3.994 0.429 0.668049
## bs(horsepower, df = 19)5 -6.723 4.072 -1.651 0.100150
## bs(horsepower, df = 19)6 -1.610 3.909 -0.412 0.680826
## bs(horsepower, df = 19)7 -5.876 4.408 -1.333 0.183889
## bs(horsepower, df = 19)8 -11.113 4.197 -2.648 0.008682 **
## bs(horsepower, df = 19)9 -6.074 4.073 -1.491 0.137249
## bs(horsepower, df = 19)10 -8.197 3.929 -2.086 0.038101 *
## bs(horsepower, df = 19)11 -3.313 5.268 -0.629 0.530034
## bs(horsepower, df = 19)12 -16.022 5.039 -3.179 0.001683 **
## bs(horsepower, df = 19)13 -9.499 4.136 -2.297 0.022541 *
## bs(horsepower, df = 19)14 -13.167 3.804 -3.461 0.000643 ***
## bs(horsepower, df = 19)15 -13.333 4.485 -2.972 0.003275 **
## bs(horsepower, df = 19)16 -13.384 4.378 -3.057 0.002504 **
## bs(horsepower, df = 19)17 -14.564 5.377 -2.709 0.007273 **
## bs(horsepower, df = 19)18 -17.734 5.539 -3.202 0.001564 **
## bs(horsepower, df = 19)19 -11.790 4.590 -2.569 0.010856 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.658 on 225 degrees of freedom
## Multiple R-squared: 0.7025, Adjusted R-squared: 0.6774
## F-statistic: 27.97 on 19 and 225 DF, p-value: < 2.2e-16#natural spline
mod_spline3nsar = lm(mpg ~ ns(horsepower, df=19),data=Data.Amerika)
summary(mod_spline3nsar)##
## Call:
## lm(formula = mpg ~ ns(horsepower, df = 19), data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -11.5848 -1.8571 -0.0835 1.8266 13.6121
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.1430 3.1678 8.884 < 2e-16 ***
## ns(horsepower, df = 19)1 -12.5712 3.3260 -3.780 0.000201 ***
## ns(horsepower, df = 19)2 8.9323 4.1235 2.166 0.031347 *
## ns(horsepower, df = 19)3 -8.9752 3.4610 -2.593 0.010131 *
## ns(horsepower, df = 19)4 -3.6048 3.6778 -0.980 0.328065
## ns(horsepower, df = 19)5 -3.2940 3.6873 -0.893 0.372634
## ns(horsepower, df = 19)6 -7.5447 3.9045 -1.932 0.054575 .
## ns(horsepower, df = 19)7 -11.3055 3.6124 -3.130 0.001982 **
## ns(horsepower, df = 19)8 -5.8213 3.6681 -1.587 0.113915
## ns(horsepower, df = 19)9 -10.2674 3.4639 -2.964 0.003362 **
## ns(horsepower, df = 19)10 0.2586 4.6586 0.056 0.955782
## ns(horsepower, df = 19)11 -17.2669 4.3223 -3.995 8.77e-05 ***
## ns(horsepower, df = 19)12 -9.0732 3.9112 -2.320 0.021249 *
## ns(horsepower, df = 19)13 -13.7884 3.7722 -3.655 0.000320 ***
## ns(horsepower, df = 19)14 -13.2823 3.6987 -3.591 0.000404 ***
## ns(horsepower, df = 19)15 -14.9049 4.0756 -3.657 0.000318 ***
## ns(horsepower, df = 19)16 -13.1652 3.8021 -3.463 0.000640 ***
## ns(horsepower, df = 19)17 -21.5057 3.0242 -7.111 1.51e-11 ***
## ns(horsepower, df = 19)18 -6.2460 7.3408 -0.851 0.395750
## ns(horsepower, df = 19)19 -19.1802 2.3140 -8.289 1.04e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.495 on 225 degrees of freedom
## Multiple R-squared: 0.7284, Adjusted R-squared: 0.7055
## F-statistic: 31.76 on 19 and 225 DF, p-value: < 2.2e-16#Buat plot bs
plot.splinebsar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Basic Spline")+
stat_smooth(method = "lm",
formula = y~bs(x, df=19 ),
lty = 1,se=F)
plot.splinebsar 
#Buat plot ns
plot.splinensar<-ggplot(Data.Amerika,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Natural Spline")+
stat_smooth(method = "lm",
formula = y~ns(x, df=19),
lty = 1,se=F)
plot.splinensar
#Secara empiris Ns-Spline
set.seed(123)
cross_valnsar <- vfold_cv(Data.Amerika,v=5,strata = "mpg")
#Mencoba dengan jumlah basis 6-12, jumlah basis 6 maka knots = 3 (dengan rumus basis = knots +3)
df <- 6:20
best_spline3nsar <- map_dfr(df, function(i){
metric_spline3nsar <- map_dfr(cross_valnsar$splits,
function(x){
modnsar <- lm(mpg ~ ns(horsepower,df=i),
data=Data.Amerika[x$in_id,])
prednsar <- predict(modnsar,
newdata=Data.Amerika[-x$in_id,])
truth <- Data.Amerika[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = prednsar
)
mae <- mlr3measures::mae(truth = truth,
response = prednsar
)
metricnsar <- c(rmse,mae)
names(metricnsar) <- c("rmse","mae")
return(metricnsar)
}
)
metric_spline3nsar
# menghitung rata-rata untuk 5 folds
mean_metric_spline3nsar <- colMeans(metric_spline3nsar)
mean_metric_spline3nsar
}
)
best_spline3nsar <- cbind(df=df,best_spline3nsar)
# menampilkan hasil all breaks
best_spline3nsar#berdasarkan rmse
best_spline3nsar %>% slice_min(rmse)#berdasarkan mae
best_spline3nsar %>% slice_min(mae)4. Komparasi Model Amerika
Setelah didapatkan hasil pendekatan dengan berbagai macam pendekatan, berikut ini adalah hasil komparasi terhadap nilai rmse dan mae nya untuk kemudian ditentukan mana pendekatan yang paling baik untuk mengetahui pola hubungan antara variabel mpg dan horsepower.
library(ggplot2)
library(tidyverse)
library(splines)
library(rsample)
library(gridExtra)
#Perbandingan semua model
nilai_metric <- rbind(mean_metric_linear1ar,
mean_metric_poly_d2ar,
mean_metric_poly_d3ar,
mean_metric_poly_d4ar,
best_tangga1ar %>% select(-1) %>% slice_min(rmse),
best_spline3nsar %>% select(-1)%>% slice_min(rmse)
)
nama_model <- c("Linear",
"Poly2",
"Poly3","Poly4","Tangga(breaks=9)",
"Natural Spline(df=19)"
)
data.frame(nama_model,nilai_metric)#Grafik Gabungan komparasi
grid.arrange(plot.linierar, plot.polinomd2ar, plot.polinomd3ar, plot.polinomd4ar, plot.tanggaar, plot.splinensar)
Dari output tersebut dapat disimpulkan bahwa regresi natural cubic spline adalah pendekatan yang paling baik dalam menjelaskan hubungan antara variabel mpg dan horsepower, hal itu terlihat dari nilai rmse nya paling kecil.
F. Eropa Visualize Data Linier
Membuat pola linier antara variabel horsepower dan mpg hanya pada mobil buatan Eropa serta menghitung nilai rmse nya dengan validasi silang.
library(tidyverse)
library(splines)
library(rsample)
#Membuat plot HP vs MpG
#Cobakan dengan model linier
mod_linear1ar = lm(mpg~horsepower,data=Data.Amerika)
summary(mod_linear1ar)##
## Call:
## lm(formula = mpg ~ horsepower, data = Data.Amerika)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.7415 -3.1643 -0.6389 2.4849 13.8357
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.476496 0.857484 40.21 <2e-16 ***
## horsepower -0.121320 0.006831 -17.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.257 on 243 degrees of freedom
## Multiple R-squared: 0.5649, Adjusted R-squared: 0.5631
## F-statistic: 315.4 on 1 and 243 DF, p-value: < 2.2e-16#Plot model linier
plot.linierer<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = T)+labs(title= "Plot Linier Eropa")+
theme_bw()
plot.linierer
#Uji Empiris Regresi Linier
set.seed(123)
cross_valer <- vfold_cv(Data.Eropa,v=10,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_linear1er <- map_dfr(cross_valer$splits,
function(x){
moder <- lm(mpg ~ horsepower,
data=Data.Eropa[x$in_id,])
preder <- predict(moder,
newdata=Data.Eropa[-x$in_id,])
truth <- Data.Eropa[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = preder
)
mae <- mlr3measures::mae(truth = truth,
response = preder
)
metricer <- c(rmse,mae)
names(metricer) <- c("rmse","mae")
return(metricer)
}
)
metric_linear1er# menghitung rata-rata untuk 10 folds
mean_metric_linear1er <- colMeans(metric_linear1er)
mean_metric_linear1er ## rmse mae
## 4.755026 3.865495Berdasarkan output plot di atas, terlihat bahwa hubungan antara variabel mpg dan horsepower cenderung tidak linier terlihat pada nilai X yang berada di ujung atas dan bawah.
G. Eropa Menentukan Model
Tentukan model non-linear terbaik untuk masing pasangan peubah (Secara visual dan/atau secara empiris) untuk mobil buatan Eropa
1. Regresi Polinomial Eropa
Pada pendekatan polinomial ini akan diuji dengan pendekatan polinomial berderajat 2,3, dan 4 untuk dilihat secara visual dan juga secara empiris dengan nilai mrse dan mae nya.
library(tidyverse)
library(splines)
library(rsample)
#Regresi polinomial derajat 2
mod_polinomial2er = lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Eropa)
summary(mod_polinomial2er)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = T), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6282 -2.7724 -0.4052 2.2773 12.9676
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 47.1311209 8.2852287 5.689 3.3e-07 ***
## poly(horsepower, 2, raw = T)1 -0.2631496 0.1994015 -1.320 0.192
## poly(horsepower, 2, raw = T)2 0.0002425 0.0011574 0.210 0.835
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.899 on 65 degrees of freedom
## Multiple R-squared: 0.4622, Adjusted R-squared: 0.4456
## F-statistic: 27.93 on 2 and 65 DF, p-value: 1.76e-09#Plot model polinom derajat 2
plot.polinomd2er<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = T)+
labs(title= "Polinom d=2 Eropa")+
theme_bw()
plot.polinomd2er
#Regresi polinomial derajat 3
mod_polinomial3er = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Eropa)
summary(mod_polinomial3er)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6346 -2.6497 -0.8141 2.1774 12.8935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.060e+01 2.781e+01 1.460 0.149
## poly(horsepower, 3, raw = T)1 -7.245e-03 1.058e+00 -0.007 0.995
## poly(horsepower, 3, raw = T)2 -2.924e-03 1.291e-02 -0.226 0.822
## poly(horsepower, 3, raw = T)3 1.241e-05 5.039e-05 0.246 0.806
##
## Residual standard error: 4.935 on 64 degrees of freedom
## Multiple R-squared: 0.4627, Adjusted R-squared: 0.4375
## F-statistic: 18.37 on 3 and 64 DF, p-value: 1.041e-08#Plot model polinom derajat 3
plot.polinomd3er<- ggplot(Data.Eropa,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 = "blue",se = T)+labs(title= "Polinom d=3 Eropa")+
theme_bw()
plot.polinomd3er 
#Regresi polinomial derajat 4
mod_polinomial4er = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Eropa)
summary(mod_polinomial4er)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.6346 -2.6497 -0.8141 2.1774 12.8935
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.060e+01 2.781e+01 1.460 0.149
## poly(horsepower, 3, raw = T)1 -7.245e-03 1.058e+00 -0.007 0.995
## poly(horsepower, 3, raw = T)2 -2.924e-03 1.291e-02 -0.226 0.822
## poly(horsepower, 3, raw = T)3 1.241e-05 5.039e-05 0.246 0.806
##
## Residual standard error: 4.935 on 64 degrees of freedom
## Multiple R-squared: 0.4627, Adjusted R-squared: 0.4375
## F-statistic: 18.37 on 3 and 64 DF, p-value: 1.041e-08#Plot model polinom derajat 4
plot.polinomd4er<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "blue",se = T)+labs(title= "Polinom d=4 Eropa")+
theme_bw()
plot.polinomd4er 
#Perbandingan Model dengan RMSE dan MAE polinom derajat 2
set.seed(123)
cross_valer2 <- vfold_cv(Data.Eropa,v=10,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d2er <- map_dfr(cross_valer2$splits,
function(x){
moder2 <- lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Eropa[x$in_id,])
preder2 <- predict(moder2,
newdata=Data.Eropa[-x$in_id,])
truth <- Data.Eropa[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = preder2
)
mae <- mlr3measures::mae(truth = truth,
response = preder2
)
metricer2 <- c(rmse,mae)
names(metricer2) <- c("rmse","mae")
return(metricer2)
}
)
metric_poly_d2er# menghitung rata-rata untuk 10 folds
mean_metric_poly_d2er <- colMeans(metric_poly_d2er)
mean_metric_poly_d2er## rmse mae
## 4.805676 3.906728#Perbandingan Model dengan RMSE dan MAE polinom derajat 3
set.seed(123)
cross_valer3 <- vfold_cv(Data.Eropa,v=10,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d3er <- map_dfr(cross_valer3$splits,
function(x){
moder3 <- lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Eropa[x$in_id,])
preder3 <- predict(moder3,
newdata=Data.Eropa[-x$in_id,])
truth <- Data.Eropa[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = preder3
)
mae <- mlr3measures::mae(truth = truth,
response = preder3
)
metricer3 <- c(rmse,mae)
names(metricer3) <- c("rmse","mae")
return(metricer3)
}
)
metric_poly_d3er# menghitung rata-rata untuk 10 folds
mean_metric_poly_d3er <- colMeans(metric_poly_d3er)
mean_metric_poly_d3er## rmse mae
## 4.832309 3.944713#Perbandingan Model dengan RMSE dan MAE polinom derajat 4
set.seed(123)
cross_valer4 <- vfold_cv(Data.Eropa,v=10,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d4er <- map_dfr(cross_valer4$splits,
function(x){
moder4 <- lm(mpg ~ poly(horsepower,4,raw = T),
data=Data.Eropa[x$in_id,])
preder4 <- predict(moder4,
newdata=Data.Eropa[-x$in_id,])
truth <- Data.Eropa[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = preder4
)
mae <- mlr3measures::mae(truth = truth,
response = preder4
)
metricer4 <- c(rmse,mae)
names(metricer4) <- c("rmse","mae")
return(metricer4)
}
)
metric_poly_d4er# menghitung rata-rata untuk 10 folds
mean_metric_poly_d4er <- colMeans(metric_poly_d4er)
mean_metric_poly_d4er## rmse mae
## 4.878846 3.9846342. Regresi Tangga Eropa
Pada pendekatan fungsi tangga ini akan dicari pada knots berapa sehingga memiliki nilai mrse dan mae yang paling kecil dan juga akan ditunjukkan visualisasi datanya.
library(tidyverse)
library(splines)
library(rsample)
#Uji Dengan regresi fungsi tangga
mod_tangga1er = lm(mpg ~ cut(horsepower,3),data=Data.Eropa)
summary(mod_tangga1er)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 3), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.147 -4.180 -1.247 2.868 15.150
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.147 1.006 30.965 < 2e-16 ***
## cut(horsepower, 3)(75,104] -4.797 1.448 -3.313 0.00151 **
## cut(horsepower, 3)(104,133] -10.667 2.012 -5.302 1.47e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.509 on 65 degrees of freedom
## Multiple R-squared: 0.3199, Adjusted R-squared: 0.299
## F-statistic: 15.29 on 2 and 65 DF, p-value: 3.618e-06#Sebaran plotnya
plot.tanggaer<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,3),
lty = 1, col = "blue",se = F)+labs(title= "Fungsi Tangga Eropa")+
theme_bw()
plot.tanggaer
mod1er <- lm(mpg ~ horsepower,
data=Data.Eropa)
View(Data.Eropa)
#Jika diuji secara empirin dengan RMSE dan MAE
set.seed(223)
cross_val1er <- vfold_cv(Data.Eropa,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.azer<-cross_val1er$splits[[1]]
#Misal ingin dicoba-coba pada interval berapa model terbaiknya, misal dicoba intervalnya 3 -10
breaks <- 2:5
best_tangga1er <- map_dfr(breaks, function(i){
metric_tangga1er <- map_dfr(cross_val1er$splits,
function(x){
training1er <- Data.Eropa[x$in_id,]
training1er$horsepower <- cut(training1er$horsepower,i)
mod1ert <- lm(mpg ~ horsepower,
data=training1er)
labs_xer <- levels(mod1ert$model[,2])
labs_x_breakser <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_xer) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_xer) ))
testing1er <- Data.Eropa[-x$in_id,]
horsepower_newer <- cut(testing1er$horsepower,c(labs_x_breakser[1,1],labs_x_breakser[,2]))
pred1er <- predict(mod1ert,
newdata=list(horsepower=horsepower_newer))
truth <- testing1er$mpg
data_eval1er <- na.omit(data.frame(truth,pred1er))
rmse <- mlr3measures::rmse(truth = data_eval1er$truth,
response = data_eval1er$pred1er
)
mae <- mlr3measures::mae(truth = data_eval1er$truth,
response = data_eval1er$pred1er
)
metricer <- c(rmse,mae)
names(metricer) <- c("rmse","mae")
return(metricer)
}
)
metric_tangga1er
# menghitung rata-rata untuk 5 folds
mean_metric_tangga1er <- colMeans(metric_tangga1er)
mean_metric_tangga1er
})
best_tangga1er <- cbind(breaks=breaks,best_tangga1er)
# menampilkan hasil all breaks
best_tangga1er#berdasarkan rmse
best_tangga1er %>% slice_min(rmse)#berdasarkan mae
best_tangga1er %>% slice_min(mae)3. Regresi Spline Eropa
Pada pendekatan regresi spline ini akan dicoba dengan menggunakan base spline dan juga naturan spline untuk kemudian dibandingkan mana yang memiliki nilai rmsep dan mae yang paling kecil.
library(tidyverse)
library(splines)
library(rsample)
#b-spline
dim(bs(Data.Eropa$horsepower,df=5))## [1] 68 5mod_spline3er = lm(mpg ~ bs(horsepower,df=5),data=Data.Eropa)
summary(mod_spline3er)##
## Call:
## lm(formula = mpg ~ bs(horsepower, df = 5), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.8438 -2.8751 -0.5414 2.4144 12.9182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.8809 2.6363 13.231 < 2e-16 ***
## bs(horsepower, df = 5)1 -0.7563 6.6537 -0.114 0.909867
## bs(horsepower, df = 5)2 -3.7651 3.2651 -1.153 0.253286
## bs(horsepower, df = 5)3 -12.7701 5.1452 -2.482 0.015790 *
## bs(horsepower, df = 5)4 -13.6222 5.3519 -2.545 0.013419 *
## bs(horsepower, df = 5)5 -18.9876 5.3309 -3.562 0.000715 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.002 on 62 degrees of freedom
## Multiple R-squared: 0.4653, Adjusted R-squared: 0.4222
## F-statistic: 10.79 on 5 and 62 DF, p-value: 1.719e-07#natural spline
mod_spline3nser = lm(mpg ~ ns(horsepower, df=5),data=Data.Eropa)
summary(mod_spline3nser)##
## Call:
## lm(formula = mpg ~ ns(horsepower, df = 5), data = Data.Eropa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.1062 -2.8126 -0.5225 2.3099 12.3513
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 35.668 2.214 16.114 < 2e-16 ***
## ns(horsepower, df = 5)1 -5.545 2.651 -2.092 0.040567 *
## ns(horsepower, df = 5)2 -9.781 3.388 -2.887 0.005345 **
## ns(horsepower, df = 5)3 -11.202 4.009 -2.794 0.006910 **
## ns(horsepower, df = 5)4 -21.297 6.192 -3.439 0.001048 **
## ns(horsepower, df = 5)5 -15.971 4.095 -3.900 0.000239 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.995 on 62 degrees of freedom
## Multiple R-squared: 0.4667, Adjusted R-squared: 0.4237
## F-statistic: 10.85 on 5 and 62 DF, p-value: 1.591e-07#Buat plot bs
plot.splinebser<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Basic Spline")+
stat_smooth(method = "lm",
formula = y~bs(x, df=5 ),
lty = 1,se=F)
plot.splinebser 
#Buat plot ns
plot.splinenser<-ggplot(Data.Eropa,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Natural Spline")+
stat_smooth(method = "lm",
formula = y~ns(x, df=5),
lty = 1,se=F)
plot.splinenser
#Secara empiris Ns-Spline
set.seed(123)
cross_valnser <- vfold_cv(Data.Eropa,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.#Mencoba dengan jumlah basis 5-8, jumlah basis 6 maka knots = 3 (dengan rumus basis = knots +3)
df <- 5:8
best_spline3nser <- map_dfr(df, function(i){
metric_spline3nser <- map_dfr(cross_valnser$splits,
function(x){
modnser <- lm(mpg ~ ns(horsepower,df=i),
data=Data.Eropa[x$in_id,])
prednser <- predict(modnser,
newdata=Data.Eropa[-x$in_id,])
truth <- Data.Eropa[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = prednser
)
mae <- mlr3measures::mae(truth = truth,
response = prednser
)
metricnser <- c(rmse,mae)
names(metricnser) <- c("rmse","mae")
return(metricnser)
}
)
metric_spline3nser
# menghitung rata-rata untuk 5 folds
mean_metric_spline3nser <- colMeans(metric_spline3nser)
mean_metric_spline3nser
}
)
best_spline3nser <- cbind(df=df,best_spline3nser)
# menampilkan hasil all breaks
best_spline3nser#berdasarkan rmse
best_spline3nser %>% slice_min(rmse)#berdasarkan mae
best_spline3nser %>% slice_min(mae)4. Komparasi Model Eropa
Setelah didapatkan hasil pendekatan dengan berbagai macam pendekatan, berikut ini adalah hasil komparasi terhadap nilai rmse dan mae nya untuk kemudian ditentukan mana pendekatan yang paling baik untuk mengetahui pola hubungan antara variabel mpg dan horsepower.
library(ggplot2)
library(tidyverse)
library(splines)
library(rsample)
library(gridExtra)
#Perbandingan semua model
nilai_metric <- rbind(mean_metric_linear1er,
mean_metric_poly_d2er,
mean_metric_poly_d3er,
mean_metric_poly_d4er,
best_tangga1er %>% select(-1) %>% slice_min(rmse),
best_spline3nser %>% select(-1)%>% slice_min(rmse)
)
nama_model <- c("Linear",
"Poly2",
"Poly3","Poly4","Tangga(breaks=5)",
"Natural Spline(df5)"
)
data.frame(nama_model,nilai_metric)#Grafik Gabungan komparasi
grid.arrange(plot.linierer, plot.polinomd2er, plot.polinomd3er, plot.polinomd4er, plot.tanggaer, plot.splinenser)
Dari output tersebut dapat disimpulkan bahwa regresi polinomial derajat 2 adalah pendekatan yang paling baik dalam menjelaskan hubungan antara variabel mpg dan horsepower, hal itu terlihat dari nilai rmse nya paling kecil.
H. Jepang Visualize Data Linier
Membuat pola linier antara variabel horsepower dan mpg hanya pada mobil buatan Jepang serta menghitung nilai rmse nya dengan validasi silang.
library(tidyverse)
library(splines)
library(rsample)
#Membuat plot HP vs MpG
#Cobakan dengan model linier
mod_linear1jp = lm(mpg~horsepower,data=Data.Jepang)
summary(mod_linear1jp)##
## Call:
## lm(formula = mpg ~ horsepower, data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.112 -2.758 -0.751 2.563 14.249
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 48.8162 2.3555 20.724 < 2e-16 ***
## horsepower -0.2300 0.0288 -7.986 1.08e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.533 on 77 degrees of freedom
## Multiple R-squared: 0.4531, Adjusted R-squared: 0.446
## F-statistic: 63.78 on 1 and 77 DF, p-value: 1.08e-11#Plot model linier
plot.linierjp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~x,lty = 1,
col = "blue",se = T)+labs(title= "Plot Linier Jepang")+
theme_bw()
plot.linierjp
#Uji Empiris Regresi Linier
set.seed(123)
cross_valjp <- vfold_cv(Data.Jepang,v=10,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_linear1jp <- map_dfr(cross_valjp$splits,
function(x){
modjp <- lm(mpg ~ horsepower,
data=Data.Jepang[x$in_id,])
predjp <- predict(modjp,
newdata=Data.Jepang[-x$in_id,])
truth <- Data.Jepang[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predjp
)
mae <- mlr3measures::mae(truth = truth,
response = predjp
)
metricjp <- c(rmse,mae)
names(metricjp) <- c("rmse","mae")
return(metricjp)
}
)
metric_linear1jp# menghitung rata-rata untuk 10 folds
mean_metric_linear1jp <- colMeans(metric_linear1jp)
mean_metric_linear1jp## rmse mae
## 4.454860 3.568692Berdasarkan output plot di atas, terlihat bahwa hubungan antara variabel mpg dan horsepower cenderung tidak linier terlihat pada nilai X yang berada di ujung atas dan bawah.
G. Jepang Menentukan Model
Tentukan model non-linear terbaik untuk masing pasangan peubah (Secara visual dan/atau secara empiris) untuk mobil buatan Jepang
1. Regresi Polinomial Jepang
Pada pendekatan polinomial ini akan diuji dengan pendekatan polinomial berderajat 2,3, dan 4 untuk dilihat secara visual dan juga secara empiris dengan nilai mrse dan mae nya.
library(tidyverse)
library(splines)
library(rsample)
#Regresi polinomial derajat 2
mod_polinomial2jp = lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Jepang)
summary(mod_polinomial2jp)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 2, raw = T), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.1949 -2.5568 -0.6165 2.2408 12.5211
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 68.401527 9.731163 7.029 7.74e-10 ***
## poly(horsepower, 2, raw = T)1 -0.710534 0.233642 -3.041 0.00323 **
## poly(horsepower, 2, raw = T)2 0.002808 0.001355 2.072 0.04169 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.439 on 76 degrees of freedom
## Multiple R-squared: 0.4823, Adjusted R-squared: 0.4687
## F-statistic: 35.4 on 2 and 76 DF, p-value: 1.365e-11#Plot model polinom derajat 2
plot.polinomd2jp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,2,raw=T),
lty = 1, col = "blue",se = T)+
labs(title= "Polinom d=2 Jepang")+
theme_bw()
plot.polinomd2jp
#Regresi polinomial derajat 3
mod_polinomial3jp = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Jepang)
summary(mod_polinomial3jp)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8602 -2.7115 -0.5224 2.1985 11.6985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.163e+01 3.494e+01 -2.050 0.04385 *
## poly(horsepower, 3, raw = T)1 4.308e+00 1.230e+00 3.503 0.00078 ***
## poly(horsepower, 3, raw = T)2 -5.498e-02 1.400e-02 -3.926 0.00019 ***
## poly(horsepower, 3, raw = T)3 2.142e-04 5.170e-05 4.142 8.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.031 on 75 degrees of freedom
## Multiple R-squared: 0.5787, Adjusted R-squared: 0.5618
## F-statistic: 34.34 on 3 and 75 DF, p-value: 4.485e-14#Plot model polinom derajat 3
plot.polinomd3jp<- ggplot(Data.Jepang,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 = "blue",se = T)+labs(title= "Polinom d=3 Jepang")+
theme_bw()
plot.polinomd3jp 
#Regresi polinomial derajat 4
mod_polinomial4jp = lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Jepang)
summary(mod_polinomial4jp)##
## Call:
## lm(formula = mpg ~ poly(horsepower, 3, raw = T), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.8602 -2.7115 -0.5224 2.1985 11.6985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.163e+01 3.494e+01 -2.050 0.04385 *
## poly(horsepower, 3, raw = T)1 4.308e+00 1.230e+00 3.503 0.00078 ***
## poly(horsepower, 3, raw = T)2 -5.498e-02 1.400e-02 -3.926 0.00019 ***
## poly(horsepower, 3, raw = T)3 2.142e-04 5.170e-05 4.142 8.92e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.031 on 75 degrees of freedom
## Multiple R-squared: 0.5787, Adjusted R-squared: 0.5618
## F-statistic: 34.34 on 3 and 75 DF, p-value: 4.485e-14#Plot model polinom derajat 4
plot.polinomd4jp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~poly(x,4,raw=T),
lty = 1, col = "blue",se = T)+labs(title= "Polinom d=4 Jepang")+
theme_bw()
plot.polinomd4jp 
#Perbandingan Model dengan RMSE dan MAE polinom derajat 2
set.seed(123)
cross_valjp2 <- vfold_cv(Data.Jepang,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d2jp <- map_dfr(cross_valjp2$splits,
function(x){
modjp2 <- lm(mpg ~ poly(horsepower,2,raw = T),
data=Data.Jepang[x$in_id,])
predjp2 <- predict(modjp2,
newdata=Data.Jepang[-x$in_id,])
truth <- Data.Jepang[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predjp2
)
mae <- mlr3measures::mae(truth = truth,
response = predjp2
)
metricjp2 <- c(rmse,mae)
names(metricjp2) <- c("rmse","mae")
return(metricjp2)
}
)
metric_poly_d2jp# menghitung rata-rata untuk 5 folds
mean_metric_poly_d2jp <- colMeans(metric_poly_d2jp)
mean_metric_poly_d2jp## rmse mae
## 4.669195 3.580362#Perbandingan Model dengan RMSE dan MAE polinom derajat 3
set.seed(123)
cross_valjp3 <- vfold_cv(Data.Jepang,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d3jp <- map_dfr(cross_valjp3$splits,
function(x){
modjp3 <- lm(mpg ~ poly(horsepower,3,raw = T),
data=Data.Jepang[x$in_id,])
predjp3 <- predict(modjp3,
newdata=Data.Jepang[-x$in_id,])
truth <- Data.Jepang[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predjp3
)
mae <- mlr3measures::mae(truth = truth,
response = predjp3
)
metricjp3 <- c(rmse,mae)
names(metricjp3) <- c("rmse","mae")
return(metricjp3)
}
)
metric_poly_d3jp# menghitung rata-rata untuk 5 folds
mean_metric_poly_d3jp <- colMeans(metric_poly_d3jp)
mean_metric_poly_d3jp## rmse mae
## 4.081574 3.074327#Perbandingan Model dengan RMSE dan MAE polinom derajat 4
set.seed(123)
cross_valjp4 <- vfold_cv(Data.Jepang,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.metric_poly_d4jp <- map_dfr(cross_valjp4$splits,
function(x){
modjp4 <- lm(mpg ~ poly(horsepower,4,raw = T),
data=Data.Jepang[x$in_id,])
predjp4 <- predict(modjp4,
newdata=Data.Jepang[-x$in_id,])
truth <- Data.Jepang[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = predjp4
)
mae <- mlr3measures::mae(truth = truth,
response = predjp4
)
metricjp4 <- c(rmse,mae)
names(metricjp4) <- c("rmse","mae")
return(metricjp4)
}
)
metric_poly_d4jp# menghitung rata-rata untuk 5 folds
mean_metric_poly_d4jp <- colMeans(metric_poly_d4jp)
mean_metric_poly_d4jp## rmse mae
## 4.320975 3.2401852. Regresi Tangga Jepang
Pada pendekatan fungsi tangga ini akan dicari pada knots berapa sehingga memiliki nilai mrse dan mae yang paling kecil dan juga akan ditunjukkan visualisasi datanya.
library(tidyverse)
library(splines)
library(rsample)
#Uji Dengan regresi fungsi tangga
mod_tangga1jp = lm(mpg ~ cut(horsepower,5),data=Data.Jepang)
summary(mod_tangga1jp)##
## Call:
## lm(formula = mpg ~ cut(horsepower, 5), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.7444 -2.8774 -0.7774 2.1891 11.7226
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 34.8774 0.7651 45.585 < 2e-16 ***
## cut(horsepower, 5)(68,84] -2.4441 1.3399 -1.824 0.072169 .
## cut(horsepower, 5)(84,100] -9.1330 1.1214 -8.144 6.91e-12 ***
## cut(horsepower, 5)(100,116] -12.9108 2.5758 -5.012 3.56e-06 ***
## cut(horsepower, 5)(116,132] -9.2441 2.5758 -3.589 0.000594 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.26 on 74 degrees of freedom
## Multiple R-squared: 0.5358, Adjusted R-squared: 0.5107
## F-statistic: 21.35 on 4 and 74 DF, p-value: 9.702e-12#Sebaran plotnya
plot.tanggajp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") +
stat_smooth(method = "lm",
formula = y~cut(x,5),
lty = 1, col = "blue",se = F)+labs(title= "Fungsi Tangga Jepang")+
theme_bw()
plot.tanggajp
mod1jp <- lm(mpg ~ horsepower,
data=Data.Jepang)
#Jika diuji secara empirin dengan RMSE dan MAE
set.seed(223)
cross_val1jp <- vfold_cv(Data.Jepang,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.azjp<-cross_val1jp$splits[[1]]
#Misal ingin dicoba-coba pada interval berapa model terbaiknya, misal dicoba intervalnya 2 -5
breaks <- 2:5
best_tangga1jp <- map_dfr(breaks, function(i){
metric_tangga1jp <- map_dfr(cross_val1jp$splits,
function(x){
training1jp <- Data.Jepang[x$in_id,]
training1jp$horsepower <- cut(training1jp$horsepower,i)
mod1jpt <- lm(mpg ~ horsepower,
data=training1jp)
labs_xjp <- levels(mod1jpt$model[,2])
labs_x_breaksjp <- cbind(lower = as.numeric( sub("\\((.+),.*", "\\1", labs_xjp) ),
upper = as.numeric( sub("[^,]*,([^]]*)\\]", "\\1", labs_xjp) ))
testing1jp <- Data.Jepang[-x$in_id,]
horsepower_newjp <- cut(testing1jp$horsepower,c(labs_x_breaksjp[1,1],labs_x_breaksjp[,2]))
pred1jp <- predict(mod1jpt,
newdata=list(horsepower=horsepower_newjp))
truth <- testing1jp$mpg
data_eval1jp <- na.omit(data.frame(truth,pred1jp))
rmse <- mlr3measures::rmse(truth = data_eval1jp$truth,
response = data_eval1jp$pred1jp
)
mae <- mlr3measures::mae(truth = data_eval1jp$truth,
response = data_eval1jp$pred1jp
)
metricjp <- c(rmse,mae)
names(metricjp) <- c("rmse","mae")
return(metricjp)
}
)
metric_tangga1jp
# menghitung rata-rata untuk 5 folds
mean_metric_tangga1jp <- colMeans(metric_tangga1jp)
mean_metric_tangga1jp
})
best_tangga1jp <- cbind(breaks=breaks,best_tangga1jp)
# menampilkan hasil all breaks
best_tangga1jp#berdasarkan rmse
best_tangga1jp %>% slice_min(rmse)#berdasarkan mae
best_tangga1jp %>% slice_min(mae)3. Regresi Spline Jepang
Pada pendekatan regresi spline ini akan dicoba dengan menggunakan base spline dan juga natural spline untuk kemudian dibandingkan mana yang memiliki nilai rmsep dan mae yang paling kecil.
library(tidyverse)
library(splines)
library(rsample)
#b-spline
dim(bs(Data.Jepang$horsepower,df=7))## [1] 79 7mod_spline3jp = lm(mpg ~ bs(horsepower,df=7),data=Data.Jepang)
summary(mod_spline3jp)##
## Call:
## lm(formula = mpg ~ bs(horsepower, df = 7), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8.7291 -2.6979 -0.2466 1.4271 10.8597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 31.9499 2.3461 13.619 <2e-16 ***
## bs(horsepower, df = 7)1 5.3774 6.3634 0.845 0.4009
## bs(horsepower, df = 7)2 5.5408 3.1509 1.758 0.0830 .
## bs(horsepower, df = 7)3 -0.8369 3.6078 -0.232 0.8172
## bs(horsepower, df = 7)4 -3.8432 3.7673 -1.020 0.3111
## bs(horsepower, df = 7)5 -9.0415 4.6610 -1.940 0.0564 .
## bs(horsepower, df = 7)6 -15.0646 6.1037 -2.468 0.0160 *
## bs(horsepower, df = 7)7 0.2628 4.6155 0.057 0.9548
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.053 on 71 degrees of freedom
## Multiple R-squared: 0.5968, Adjusted R-squared: 0.557
## F-statistic: 15.01 on 7 and 71 DF, p-value: 7.304e-12#natural spline
mod_spline3nsjp = lm(mpg ~ ns(horsepower, df=7),data=Data.Jepang)
summary(mod_spline3nsjp)##
## Call:
## lm(formula = mpg ~ ns(horsepower, df = 7), data = Data.Jepang)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.2418 -2.2803 -0.4509 1.4414 10.7903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 32.119 2.162 14.855 < 2e-16 ***
## ns(horsepower, df = 7)1 1.731 2.489 0.696 0.489
## ns(horsepower, df = 7)2 2.796 3.800 0.736 0.464
## ns(horsepower, df = 7)3 -6.496 5.675 -1.145 0.256
## ns(horsepower, df = 7)4 -2.724 3.135 -0.869 0.388
## ns(horsepower, df = 7)5 -17.361 3.665 -4.737 1.08e-05 ***
## ns(horsepower, df = 7)6 -1.522 5.697 -0.267 0.790
## ns(horsepower, df = 7)7 -5.532 3.667 -1.509 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.08 on 71 degrees of freedom
## Multiple R-squared: 0.5915, Adjusted R-squared: 0.5512
## F-statistic: 14.69 on 7 and 71 DF, p-value: 1.134e-11#Buat plot bs
plot.splinebsjp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Basic Spline")+
stat_smooth(method = "lm",
formula = y~bs(x, df=7 ),
lty = 1,se=F)
plot.splinebsjp 
#Buat plot ns
plot.splinensjp<-ggplot(Data.Jepang,aes(x=horsepower, y=mpg)) +
geom_point(alpha=0.55, color="black") + labs(title= "Natural Spline")+
stat_smooth(method = "lm",
formula = y~ns(x, df=7),
lty = 1,se=F)
plot.splinensjp
#Secara empiris Ns-Spline
set.seed(123)
cross_valnsjp <- vfold_cv(Data.Jepang,v=5,strata = "mpg")## Warning: The number of observations in each quantile is below the recommended
## threshold of 20. Stratification will be done with 3 breaks instead.#Mencoba dengan jumlah basis 5-8, jumlah basis 6 maka knots = 3 (dengan rumus basis = knots +3)
df <- 5:15
best_spline3nsjp <- map_dfr(df, function(i){
metric_spline3nsjp <- map_dfr(cross_valnsjp$splits,
function(x){
modnsjp <- lm(mpg ~ ns(horsepower,df=i),
data=Data.Jepang[x$in_id,])
prednsjp <- predict(modnsjp,
newdata=Data.Jepang[-x$in_id,])
truth <- Data.Jepang[-x$in_id,]$mpg
rmse <- mlr3measures::rmse(truth = truth,
response = prednsjp
)
mae <- mlr3measures::mae(truth = truth,
response = prednsjp
)
metricnsjp <- c(rmse,mae)
names(metricnsjp) <- c("rmse","mae")
return(metricnsjp)
}
)
metric_spline3nsjp
# menghitung rata-rata untuk 5 folds
mean_metric_spline3nsjp <- colMeans(metric_spline3nsjp)
mean_metric_spline3nsjp
}
)
best_spline3nsjp <- cbind(df=df,best_spline3nsjp)
# menampilkan hasil all breaks
best_spline3nsjp#berdasarkan rmse
best_spline3nsjp %>% slice_min(rmse)#berdasarkan mae
best_spline3nsjp %>% slice_min(mae)4. Komparasi Model Jepang
Setelah didapatkan hasil pendekatan dengan berbagai macam pendekatan, berikut ini adalah hasil komparasi terhadap nilai rmse dan mae nya untuk kemudian ditentukan mana pendekatan yang paling baik untuk mengetahui pola hubungan antara variabel mpg dan horsepower.
library(ggplot2)
library(tidyverse)
library(splines)
library(rsample)
library(gridExtra)
#Perbandingan semua model
nilai_metric <- rbind(mean_metric_linear1jp,
mean_metric_poly_d2jp,
mean_metric_poly_d3jp,
mean_metric_poly_d4jp,
best_tangga1jp %>% select(-1) %>% slice_min(rmse),
best_spline3nsjp %>% select(-1)%>% slice_min(rmse)
)
nama_model <- c("Linear",
"Poly2",
"Poly3","Poly4","Tangga(breaks=5)",
"Natural Spline(df7)"
)
data.frame(nama_model,nilai_metric)#Grafik Gabungan komparasi
grid.arrange(plot.linierjp, plot.polinomd2jp, plot.polinomd3jp, plot.polinomd4jp, plot.tanggajp, plot.splinensjp)
Dari output tersebut dapat disimpulkan bahwa regresi polinomial derajat 3 adalah pendekatan yang paling baik dalam menjelaskan hubungan antara variabel mpg dan horsepower, hal itu terlihat dari nilai rmse nya paling kecil.
I. Index Horsepower/mpg
Untuk menentukan mobil yang diproduksi dari wilayah mana yang memiliki index tertinggi, berikut ini disajikan dengan analisa anova
library(ISLR)
library(dplyr)
library(car)
library(agricolae) # paket untuk uji lanjut
library(ggpubr) # paket menyusun grafik
attach(Auto)## The following objects are masked from Auto (pos = 15):
##
## acceleration, cylinders, displacement, horsepower, mpg, name,
## origin, weight, year## The following object is masked from package:ggplot2:
##
## mpg#Memanggil data yang akan dismulasikan
Data.Amerika<-Auto %>% filter(origin==1)
Data.Eropa<-Auto %>% filter(origin==2)
Data.Jepang<-Auto %>% filter(origin==3)
Data.komparasi<- Auto%>% mutate("HP/mpg" = horsepower/Auto$mpg)
Data.Anova<-Data.komparasi%>% select ("HP/mpg")
Perlakuan<-Data.komparasi%>% select ("origin")
Data.test.anova<- cbind(Perlakuan, Data.Anova)
# lm: kode instruksi untuk penyusunan model linier dalam analisa data
model = lm(Data.test.anova$`HP/mpg` ~ Data.test.anova$origin, data = Data.test.anova)
# anova: kode instruksi untuk analisa ragam satu jalur
anova(model) # HSD.test: kode instruksi untuk uji lanjut - Tukey HSD)
(HSD.test(model, "Data.test.anova$origin", unbalanced = FALSE)) ## $statistics
## MSerror Df Mean CV
## 12.50178 390 5.542063 63.79909
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Data.test.anova$origin 3 3.327215 0.05
##
## $means
## Data.test.anova$`HP/mpg` std r Min Max Q25 Q50
## 1 7.065061 4.274042 245 1.641026 21.500000 3.739130 5.434783
## 2 3.209041 1.507387 68 1.083521 8.209877 2.240152 2.843488
## 3 2.827014 1.168433 79 1.394850 6.100000 1.938839 2.314815
## Q75
## 1 10.000000
## 2 3.902174
## 3 3.879167
##
## $comparison
## NULL
##
## $groups
## Data.test.anova$`HP/mpg` groups
## 1 7.065061 a
## 2 3.209041 b
## 3 2.827014 b
##
## attr(,"class")
## [1] "group"Dapat disimpulkan bahwa mobil yang paling baik adalah buatan Amerika, hal itu terlihat dari hasil uji Tukey yang menunjukkan bahwa mobil yang berasal dari wilayah Amerika berbeda nyata dengan wilayah lain, dengan nilai mean 7,06.